Using quantitative magnetic resonance imaging to track cerebral damage in multiple sclerosis: a longitudinal study

2.1 Population

Seventeen patients, recruited at the specialized MS outpatient clinic of the CHU Liège, Belgium, with a diagnosis of MS according to the McDonald criteria 2010 [27], were gathered from two studies: ten of them were part of the work reported by Lommers et al. 2019 [24], the other seven were recruited from another MS study taking place at the GIGA Cyclotron Research Centre – In Vivo Imaging (Liège, Belgium) [28]. For the first study (10 subjects), the inclusion criteria were (1) age between 18 and 65 years ; (2) Expanded Disability Status Scale (EDSS) inferior or equal to 6.5 ; (3) absence of relapse within the previous four weeks ; Both RRMS and PMS patients were recruited. The second study (7 subjects) differs a bit as it comprises only RRMS patients, and the inclusion criteria were (1) age between 18 and 45, (2) EDSS between 0 and 4, (3) absence of relapse for at least 6 months prior to the study, (4) disease duration was below or equal to 5 years. For both studies, compatibility with MRI and absence of other neurological/psychiatric diseases were required. These studies were approved by the local ethics committee (approval numbers B707201213806 and B707201835630, respectively). All patients were followed up and scanned twice on the same 3T MRI scanner, every 1 to 3 years. For each of the 17 MS patients, data from two MRI sessions were available, at T0 and T1. This cohort included 11 RRMS and 6 (primary and secondary) PMS patients. Thirteen were receiving disease-modifying treatments (DMTs). The patients’ median age was 36 years (range: 25-65) and the median time interval between two scans was 30 months (range: 14-61). Demographic data appears in Table 1. Extended individual information appears in Supplementary data.

2.2 MR image acquisition

MRI data were acquired on a 3T whole-body MRI-scanner (Magnetom Prisma, Siemens Medical Solutions, Erlangen, Germany). The whole-brain MRI acquisitions included a multi-parameter mapping protocol (MPM), from which one can simultaneously estimate (semi)quantitative maps of magnetization transfer saturation (MTsat), proton density (PD), transverse relaxation (R1) and effective longitudinal relaxation (R2*). This protocol arising from an international collaborative effort [10, 26], has already been used to study brain microstructure in various conditions including normal aging [26, 9, 30], brain tumor [31], Parkinson’s disease [32-34] as well as MS. It consists of three co-localized 3D multi-echo fast low angle shot (FLASH) acquisitions at 1mm_³ resolution and two additional calibration sequences to correct for inhomogeneities in the RF transmit field [35, 36]. The FLASH datasets were acquired with predominantly PD, T1 and MT weighting, referred to in the following as PDw, T1w and MTw, at multiple echo times. All three had high bandwidth to minimize off-resonance and chemical shift artifacts. Volumes were acquired in 176 sagittal slices using a 256×224 voxel matrix. GRAPPA parallel imaging was combined with partial Fourier acquisition to speed up acquisition time to approximately 20 min. An additional FLAIR sequence was recorded with spatial resolution 1mm_³ and TR/TE/TI=5000ms/516ms/1800ms. Extra B1 field mapping images (transmit B1+ and receive B1-fields) were also acquired to reduce spatial inhomogeneities related to B1 effect. This was essential for proper quantification of T1 (or R1=1/T1) in particular. Finally, B0 field mapping images, corresponding to both magnitude images and pre-subtracted phase image, were acquired for image distortions corrections. A summary of the acquisition parameters appears in Supplementary data.

2.3 MR image processing

