PREVAIL: Predicting Recovery through Estimation and Visualization of Active and Incident Lesions.

OBJECTIVE: The goal of this study was to develop a model that integrates imaging and clinical information observed at lesion incidence for predicting the recovery of white matter lesions in multiple sclerosis (MS) patients.
METHODS: Demographic, clinical, and magnetic resonance imaging (MRI) data were obtained from 60 subjects with MS as part of a natural history study at the National Institute of Neurological Disorders and Stroke. A total of 401 lesions met the inclusion criteria and were used in the study. Imaging features were extracted from the intensity-normalized T1-weighted (T1w) and T2-weighted sequences as well as magnetization transfer ratio (MTR) sequence acquired at lesion incidence. T1w and MTR signatures were also extracted from images acquired one-year post-incidence. Imaging features were integrated with clinical and demographic data observed at lesion incidence to create statistical prediction models for long-term damage within the lesion. VALIDATION: The performance of the T1w and MTR predictions was assessed in two ways: first, the predictive accuracy was measured quantitatively using leave-one-lesion-out cross-validated (CV) mean-squared predictive error. Then, to assess the prediction performance from the perspective of expert clinicians, three board-certified MS clinicians were asked to individually score how similar the CV model-predicted one-year appearance was to the true one-year appearance for a random sample of 100 lesions.
RESULTS: The cross-validated root-mean-square predictive error was 0.95 for normalized T1w and 0.064 for MTR, compared to the estimated measurement errors of 0.48 and 0.078 respectively. The three expert raters agreed that T1w and MTR predictions closely resembled the true one-year follow-up appearance of the lesions in both degree and pattern of recovery within lesions.
CONCLUSION: This study demonstrates that by using only information from a single visit at incidence, we can predict how a new lesion will recover using relatively simple statistical techniques. The potential to visualize the likely course of recovery has implications for clinical decision-making, as well as trial enrichment.

Introduction

Multiple sclerosis (MS) is an inflammatory disease of the central nervous system, which is typically characterized by demyelinating lesions that occur in the brain and spinal cord. These lesions evolve dynamically from actively inflamed tissue over a period of months to more stable demyelinated regions of acute long-term axonal injury (Lassmann, 2013, Lassmann et al., 2007). A competing process of remyelination is also known to occur to varying degrees in patients, and has been documented in both relapsing-remitting and progressive cases (Patrikios et al., 2006, Bramow et al., 2010). Both the destructive and remyelinating processes are known to progress through the disease course (Frischer et al., 2015), and are associated with disability and morbidity. As therapeutics designed to promote tissue repair and remyelination are being developed, sensitive markers for in vivo assessment of these processes are increasingly important for studying therapeutic efficacy and patient management.

Magnetic resonance imaging (MRI) is a commonly used technique for identifying lesions, particularly in the white matter of the brain (Radü & Sahraian, 2008). The presence of new and active lesions is a key factor in the diagnosis and monitoring of MS, and several MRI sequences have been demonstrated to be effective in measuring the severity of these lesions (Polman et al., 2011, Sweeney et al., 2016, Sweeney et al., 2013, Pike et al., 2000). In recent years, successful attempts have been made to utilize quantitative methods in concert with MRI for the study of tissue damage in lesions. These techniques have included the use of advanced quantitative MRI sequences including T1 mapping (Larsson et al., 1989, Vrenken et al., 2006), magnetization transfer imaging (van Waesberghe et al., 1998, van Waesberghe et al., 1999), and diffusion tensor imaging (Narayanan et al., 1997, Werring et al., 1999, Filippi et al., 2001), as well as statistical techniques for modeling tissue damage using conventional MRI (Shinohara et al., 2011, Mejia et al., 2015, Reich et al., 2015) and the development of time-series models to examine lesion activity (Sweeney et al., 2016, Meier et al., 2007, Meier and Guttmann, 2003, Meier and Guttmann, 2006).

