Using quantitative magnetic resonance imaging to track cerebral alterations in multiple sclerosis brain: A longitudinal study

Abstract Introduction Quantitative MRI quantifies tissue microstructural properties and supports the characterization of cerebral tissue damages. With an MPM protocol, 4 parameter maps are constructed: MTsat, PD, R1 and R2*, reflecting tissue physical properties associated with iron and myelin contents. Thus, qMRI is a good candidate for in vivo monitoring of cerebral damage and repair mechanisms related to MS. Here, we used qMRI to investigate the longitudinal microstructural changes in MS brain. Methods Seventeen MS patients (age 25–65, 11 RRMS) were scanned on a 3T MRI, in two sessions separated with a median of 30 months, and the parameters evolution was evaluated within several tissue classes: NAWM, NACGM and NADGM, as well as focal WM lesions. An individual annual rate of change for each qMRI parameter was computed, and its correlation to clinical status was evaluated. For WM plaques, three areas were defined, and a GLMM tested the effect of area, time points, and their interaction on each median qMRI parameter value. Results Patients with a better clinical evolution, that is, clinically stable or improving state, showed positive annual rate of change in MTsat and R2* within NAWM and NACGM, suggesting repair mechanisms in terms of increased myelin content and/or axonal density as well as edema/inflammation resorption. When examining WM lesions, qMRI parameters within surrounding NAWM showed microstructural modifications, even before any focal lesion is visible on conventional FLAIR MRI. Conclusion The results illustrate the benefit of multiple qMRI data in monitoring subtle changes within normal appearing brain tissues and plaque dynamics in relation with tissue repair or disease progression.