All data processing was performed in Matlab (The MathWorks Inc., Natick, MA, USA) using SPM12 (www.fil.ion.ucl.ac.uk/spm) and three additional dedicated SPM extensions: the Lesion Segmentation Tool (LST) version 1.2.3 (www.statisticalmodelling.de/lst.html) [37], the “quantitative MRI and in vivo histology using MRI” toolbox (hMRI, http://hmri.info) [10], and “US-with-Lesion” tool (USwL, https://github.com/CyclotronResearchCentre/USwLesion). Quantitative maps - MTsat, PD, R1 and R2* - were estimated using the hMRI toolbox. T1w, PDw and MTw images acquired at multiple echo times (TE) were extrapolated to TE=0 to increase signal-to-noise ratio and remove the otherwise remaining R2* bias [10, 24, 38]. The TE=0 extrapolated MTw, PDw and T1w images were used to calculate MT saturation, R1 and apparent signal amplitude A* maps. PD map was derived from A* map, which is proportional to proton density. All quantitative maps were corrected for inhomogeneities from local RF transmit field (B1+), and R1 quantitative maps were further corrected for imperfect RF spoiling using the strategy of Preibisch and Deichmann [39]. The receive bias field map (B1-) was used to correct PD maps for instrumental biases. The R2* map was estimated from all three multi-echo series (MTw, PDw and R1w) using the ESTATICS model [38].

After generating quantitative maps using the hMRI toolbox for all sessions, spatial preprocessing involved the following steps (Figure 1): within-patient registration brought the two serial MR data sets into the individual T0 space, using the longitudinal registration tool from SPM [40]. For each individual patient, a preliminary WM lesion mask was generated based on FLAIR and T1w images by the lesion growth algorithm implemented in the LST toolbox [37], followed by manual corrections by an MS expert (EL) to remove aberrant/artefactual lesion detections [24]. The images were then segmented using the USwL toolbox, which consists of an extended version of the traditional Unified Segmentation (US) algorithm [41] and includes an additional tissue class representing the WM lesion(s). The USwL method internally generates a subject-specific extended set of tissue probability maps (TPM) [42]: an extra tissue class, based on the smoothed preliminary lesion mask warped into template space (using cost function masking during normalization [43]), is added to account for the lesion, and the original white matter prior map is updated accordingly [44]. The grey matter TPM was not updated due to a very low number of lesions present in the cortical ribbon. Multi-channel segmentation was conducted, using MTsat, PD, R1 and FLAIR images. This pipeline did not use the PD-, T1- and MT-weighted images acquired for the MPM maps construction, but the parametric maps themselves instead. In this way, voxels do not depict MR intensities but rather physical quantitative parameters. The method generated the segmented tissue classes (a posteriori tissue, including lesion, probability maps), as well as spatial warping into standard template space. The preliminary lesion mask was used as input for the first session data (at T0) then the a posteriori lesion map generated at this initial step served as prior to the subsequent session (at T1). Segmentation teased out the different tissue classes of interest: NAWM, NACGM and NADGM, as well as WM lesions. To analyze the microstructure within those tissue classes, a posteriori tissue maps were binarized and tissue-specific independent masks were constructed: each voxel is assigned to one single tissue class with the highest probability for that voxel (provided that this probability was above 0.2). The lesion binary mask was further cleaned for lesions <10mm_³ which likely resulted from segmentation errors. Finally, binarized tissue class masks were in turn applied on the MPM maps to extract voxel values inside them.

• Download figure
• Open in new tab

Figure 1:

Chartflow of data creation and processing (see text). MPM maps were created with the hMRI-toolbox, FLAIR images were directly acquired for both sessions (T0 and T1). A preliminary mask was constructed based on T0 FLAIR. All images (MPM and FLAIR, T0 and T1) were co-registered to the MPM T0 space. Segmentation using USwL allowed to isolate the different tissue classes.

2.4 Brain volume change

Volumetric changes were investigated using the USwL a posteriori tissue probability maps. The following measures of brain volume were computed for each session of each participant: (1) Total intra-cranial volume (TIV) = volume (NAWM + GM + CSF + lesions), (2) brain parenchymal fraction (BPF) = volume (NAWM + GM + lesions)/TIV, (3) GM fraction (GMF) = volume (GM)/TIV, and (4) lesion fraction (LF) = volume (lesion)/TIV. The percentage of change between both scanning sessions was evaluated for each volumetric measurement, then annualized changes were computed by dividing these measures by scan intervals (in years). Results were directly analysed with a t-test (testing if significantly different from 0 at p < .05), but also in the same way as the normal appearing tissues MR parameters in relation to the patients’ clinical status (see next section).

2.5 Analysis of normal appearing tissues

The median value of quantitative MRI parameters was extracted from the three normal appearing tissues (NAWM, NACGM and NADGM), and an individual annual rate of change (ARoC) was computed for each parameter in each tissue class, based on the initial and final values and accounting for the time interval (in years) between scans. This rate of change in qMRI parameters served as dependent variable in a general linear model testing the effect of clinical status: Y is the ARoC for a qMRI parameter and tissue class, β’s are the regression parameters corresponding to the associated regressor (with β₀ the intercept), and ε the residuals. X_status is a binary categorical variable representing the patient’s disease activity status: a status score of 1 was assigned to patients stable or improving from T0 to T1.

This patient status X_status was derived from one score of disease activity: NEDA-3 (No Evidence of Disease Activity [45]), a composite of three related measures of disease activity. A score of 0 was assigned in the presence of new clinical relapses (only concerning RRMS patients) and/or MRI activity (new or enlarged lesions visible on FLAIR T2 or Gadolinium-enhanced images) and/or sustained disability progression over six months based on Expanded Disability Status Scale (EDSS). For both RRMS and PMS patients, disability progression was defined as a 1.0-point increase if the EDSS score was ≤ 4.0 at baseline and as a 0.5 point increased if the baseline EDSS score was > 4.0. The threshold of 4.0 was proposed in this study because it is considered as a milestone regarding ambulatory performance.

NEDA-3 and was evaluated at mid- and end-scanning interval, and a final status score of 0 was given only to patients for which disease activity or progression was noted in both cases, indicating a clear progression of the disease over the whole interscan interval.

The influence of several clinical measurements such as 25 FWT, 9HPT and SDMT was also considered to refine the evaluation of disease activity. However complete data were lacking for several patients. Moreover, when available, these additional clinical parameters did not modify the final X_status. Longitudinal clinical information allowing to derive the disease activity status for each subject appears in Supplementary data.

Permutation tests were employed for inferences [46]. R-squared value was tested against computed statistics after permutation of the data. For a number n of permutations, the X_status. values were randomly shuffled (constructing a new regressor written ), tested against the unchanged response Y, and generating each time a permuted R-squared value (noted R, R_obs being the true R-squared value computed without permutation of the data). The condition is verified at each permutation. After n permutations (with n = 5000 in this study), a p-value was computed based on the following formula: which estimates the probability of obtaining R_obs under the null hypothesis that Y is not correlated to X_status. The null hypothesis is rejected if p < .05 FDR-corrected for multiple comparisons [47], for the 12 tests performed (3 tissue classes with 4 qMRI parameters).

Two-tailed t-tests were applied post-hoc on the significant results of permutation tests to compare the ARoC distribution between disease status, i.e., X_status = 0 against X_status = 1. Inferences were conducted at a significance level of .05.

The same pipeline was applied to the brain volumetric changes (BPF, GMF and LF) to test their correlation to the disease activity status.

2.6 Analysis of lesions and peripheral tissues

For white matter lesions analysis, we did not use ARoC but exploited directly the qMRI parameters voxel values. Importantly, with USwL, the prior lesion mask is only used in a probabilistic way and the estimated posterior lesion map, obtained using MTsat, PD, R1 and FLAIR images, typically showed more extended lesion than clinically visible on the FLAIR image alone.

Therefore, we separated focal lesions detected on FLAIR images, with LST and visual inspection, from their peripheral regions detected on qMRI maps. Two different peripheral regions were considered: one for each time point (T0 and T1). Therefore, at T0, three distinct lesion-related regions were isolated:

• The lesions, as clinically defined, corresponding to hyperintensity on conventional FLAIR MR image acquired at T0 (referred to as ‘focal FLAIR lesion’).

• The peripheral region detected on qMRI maps at T0, not including the FLAIR lesion (referred to as ‘initial peripheral lesion’).

• The peripheral region, detected on qMRI maps at follow up (appearing hyperintense on FLAIR images or not); this was computed by masking out the T1 lesion mask with the T0 lesion mask (referred to as ‘later peripheral lesion’). This region allows us to determine whether its microstructure at T0 forebodes a full-blown plaque, detectable during follow up.

The three areas were compared between each other and with NAWM, in order to characterize them on a microstructural basis (Figure 2). For an accurate lesion-by-lesion analysis, only enlarging lesions, i.e., present in the three masks, were considered for these comparisons.

• Download figure
• Open in new tab

Figure 2:

Schematic illustration of the NAWM and 3 lesions-related areas: focal FLAIR lesion (dark gray area), initial peripheral lesion detected at T0 (medium gray area), later peripheral lesion detected at T1 (dashed, left, and light gray, right, area).

NAWM region consisted of all white matter voxels which did not belong to any of the three lesion-related regions. The four areas are not overlapping as no voxel could belong to more than one class at the same time.

For all participants, MTsat, PD, R1 and R2* median values were extracted from each lesion area, considering lesions individually (between 2 and 66 measurements per subject). Similarly, the median qMRI values within NAWM were also extracted (one measurement per subject). These values were extracted from T0 and T1 scans separately. Statistical analyses were performed in SAS 9.4 (SAS Institute, Cary, NC). None of the qMRI parameter was normally distributed, therefore we applied a log transformation on each of them prior to statistical analysis. For each qMRI parameter, a separate Generalized Linear Mixed Model (GLMM) tested the effect of areas (NAWM and the three lesion-related areas), and time points (T0 and T1), as well as their interaction (i.e., area*time), on the median qMRI parameter value, with a first-order autoregressive variance/covariance model and participants as a random factor (intercept). The degrees of freedom were estimated using Kenward-Roger’s method. Statistical significance was estimated at p < .05 after adjustment for multiple comparison using Tukey’s procedure.

3.1 Volume changes

BPF annually decreased between T0 and T1 by -0.67 ± 1.12% (significantly different from zero; paired-sample t-tests; t 16) = 2.57; p = .0204) whereas LF increased by 22.88 ± 26.13% (t 16) = 3.70; p = .0019). GMF non-significantly decreased by -0.30 ± 1.44%.