Specifically, much research has engaged with the apparent paradox related to the lack of coherence between the presence of lesions and clinical disease measures (Barkhof, 2002). One recent study retrospectively related the longitudinal behavior of lesions, as opposed to simply their presence, to clinical covariates and treatment status (Sweeney et al., 2016). Significant relationships between treatment and longitudinal behavior indicated that receiving disease-modifying therapy or steroids was associated with a better healing trajectory within lesion tissue. These findings signify the presence of potentially important relationships between the repair processes in the brain, therapeutics, and disability.

Unfortunately, today there is still relatively little that can be determined in advance about the way specific lesions will recover, or the degree to which they may be responsive to treatment. The ability to visually examine the likely course of recovery for a given incident lesion would have the potential to be useful in several settings. Specifically, such visualizations could be a beneficial tool for physicians, providing them important supplemental information when making treatment decisions. Additionally, knowledge of how patients' brains are likely to recover from lesion damage could be beneficial in clinical trials, for which advanced knowledge of lesion characteristics could inform recruitment enrichment and trial design.

To build on the previous work, and to address the needs outlined above, the current study attempted to develop a statistical model that would be capable of prospectively predicting how lesions would heal over the course of a year. In this paper, we discuss the development of such prediction models for two outcome MRI modalities, we present statistical and clinical measures of validity and prediction accuracy, and we discuss the implications and potential next steps of this line of research.

Methods

Image acquisition and preprocessing

Details of the image acquisition and preprocessing have been previously published (Sweeney et al., 2016) and are summarized in this section. Whole-brain two-dimensional T2-weighted FLAIR, PD, T2, and three-dimensional T1-weighted volumes were acquired in a 1.5 tesla (T) MRI scanner (Signa Excite HDxt; GE Healthcare, Milwaukee, Wisconsin) using the body coil for transmission. The 2D FLAIR, PD, and T2 volumes were acquired using fast-spin-echo sequences, and the 3D T1 volume was acquired using a gradient-echo sequence. All scanning parameters were clinically optimized for each acquired image.

For image preprocessing, we used Medical Image Processing Analysis and Visualization (http://mipav.cit.nih.gov) and the Java Image Science Toolkit (http://www.nitrc.org/projects/jist) (Lucas et al., 2010). All images for each subject at each visit were interpolated to a voxel size of 1 mm3 and rigidly co-registered longitudinally and across sequences to a template space (Fonov et al., 2011). To coregister the T1 images across study visits, a two-step procedure was applied: first, subject-specific templates were generated by averaging after rigid alignment of the T1 images to the MNI template. Second, all T1 images were then realigned to the subject-specific templates. Finally, the additional MRI sequences were aligned to the T1 images within each study visit and this transformation was composed with the T1-based transformation to the subject-specific template.

Extracerebral voxels were removed using a skull-stripping procedure (Carass et al., 2007) and the brain was automatically segmented using the T1 and FLAIR images (Shiee et al., 2010) to produce a mask of normal-appearing white matter (NAWM), or white matter excluding lesions. Intensity normalization was then conducted using z-scoring based on the mean and variance of the variability in the NAWM (Shinohara et al., 2011, Shinohara et al., 2014). After preprocessing, studies were manually quality controlled by a researcher with over five years' experience with structural MRI (EMS) and studies with motion or other artifacts were removed.

Patient demographics

For this study, 60 subjects diagnosed with MS were scanned between 2000 and 2008 on a monthly basis over a period of up to 5.5 years (mean = 2.2 years, sd = 1.2) as part of a natural history study at the National Institute of Neurological Disorders and Stroke in Bethesda, Maryland. To be included in the analysis, subjects were required to meet certain pre-specified inclusion criteria. Specifically, only subjects with at least one new lesion during the observation period were included, and these subjects were required to have been rescanned at least twice 360 days after lesion incidence. 32 subjects met these criteria and were included in the analyses. The 32 subjects ranged from 18 to 60 years of age, with a mean age of 37 years (sd = 9). Of the 32 subjects, 11 were male and 21 were female. The majority of the subjects (n = 27) were diagnosed with relapsing-remitting MS, and the remaining five were characterized as secondary-progressive. Subjects were either untreated or treated with a variety of disease-modifying therapies during the observation period, including both FDA-approved therapies (Avonex, Betaseron, Daclizumab, and Rebif) and experimental therapies.

Prediction model

Outcomes

The outcomes of interest in this study were 1) normalized T1-weighted voxel intensity (nT1w) (Shinohara et al., 2014) and 2) MTR voxel intensity approximately one-year post-incidence, and is denoted by Y(v) for subject i in voxel v. Due to the noise inherent in both sequences, outcome variables were created by averaging the intensity of each voxel at the visit immediately following the 360-day cutoff (referred to as the one-year visit), the visit prior to the one-year visit (mean = 10.6 months from incidence, sd = 1.3 months), and the visit following the one-year visit. Because no change is expected in the lesion after that length of time, this average only reduced variability due to measurement error (Meier et al., 2007). Thus, the average score represents a more precise estimate of true voxel intensity than the one-year visit intensity alone.