INTRODUCTION
Multiple sclerosis (MS) is a chronic autoimmune disease of the central nervous system (CNS). Plaques are the pathological hallmark of MS. They are spread in acute, focal, disseminated, and recurrent way throughout the CNS and harbor variable degrees of inflammation, demyelination, gliosis, and axonal injury (Lassmann, 2013;Trapp et al., 1998). Plaques are not restricted to the white matter (WM) but are also present in the cortex and deep gray matter (GM) (Haider et al., 2016;Kutzelnigg et al., 2005;Popescu et al., 2012).
Over and above focal WM lesions, an early, diffuse, and chronic inflammation within the normal appearing white matter (NAWM) and gray matter (NAGM) is ultimately responsible for diffuse neuroaxonal loss and neurodegeneration, which is deemed responsible for a progressive accumulation of disability (Frischer et al., 2009;Haider et al., 2016;Kutzelnigg et al., 2005).
By contrast, effective repair mechanisms can occur within focal lesions but probably also in normal appearing brain tissue (NABT) (Brown et al., 2014). However, our understanding of these complex processes is still fragmentary. The difficulty of acquiring histopathological data on MS patients at various stages of the disease makes it challenging to describe the time course of injury and potential repair mechanisms in MS. Consequently, there is a need for new imaging techniques to improve the in vivo monitoring of brain damages formation, progression and repair in MS (Wang et al., 2019).
Conventional MRI (cMRI) readily depicts focal WM lesions on T2/FLAIR sequences and is able to distinguish between acute and allegedly chronic lesions. T2-hyperintensities in cMRI constitute the keystone of McDonald diagnostic criteria (Thompson et al., 2018) and also make an important contribution to the monitoring of WM lesion burden. Unfortunately, cMRI sequences do not sensitively detect cortical lesions and diffuse changes in NABT, due to a rather low sensitivity of cMR imaging for cortical lesions, mixed contrast weight, and an overall limited histopathological specificity within cerebral tissues.
Quantitative MRI (qMRI) potentially overcomes these limitations by quantifying physical microstructural properties of cerebral tissue in standardized units. qMRI is more sensitive but also more specific to microstructural properties of CNS tissues. Magnetization transfer ratio (MTR) was regularly linked to cerebral macromolecular content detected by a greater percentage loss of magnetization in voxels with a higher myelin content and axons density Schmierer et al., 2004;Tabelo et al., 2019). Postmortem studies comparing the relative contribution of these two factors indicate that myelin has a stronger and more direct influence on MTR than the axonal density, which is considered as a T1-dependent effect. Tissue water content (inflammation, edema, etc.), another T1-dependent effect, also accounts for MTR variability (Mottershead et al., 2003;Schmierer et al., 2004;Van Waesberghe et al., 1999). However, the MT saturation (MTsat) map offers a measure which, unlike MTR, is minimally affected by longitudinal relaxation and B1 mapping inhomogeneities (Lema et al., 2017), increasing its sensitivity to myelin content.
Moreover, the brain contrast to noise ratio is larger for the MTsat map than for MTR, thus improving brain tissue segmentation in healthy subjects Schmierer et al., 2004). R2* was usually linked to iron and myelin contents, as paramagnetic iron and diamagnetic myelin generate microscopic field gradients in the CNS, thus shortening T2* and increasing R2* (1/T2*). Orientation and density of myelin fibers are also a determining factor of R2* values (Bagnato et al., 2018;Hametner et al., 2018;Stüber et al., 2014). Iron is probably a key factor in MS monitoring as it was shown that aberrant iron metabolism occurs in the course of the disease (Stankiewicz & Neema et al., 2014).
Regarding focal WM plaques, qMRI emerges as an appealing biomarker to describe the dynamic processes of demyelination and remyelination. For instance, MTR was shown to sharply decrease within gadolinium enhancing lesions before recovering during the subsequent months (Dousset et al., 1998;Elskamp et al., 2010;Levesque et al., 2010), and within NAWM days to weeks before the formation of a new active lesion (Fazekas et al., 2002;Filippi et al., 1998).
Because each qMRI parameter is differently sensitive to histologically measured iron and myelin contents, this approach might become a fundamental tool for longitudinal in vivo monitoring of MS lesions and NABT evolution at the tissue microstructural level.
In this longitudinal study, we investigate the evolution of four simultaneously acquired qMRI parameters (MTsat, PD, R1, R2*) within NABT and WM lesions of 17 MS patients-relapsing remitting (RRMS) and progressive MS (PMS)-who were scanned two times with at least a 1-year interval, following the same multiparameter mapping (MPM) protocol at 3 Tesla (Draganski et al., 2011;Tabelo, 2019

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 multiparameter mapping protocol (MPM), from which one can simultaneously estimate (semi)quantitative maps of magnetization transfer saturation (MTsat), proton density (PD), longitudinal relaxation (R1), and effective transverse relaxation (R2*). This protocol arising from an international collaborative effort (Draganski et al., 2011;Tabelo, 2019) has already been used to study brain microstructure in various conditions including normal aging (Carey et al., 2018;Draganski et al., 2011;Thompson et al., 2018), brain tumor (Reuter et al., 2020), Parkinson's disease (Depierreux et al., 2021;Klein et al., 2018;Nürnberger et al., 2017), as well as MS. It consists of three colocalized 3D multiecho fast low angle shot (FLASH) acquisitions at 1 mm 3 resolution and two additional calibration sequences to correct for inhomogeneities in the RF transmit field (Lutti et al., 2010;Lutti et al., 2012). 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 1 mm 3 and TR/TE/TI = 5000 ms/516 ms/1800 ms.
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 presubtracted phase image, were acquired for image distortions corrections. A summary of the acquisition parameters appears in Supplementary data.
Note that these MR sequences at 3 Tesla are not sensitive to cortical lesion as described in Filippi et al. (2013) and Hulst and Geurts (2011) although a few lesions at the corticosubcortical border were detected. Quantification of cortical parameters is thus confounded by voxels potentially located within cortical lesions.

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 Segmenta- Baseline EDSS, median (range) 2.5 (1-6.5) Baseline number of relapses, median (range) RRMS: 2 (1-5) PMS: N/A Disease-modifying treatment RRMS: first line, n: 5 second line, n: 6 PMS: ocrelizumab, n: 2 none, n: 4 F I G U R E 1 Chart flow of data creation and processing (see text). MPM maps were created with the hMRI toolbox, and 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 coregistered to the MPM T0 space. Segmentation using USwL allowed to isolate the different tissue classes.
spoiling using the strategy of Preibisch and Deichmann (2009). The receive bias field map (B1-) was used to correct PD maps for instrumental biases. The R2* map was estimated from all three multiecho series (MTw, PDw, and R1w) using the ESTATICS model .
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 (Ashburner & Ridgway, 2013). 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 (Schmidt et al., 2012), followed by manual corrections by an MS expert (EL) to remove aberrant/artifactual lesion detections (Lommers et al., 2019).
The images were then segmented using the USwL toolbox, which consists of an extended version of the traditional Unified Segmentation (US) algorithm (Ashburner & Friston, 2005) and includes an additional tissue class representing the WM lesion(s). The US-with-lesion method internally generates a subject-specific extended set of tissue probability maps (TPM) (Lorio et al., 2016): an extra tissue class, based on the smoothed preliminary lesion mask warped into template space (using cost function masking during normalization; Andersen et al., 2010), is added to account for the lesion, and the original white matter prior map is updated accordingly (Moon et al., 2002). The gray matter TPM was not updated due to a very low number of lesions present in the cortical ribbon. Multichannel 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 < 10 mm 3 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.

Brain volume change
Volumetric changes were investigated using the USwL a posteriori tissue probability maps. The following measures of brain volume were

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: where Y is the ARoC for a qMRI parameter and tissue class, β's are the regression parameters corresponding to the associated regressor (with 0 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 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 Table 2.
Additional clinical information concerning annual relapse rate and treatment administration appears in Supplementary material. Permutation tests were employed for inferences (Anderson, 2001).
R 2 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 X status ), tested against the unchanged response Y, and generating each time a permuted R 2 value (noted R , R obs being the true R 2 value computed without permutation of the data). The condition X status ≠ X status is verified at each permutation. After n permutations (with n = 5000 in this study), a p value was computed based on the following formula:  (Benjamini & Hochberg, 1995), 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, that is, 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.

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 segmentation, 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 segmentation 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, pertaining to hyperintensity on the conventional FLAIR MR image acquired at T0. These are referred to as "focal FLAIR lesion." • The peripheral region detected on qMRI maps at T0, at the borders of (but not including) the focal FLAIR lesion. Those are referred to as "initial peripheral lesion." • The peripheral region, detected on qMRI maps at follow-up, bordering (but not including) the initial peripheral lesion, further referred to as "later peripheral lesion." This was computed by masking out the T1 lesion mask with the T0 lesion mask. This region allows us to determine whether its microstructure at T0 forebodes a fullblown plaque, detectable during follow-up. Those sometimes appear hyperintense on FLAIR images.
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, that is, present in the three masks, were considered for these comparisons.
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 F I G U R E 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). 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.