3.2 Analysis of normal appearing tissues

As expected, changes in MTsat and R2* within NABT between T0 and T1 varied across subjects (Figure 3). PD and R1 exhibited similar behaviors, see Supplementary data. At the group level, with the GLM analysis and permutation inference, we observed that the ARoC of MTsat and R2* positively regressed with disease status as follows (Table 2): MTsat in NAWM and NACGM and R2* in NAWM significantly increased in patients who fare well (X_status = 1).

• Download figure
• Open in new tab

Figure 3:

Line plots illustrating individual ARoC’s for MTsat (left) and R2* (right) in NAWM. Each line corresponds to one subject. Dotted lines represent increasing rates.

Post-hoc t-tests applied on these significant results for a clearer illustration of the difference in disease status (Figure 4) were all significant at a level of .05.

• Download figure
• Open in new tab

Figure 4:

Violin plots of significant change rates in microstructure with respect to X_status. From left to right: MT in NAWM, MT in NACGM, R2* in NAWM. * P < .05.

Regarding BPF and LF, their correlation to the disease activity status was not significant (Table 2), suggesting that qMRI parameters are more sensitive to subtle microstructural changes in NABT over time than global morphological measurements

3.3 Analysis of lesion microstructure

The number of enlarging WM lesions between T0 and T1 varied from 2 to 66 across patients, for a total of 741 identified enlarging lesions among all subject, corresponding on average among patients to 63% (±31%) of the amount of initial focal lesions. The number of enlarging lesions did not significantly differ between patients’ disease status groups (t 15) = .244, p = .811). GLMMs found a significant effect of areas (3 lesional regions and NAWM) for MT, R1, R2* and PD median (MT:F₃ = 35.34, p < .0001, PD:F₃ = 68.03, p < .0001, R1:F₃ = 40.26, p < .0001, R2*:F₃ = 32.32, p < .0001). By contrast, neither time effect (T0 vs T1; MT:₃ = 0.36, p = .5481, PD:F₃ = 1.20, p = .2735, R1:F₃ = 2.05, p = .1520, R2*:F₃ = 2.86, p = .0911), nor the area*time interaction (MT: = 0.09, p = .9671, PD: = 0.14, p = .9346, R1:F₃ = 0.14, p= .9331, R2*:F₃ = 0.40, p = .7565) were significant, suggesting the microstructural stability of the initial lesion core. Post-hoc tests confirmed significant differences between the four tissue areas. At T0 and T1, MTsat, R1 and R2* values were significantly larger in the initial peripheral lesion than FLAIR lesion, in the later peripheral lesion than the initial one, and in the NAWM than later peripheral lesion. The reverse was observed for PD. The significant difference in parameters between initial and later peripheral lesion at T0 suggests that subtle microstructural changes appear in the periphery of the initial lesion, months before their detection as focal FLAIR lesions at T1. Adjusted p-values appear in Figure 5. Detailed statistical results of the GLMM’s appear in Supplementary data.