Predictors

A dataset made up of scan data and relevant demographic variables was created to predict the one-year post-incidence voxel intensities. For each voxel, this included the MTR as well as the nFLAIR, nPD, nT2w, and pre- and post-contrast nT1w intensities at incidence, denoted by (v). After applying a 3D Gaussian smoother with variance parameter 3 mm and width 5 mm, each voxel's blurred intensities on the five scan modalities, , were also included, as well as the distance, in number of voxels, from the voxel to the nearest boundary of the lesion, d(v), and the size, in number of voxels, of the lesion, s(v). Additional predictors included were the patient's age, sex, disease subtype, expanded disability status score (EDSS; (Kurtzke, 1983)), disease-modifying treatment status (treated versus untreated, with use of one or more therapies counting as treated), and steroid status (receiving steroids versus not on steroids) at the time of lesion incidence.

Model creation and validation

All statistical modeling was conducted in the R statistical environment (R: a Language and Environment for Statistical Computing, 2015). Separate linear regression models for nT1w intensity one-year post-incidence and MTR intensity one-year post-incidence were created using all of the variables in the dataset, as well as the interaction between voxel distance to lesion boundary and lesion size:

Predictions were obtained using leave-one-lesion-out cross-validation. In this cross-validation technique, the prediction for each voxel in a given lesion is made using a model trained on all of the data except those from that lesion. As a result, 397 models were developed, each excluding one of the 397 lesions in the dataset. This method ensures that the prediction for each voxel is not influenced by the true outcome of that voxel. Secondary cross-validation was also performed using a leave-one-subject-out technique, but due to the small sample size and large number of variables in the model, performance was assessed on the predictions obtained by the leave-one-lesion-out procedure.

Performance assessment

The performance of the model was assessed in two ways. First, prediction accuracy was measured quantitatively by calculating the root mean square error (RMSE) of both the nT1w and MTR predicted intensities for estimating the average one-year intensity outcomes. This measure gives an estimate of the average difference between the predicted intensity and observed intensity of the voxels in the dataset. Because the RMSE is dependent on the scale of the outcome, the measurement error of voxel intensity was estimated for comparison. The measurement errors for T1w and MTR were estimated by calculating the RMSE of the voxel intensities at the pre-one-year visit for estimating the average one-year intensity outcomes. This directly compares the accuracy of the model's prediction using information at incidence to the accuracy of a scan taken at approximately 11 months (mean = 10.6, sd = 1.3 months), for predicting the average intensity of a voxel at one year.

The second method for assessing accuracy was a rater study. Three board certified MS clinicians (two neurologists and one neuroradiologist) with between 5 and 12 years of research experience in MS participated in this validation. To assess the prediction performance visually, the raters were asked to individually score how similar the model-predicted one-year appearance was to the observed one-year appearance for a random sample of 100 lesions. Raters viewed images of the lesion voxels' intensities at incidence, their predicted intensities one-year post-incidence, and their observed intensities one-year post-incidence.