Analysis of normal appearing tissues
As expected, changes in MTsat and R2* within normal appearing brain tissues (NABT) between T0 and T1 varied across subjects (Figure 3) ( Table 3): MTsat in NAWM and NACGM and R2* in NAWM significantly increased in patients who fare well (X status = 1).
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.
Regarding BPF and LF, their correlation to the disease activity status was not significant (

Analysis of lesion microstructure
The number of enlarging WM lesions between T0 and T1 varied from At times 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 GLMMs appear in Supplementary data.

DISCUSSION
This longitudinal study followed up volumetric data and qMRI brain metrics (MTsat, PD, R1, R2*) in 17 patients with multiple sclerosis F I G U R E 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.
for a median time interval of 30 months. The main results are threefold.
First, the microstructure of normal appearing brain tissues changes over time and these modifications concur with, and potentially drive, clinical evolution. This critical finding suggests that repair mechanism and edema resorption can be monitored in vivo. Second, the microstructure within WM plaques is remarkably heterogeneous.
Importantly, at their periphery, microstructural alterations foreshadow their expansion, as detected by conventional MRI. Third, as expected, we observed a small but significant brain atrophy and lesion load increase with time.

Quantitative MRI parameter time course within NABT
In this study, we used a multiparameter mapping protocol that was gradually optimized and validated for multicentric studies (Leutritz et al., 2020). 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 (Andica et al., 2019;Bonnier et al., 2014;Engström et al., 2014;Gracien et al., 2016;Lommers et al., 2019;Lommers et al., 2021;Neema et al., 2007;Reitz et al., 2017;Stevenson et al., 2000). By contrast, longitudinal analyses of multiparameter qMRI data are scarce. A progressive shortening of T2/T2* (Bonnier et al., 2017) or increase in R2* (Elkady et al., 2018;Elkady et al., 2019;Khalil et al., 2015) 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 . MTR progressively decreases in NAWM of MS patients over 1 (Laule et al., 2003) or 2 years (Hayton et al., 2012). These abnormalities tend to be more pronounced in progressive phenotypes (Rocca et al., 1999) 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 both MTsat and R2* correlate with myelin content Carey et al., 2018;Hametner et al., 2018;Mangeat et al., 2015;Schmierer et al., 2004;Weiskopf et al., 2015), 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 (Lommers et al., 2019). 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 (Gracien et al., 2016;Lommers et al., 2019;Neema et al., 2007) and longitudinal  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 (Dutta & Trapp, 2007;Frischer et al., 2015;Lassmann et al., 2001). 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 (Frischer et al., 2015;Lucchinetti et al., 2000). This accounts for the large inter-and intraindividual 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 (Dousset et al., 1998;Elskamp et al., 2010;Levesque et al., 2010). Likewise, reduction of MTR within NAWM, days to weeks before the formation of a new active lesion, was also demonstrated (Fazekas et al., 2002;Filippi et al., 1998), and long-term MTR changes in WM plaques were observed in relation with disease progression (Rocca et al., 1999;Zheng et al., 2018). The present study broadens the quantitative characterization of plaque dynamics, in keeping with previous longitudinal studies (Bonnier et al., 2017;Chawla et al., 2018). 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 (Lieury et al., 2014). 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 (Frischer et al., 2015;Kuhlmann et al., 2017;Lassmann et al., 2007;Lucchinetti et al., 2000), 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 postgadolinium sequences. If confirmed on larger population samples, this finding might significantly modify treatment management in MS patients.
Another plausible hypothesis explaining the progressive decrease of R2* in initial and later peripheral regions is that iron-containing macrophages could be removing iron from the lesions through perivascular drainage into the extracellular compartment. Previous neuropathological studies have reported an iron loss at the edges of a subset of MS lesions, depending on their type (active, inactive, smoldering, etc.) as well as the patient's age and disease duration Popescu et al., 2017). Due to the limited sensitivity of R2* to local iron concentration as compared, for example, to the combined use of R2* and quantitative susceptibility mapping (QSM) , validating this theory would require additional measures, which can better describe iron dynamics in MS lesions and NAWM.

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 (Bermel & Bakshi, 2006;De Stefano et al., 2010;Eshaghi et al., 2018;Zivadinov et al., 2008). In this study, the annual brain percentage volume loss at the group level was 0.67%, which is in line with previous publications (De Stefano et al., 2015). We also showed a significant increase in lesion fraction.  (Bermel & Bakshi, 2006;Zivadinov et al., 2008).
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 (Enzinger et al., 2015;. This highlights the lack of specificity and sensitivity of volumetric measurements, at least at the individual level. It can appear odd that brain atrophy progresses in parallel to repair mechanisms, as suggested by qMRI parameters. However, BPF reduction is minimal and is not significant (see Table 3) between T0 and T1.
One should keep in mind that cortical atrophy is an irreversible phenomenon. Given the inter-and intraindividual heterogeneity of MS progression, it is possible that patients who have undergone neuronaxonal loss at some point in the disease might be able to remyelinize their remaining axons, hopefully through therapeutic intervention or lifestyle changes. Besides, axonal remyelination is not always effective.
Here we showed that variations in MTsat and R2* correlated to the disease activity status, but our clinical evaluation based on EDSS is undoubtedly imprecise. Once again, the size and heterogeneity of our cohort limits the interpretation of such results.

Study limitation
As mentioned here above, 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 (Bodini et al., 2016;Frischer et al., 2015;Lassmann, 2013;Lassmann et al., 2001;Lucchinetti et al., 2000;Patrikios et al., 2006). 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.
A second 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 analyzed tissue microstructure could constitute a solution for comparison with MS patients. For example, in Bonnier et al. (2017), the control group did not show any significant differences regarding T1, T2*, or MTR measurements over 2 years, and the median age of their group is quite similar to ours (34.3 vs. 36 years). Also, in Elkady et al. (2018), 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 aging studies detected differences over much larger time periods Draganski et al., 2011;.

CONCLUSION
These preliminary results highlight the relevance of multiple qMRI data in the monitoring of MS disease, highlighting subtle changes within NABT and plaque dynamics in relation with repair or disease progression. Of course, large-scale longitudinal study would be needed to reproduce these findings and better exploit the full potential of qMRI parameters.

ACKNOWLEDGMENTS
The authors are particularly thankful for the patients who took part in this study.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.