• Download figure
• Open in new tab

Figure 5:

Microstructural parameters in NAWM and the 3 lesion-related areas, for each scanning time T0 and T1. P-values were obtained with post-hoc tests on the tissue area effect. * P < .05.

Quantitative MRI parameter time course within NABT

In this study, we used a multiparameter mapping protocol that was gradually optimized and validated for multi-centric studies [48]. It provides high-resolution maps of multiple qMRI parameters from data acquired during a single scanning session of acceptable duration. A number of cross-sectional studies using a combination of MT, R1, R2* or PD parameters reported significant changes in the microstructure of NABT in MS [24, 49-55]. By contrast, longitudinal analyses of multiparameter qMRI data are scarce. A progressive shortening of T2/T2* [56] or increase in R2* [57-59] was reported within the basal ganglia, suggesting increased of myelin and/or iron contents as well as edema resorption. Likewise, PD and T1 increased within a year, suggesting a demyelination and/or axonal loss [60]. MTR progressively decreases in NAWM of MS patients over one [61] or two years [62]. These abnormalities tend to be more pronounced in progressive phenotypes [63] and were associated to a slow, diffuse and global myelin pathology.

Here, we showed that MTsat within NAWM and NACGM and R2* values within NAWM increase in clinically stable or improving patients. Because MTsat and R2* both correlate with myelin content [11, 29, 30, 64-66], our results suggest repair mechanisms within NABT of patients who are responding to disease modifying treatments, despite the initial myelin/axonal loss and independently from WM focal lesion evolution. Such increases could also be explained by an edema/inflammation resorption, but less likely than myelin/axonal density changes since MTsat is the least dependent to water content among the four qMRI parameters. These results echo cross-sectional analyses showing that healthy controls (HC) have higher MTsat and R2* values within the same tissue classes compared to MS patients [24]. Annual rates of change of R1 and PD within NABT were not significantly associated with the individual clinical status in this study, although R1 reduction within NABT has already been reported in cross sectional [24, 52, 53] and longitudinal [60] studies comparing MS subjects to HC.