For each lesion in the sample, raters were asked to determine “overall how well the prediction reflects the appearance of the lesion after one year,” “how well the degree of recovery in the prediction reflects the degree of recovery after one year,” and “how well the pattern of recovery in the prediction reflects the pattern of recovery after one year.” Each question was asked for both the nT1w and MTR predictions, resulting in six ratings per lesion. Ratings were given on a 1-to-4 scale, with labels of “1 - Failed miserably,” “2 - Some redeeming features,” “3 - Passed with minor errors,” and “4 - Passed.” Raters were broadly instructed that “Failed miserably” indicated no correspondence whatsoever between prediction and observed images, “Some redeeming features” indicated some correspondence, “Passed with minor errors” indicated correspondence, and “Passed” indicated excellent correspondence. Images were observed privately and ratings were given independently, with no discussion by raters occurring during the rating process.

Results

RMSE

Using this measure, both cross-validated models performed well. The overall RMSE of the T1w prediction was 0.95 (mean of lesion RMSEs = 0.89, sd = 0.35), as compared to a measurement error in T1w scans of 0.48. This indicates that the correspondence between the predicted intensity and the observed intensity of each voxel was only slightly (approximately one-half standard deviation unit) lower than the correspondence between the voxel intensity of a scan taken at approximately 11-months and the true 12-month observed intensity. The RMSE of the MTR prediction was 0.064 (mean of lesion RMSEs = 0.057, sd = 0.023), as compared to a measurement error in MTR scans of 0.078, where the average MTR outcome value was.36. This demonstrates that the correspondence between the model-predicted intensity and the observed intensity of each voxel was better than the correspondence between the voxel intensities approximately one month apart. This is consistent with the literature, as MTR has been demonstrated to be more noisy than T1w (Reich et al., 2015). Prediction images demonstrating above-average, average, and below-average performance based on the RMSE measure are presented in Fig. 1, Fig. 2.

Fig. 1

Representative predictions for the nT1w model. Images are axial slices of lesions, with rows representing three example lesions with varying levels of predictive accuracy. For nT1w intensities, red areas represent hyperintensity and blue areas represent hypointensity, with 0 (white) representing the intensity of normal-appearing white matter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Representative predictions for the MTR model. Images are axial slices of lesions, with rows representing three example lesions with varying levels of predictive accuracy. For MTR, intensities range from 0 to 1, with a mean of approximately 0.5 for normal-appearing white matter.

Rater study

Consistency between raters

Of the 100 lesions included in the rater study, 4 were excluded for segmentation errors. Six ratings were obtained for each image, resulting in 576 distinct scores for each rater. In 17% (n = 98) of these cases, all three raters assigned the same score. In another 54% (n = 309), two raters assigned the same score and the third gave a score either one unit above or one unit below the score given by the other two. Thus, in 71% of the ratings, the raters all assigned scores within one unit of each other.

Ratings

Scores were averaged across the three raters to obtain mean ratings for each image. Using these mean ratings, both the nT1w and MTR prediction models performed well. For both models, the median score for the overall similarity between the predicted image and the true lesion image was 3.0, corresponding to a rating of “Passed with minor errors.” For the specific rating of similarity of degree of recovery between the predicted and observed images, the median for both nT1w and MTR was 2.6. This score was partway between the rating “Some redeeming features” and the rating “Passed with minor errors.” For the similarity of “pattern of recovery” between the predicted and observed images, the median for both nT1w and MTR was 3.00 (Fig. 3). Overall ratings of accuracy were significantly (p < 0.001) negatively associated with a lesion's RMSE (r = − 0.60, r = − 0.56). Additionally, ratings of degree were more associated with RMSE (r = − 0.71, r = − 0.68) than ratings of pattern (r = − 0.52, r = − 0.53), but all associations were statistically significant (p < 0.001).

Fig. 3

Scores for the six rater study questions, averaged across the three raters. The three rows show the distributions of the ratings of overall accuracy, accuracy of the degree of healing, and accuracy of the pattern of healing, respectively. Plots in the first column are distributions of ratings of the nT1w prediction images, and plots in the second column are distributions of ratings of the MTR prediction images.