Lesion microstructure

Focal inflammatory demyelinating lesions have been extensively characterized and are traditionally classified as active, chronic active (smoldering) or inactive plaques according to the presence and distribution of plaque-infiltrating macrophages/microglia [67-69]. Focal WM pathology is a constantly evolving process including episodes of demyelination and remyelination but also accumulation of irreversible axonal damage. Age, disease duration, clinical phenotype as well as disease modifying treatment all contribute to the dynamic nature of focal WM pathology [68, 70]. This accounts for the large inter- and intra-individual heterogeneity of MS, which conventional MRI is largely unable to capture. By contrast, quantitative MRI parameters are sensitive to myelin, axonal as well as iron contents and appear as promising markers of plaque dynamics. For instance, MTR was shown to sharply decrease within gadolinium enhancing lesions before recovering during the subsequent months [71-73]. Likewise, reduction of MTR within NAWM, days to weeks before the formation of a new active lesion was also demonstrated [74, 75], and long-term MTR changes in WM plaques were observed in relation with disease progression [63, 76]. The present study broadens the quantitative characterization of plaque dynamics, in keeping with previous longitudinal studies [56, 77]. Two important findings emerge from the results. First, qMRI refines lesion segmentation, as compared to the processing based on the sole FLAIR image. In consequence, the initial lesion revealed by qMRI is typically wider that the plaque detected in FLAIR. Its periphery is characterized by a decrease in MTsat and R2* as compared to NAWM, suggesting an incipient demyelination, reminiscent of the so-called ‘periplaques’ [78]. Moreover, MTsat, R2* and R1 values progressively decrease from NAWM to plaque core, suggesting a centripetal loss of myelin content. Second, plaque microstructure is altered in plaque periphery before any observable change in conventional MRI signals. This finding suggests, in keeping with neuropathological observations [68, 70, 79, 80] that subclinical ongoing inflammation and/or demyelination takes place in the periphery of an active plaque, well before it is detectable on FLAIR or T1 post-gadolinium sequences. If confirmed on larger population samples, this finding might significantly modify treatment management in MS patients.

Oddly enough, plaque qMRI parameters did not significantly change across time. Because iron concentration increases within chronic active or smoldering lesions [81, 82], we were expecting a progressive increase in R2* value. This negative result might be due to the small sample size, the short period of follow up or the limited sensitivity of R2* to local iron concentration as compared to the combined use of R2* and quantitative susceptibility mapping (QSM) [18].

Volumetric Data

CNS atrophy occurs in all stages of MS, since the preclinical phase of the disease and progresses throughout its course, at a much higher rate than one reported in normal aging [83-86]. In this study, the annual brain percentage volume loss at the group level was 0.67%, which is in line with previous publications [87]. We also showed a significant increase in lesion fraction. Volumetric data (ARoC’s) were highly variable across subjects: changes in BPF range from -2.52 to 1.17% and LF from -0.78 to 103.06%. This variability arises from a large number of factors which do not necessarily relate to MS: age, disease duration, disease phenotype, disease modifying treatment, circadian rhythm, hydration… [85, 86]. Moreover, annual changes in brain parenchymal fraction as well as lesion fraction only partially correlated to patients’ disease status, in accordance with a large amount of publications [60, 88]. This highlights the lack of specificity and sensitivity of volumetric measurements, at least at the individual level.

Study limitation

As mentioned earlier, the small size and heterogeneous aspect of the present dataset constitute major limitations of this study. Indeed, it is composed of only 17 patients, with a rather broad range of characteristics such as age, disease duration, disease phenotype, disease modifying treatment, etc., which are known to influence the disability state of the patient and thus their ability to put together repair mechanisms within cerebral tissues [1, 68-70, 89, 90]. In addition, the time interval between two scanning sessions varied rather widely across patients (between 14 and 61 months), although it was brought back to an annual rate where possible. All of these parameters were imposed by standard clinical follow up. Therefore, these results should not be over-interpreted but are nevertheless promising and call for a replication with a larger and more homogeneous or controlled set of MS patients. Larger longitudinal studies are currently being held and will probably confirm these preliminary results.

The second significative limitation is the absence of longitudinal MRI data acquired in a control group of healthy subjects. However, we considered that literature of longitudinal studies of healthy subjects that analysed tissue microstructure could constitute a solution for comparison with MS patients. For example, in Bonnier et al. (2017) [56], the control group did not show any significant differences regarding T1, T2* or MTR measurements over two years, and the median age of their group is quite similar to ours (34,3 vs 36 years). Also, in Elkady et al. (2018) [57], they found no longitudinal R2* effect in their control groups, even with an age range superior to ours. Moreover, the median age of our population (< 60 years), as well as the short period between two scanning sessions (median of 14 months), suggests that microstructural alterations would not be noticeable in a healthy participants group, as many quantitative ageing studies detected differences over much larger time periods [26, 29, 91].