Discussion

In this paper, we developed models to estimate the appearance of white matter MS lesions one year after their incidence. Both the model for normalized T1w voxel intensity and the model for MTR voxel intensity produced accurate, cross-validated predictions of intensity at one-year post-incidence. This accuracy was measured quantitatively using RMSE for estimating the true one-year intensity with the statistically predicted one-year intensity. RMSE for the predicted intensity was compared to the measurement error of voxel intensity, operationalized as the RMSE for estimating intensity at a one-year post-incidence visit with the intensity at the visit immediately preceding it. In the model for nT1w intensity, the prediction accuracy was approximately one-half standard deviation unit larger than the value of measurement error, indicating that the nT1w model was able to predict one-year intensities almost as accurately as these can be assessed directly by re-acquiring the image. For the model of MTR intensity, the prediction accuracy surpassed that of measurement error. The accuracy of the model was also confirmed by three board-certified MS clinicians, who rated the correspondence between images produced by the model and images of true one-year intensities.

During the development of the models, flexible spline modeling was explored to account for possible nonlinear relationships. This method did not improve the predictions, and thus the simpler linear regression approach was used for the final models. It is also worth noting that voxels containing vasogenic edema, as described previously (Sweeney et al., 2016), were included in the data for both model development and accuracy testing. Though edema behaves differently than structurally damaged lesion tissue, and its healing pattern is potentially easier to predict, it was determined that for the purpose of these predictions it was appropriate to leave edema voxels in the data. This is largely due to the fact that in many cases, differentiating between edema and lesion is difficult without examining how the tissue changes over the course of several weeks (Sweeney et al., 2016). Since our model was developed for the purpose of predicting healing using only information obtained at the time of lesion incidence, using data from later visits to categorize and remove edema would not accurately represent the way the model would be used clinically. However, voxels were categorized as edema or lesion after the models had been developed for performance assessment, and we found that the prediction accuracy was comparable between edema voxels and lesion voxels for both nT1w (RMSE = 0.80, RMSE = 1.04) and MTR models (RMSE = 0.061, RMSE = 0.066).

There were some limitations to the predictions developed in this study. Spatial covariance was not explicitly modeled, though spatial relationships were loosely accounted for by including the voxel intensities of spatially smoothed images. Additionally, the measure of voxel location used (distance to the boundary) likely could be improved upon as well, as it is unable to account for differing lesion shapes. However, in spite of the relative simplicity of the spatial and locational variables, the models were still able to achieve good accuracy in these areas. Both T1w and MTR predictions had a median rating of “Passed with minor errors” with respect to the pattern of healing.

An additional limitation arises in the interpretation of regression coefficients. For both T1w and MTR models, treatment with steroids was associated with lower voxel intensity at the one-year follow-up, indicating worse healing in patients on steroids. This is most likely an effect of the observational nature of this study, which suggests that the steroids variable is capturing an aspect of disease severity or activity that was not accounted for completely by EDSS. Including other characteristics in the model could take steps to address this phenomenon, but ideally future work would test this model in the context of clinical trial data in order to obtain a clearer sense of how treatment impacts the recovery of incident lesions.

A strength of the current study is the novel integration of neuroimaging with demographic and clinical data to predict how lesion tissue will heal in MS patients. We believe that accurate predictions of this sort may have several important applications in MS treatment and research. In MS research, clinical trials may benefit from the ability to predict which patients are more likely to be responsive to treatment. Clinically, this tool could enable physicians to view a prediction of the likely course of healing of patients' incident lesions in order to make more informed and personalized treatment decisions. As such, future research will focus on the refinement of the treatment and steroid effect estimation in the model through the use of clinical trial data, with the goal of facilitating the direct comparison of predicted recovery when treated with various disease modifying therapies, steroids, or when left untreated. This would provide clinicians and researchers with previously unavailable information about the course a lesion is likely to take, and would allow for greater personalization of treatment decisions, as well as better informed and more powerful study designs.