Voxel‐Based quantitative MRI reveals spatial patterns of grey matter alteration in multiple sclerosis

Despite robust postmortem evidence and potential clinical importance of gray matter (GM) pathology in multiple sclerosis (MS), assessing GM damage by conventional magnetic resonance imaging (MRI) remains challenging. This prospective cross‐sectional study aimed at characterizing the topography of GM microstructural and volumetric alteration in MS using, in addition to brain atrophy measures, three quantitative MRI (qMRI) parameters—magnetization transfer (MT) saturation, longitudinal (R1), and effective transverse (R2*) relaxation rates, derived from data acquired during a single scanning session. Our study involved 35 MS patients (14 relapsing–remitting MS; 21 primary or secondary progressive MS) and 36 age‐matched healthy controls (HC). The qMRI maps were computed and segmented in different tissue classes. Voxel‐based quantification (VBQ) and voxel‐based morphometry (VBM) statistical analyses were carried out using multiple linear regression models. In MS patients compared with HC, three configurations of GM microstructural/volumetric alterations were identified. (a) Co‐localization of GM atrophy with significant reduction of MT, R1, and/or R2*, usually observed in primary cortices. (b) Microstructural modifications without significant GM loss: hippocampus and paralimbic cortices, showing reduced MT and/or R1 values without significant atrophy. (c) Atrophy without significant change in microstructure, identified in deep GM nuclei. In conclusion, this quantitative multiparametric voxel‐based approach reveals three different spatially‐segregated combinations of GM microstructural/volumetric alterations in MS that might be associated with different neuropathology.

Importantly most of these studies did not characterize the spatial distribution of GM microstructural alterations when considering several quantitative parameters. In this paper, we precisely provide a whole-brain voxel-based quantification (VBQ) of three qMRI parameters-MT saturation, R1 and R2*, derived from data acquired in a single MR session-and assess the spatial distribution of their changes in a cross-sectional study which contrasted MS patients to HC. Potential GM atrophy was also investigated by a concurrent voxel-based morphometry (VBM) analysis.
Our study aimed at characterizing the spatial distribution of microstructural and volumetric GM alterations induced by MS at the regional level (i.e., voxel-wise). Results are presented as different spatial combinations of atrophy and microstructural damages in cortical and deep GM involvement, in MS compared with healthy controls.

| Participants
Seventy-two participants were initially included in the study (Lommers et al., 2019) (Lommers et al., 2019). This study was approved by the local ethic committee (B707201213806) and written informed consent was obtained from each participant.

| MR image acquisition and spatial processing
MRI data were acquired on either of the following 3 T MRI-scanners: Magnetom Allegra and Magnetom Prisma, Siemens Medical Solutions, Erlangen, Germany. The whole-brain MRI acquisitions included a multiparameter mapping (MPM) protocol that has been gradually optimized and validated for multi-centric acquisitions (Leutritz et al., 2020;Tabelow et al., 2019;Weiskopf et al., 2013). It consists of three co-localized series of 3D multi-echo fast low angle shot (FLASH) acquisitions at 1 × 1 × 1 mm 3 resolution and two additional calibration sequences to correct for inhomogeneities in the RF transmit field (Lutti et al., 2012). The FLASH data sets were acquired with predominantly proton density (PD), T1, and MT weighting, referred to in the following as PDw, T1w and MTw echoes. Volumes were acquired in 176 sagittal slices using a 256 × 224 voxel matrix. Details of the MPM protocol used for this study are available as supplementary data. An additional FLAIR sequence was recorded with spatial resolution 1 × 1 × 1 mm 3 and TR/TE/TI = 5,000 ms/516 ms/1800 ms.
All data analyses and processing were performed in Matlab (The MathWorks Inc., Natick, MA) using SPM12 (http://www.fil.ion.ucl.ac. uk/spm) and its extensions. MT saturation, R1 and R2* quantitative maps were estimated using the hMRI toolbox (http://hmri.info/) as previously described (Tabelow et al., 2019). Briefly, echoes for T1w, PDw, and MTw were extrapolated to TE = 0 to increase the signal-tonoise ratio and get rid of the otherwise remaining R2* bias (Tabelow et al., 2019). The resulting MTw and T1w (TE = 0) images were used to calculate MT saturation and R1 quantitative maps. To maximize the accuracy of the R1 and MT saturation maps, inhomogeneity in the flip angle was corrected by mapping the B1 transmit field according to the procedure detailed in (Lutti et al., 2012). In addition, intrinsically imperfect spoiling characteristics were accounted for and corrected in R1 map, using the approach described previously (Preibisch et al., 2009). The MT saturation map differs from the commonly used MT ratio (MTR, percent reduction in steady state signal) by explicitly accounting for spatially varying T1 relaxation time and flip angles. MT saturation shows a higher brain contrast to noise ratio than the MTR, leading to improved and more robust segmentation in healthy subjects (Helms, Dathe, & Dechent, 2010). The R2* map was estimated from all three multi-echo series using the ESTATICS model (Tabelow et al., 2019). Example whole brain maps are shown in Figure 1. Note that these MR sequences at 3 T are not sensitive enough to detect focal cortical lesions, as previously described (Hulst & Geurts, 2011).
Quantification of cortical parameters is thus possibly confounded by voxels located within cortical plaques.
MR images multi-channel segmentation and normalization were performed with the standard "unified segmentation" (US) approach for the HC and its US-with-Lesion extension, accounting for WM lesions, for the MS patients. This part of the processing is largely detailed in a previous publication (Lommers et al., 2019). Briefly, for each MS patients, a preliminary lesion mask was derived from the FLAIR images and used to update the "tissue probability maps" with an extra lesion tissue class limited to the WM. This patient specific extended TPM was then used in the US tool, therefore accounting for the usual brain and head tissue plus the lesions (Phillips, Lommers, & Pernet, 2016). Individual lesion fraction (LF, ratio of WM lesion load to total intracranial volume) was computed afterwards from the segmented tissue classes. For VBM analyses, GM probability map (including cortical and deep GM) were spatially warped to standard space, modulated by the Jacobian determinants of the deformations, and smoothed with an isotropic Gaussian kernel (6 mm full width at half maximum-FWHM). For VBQ analyses, the 3 quantitative maps were normalized using the subject-specific deformation field but without modulation. A tissue weighted smoothing (3 mm FWHM isotropic) yielded smoothed tissue-specific multiparameter maps which optimally preserved quantitative parameter values within each tissue class (Draganski et al., 2011). Detailed analysis of the influence of spatial deformations onto quantitative parametric values proved this method to be largely insensitive to volumetric changes (i.e., atrophy) (Salvoni et al., 2019). Finally, a GM mask was generated: the smooth modulated warped individual GM, WM and CSF maps were averaged across all subjects and the GM mask included voxels for which mean GM probability was larger than that of WM or CSF and exceeded 20% (Callaghan et al., 2014).

| Statistical analyses
Whole-GM voxel-wise VBM and VBQ statistical analyses, explicitly using the GM mask, were carried out using a multiple linear regression model embedded in the general linear model framework of SPM12.
MRI data were analyzed in a factorial design, with the 2 different scanners as one factor and the group (MS vs. HC) as the second factor (Stonnington et al., 2008). Age, gender and total intracranial volume were entered as covariates of no interest. Differences between MS patients and HC as well as interactions between groups and scanners were tested by separate F-tests for each quantitative parameter (MT, R1, R2*) and volume.
Post hoc t tests explored significant effects. Cluster-level inferences were conducted at p < .05 after family-wise error rate (FWER) correction for multiple comparisons across the whole GM (p < .0001 uncorrected cluster-defining threshold). These 2-sample t-tests identified significant group effects, over and above the normal spatially heterogeneous distribution of quantitative parameters (Deistung et al., 2013) and accounting for potential unequal variance across groups.
In the patient population, three F tests looked for significant voxel-wise regression between each qMRI parameter and clinical scores (EDSS, motor and cognitive composite scores) as well as lesion fraction. Significance threshold was set at p < .05 FWER corrected at cluster level (p < .0001 uncorrected cluster-defining threshold).

| RESULTS
Compared with HC, we identified significant loco-  supposedly because of their intrinsically heavy metabolic load (Calabrese et al., 2015) and their frequent involvement in focal WM inflammation. Finally, because of their numerous folds, these regions are more exposed to cerebrospinal fluid (CSF) stasis, supporting the hypothesis that soluble factors produced in the CSF by lymphocytes influence subpial demyelination, particularly in patients with progressive MS (Magliozzi et al., 2018).

| Pattern 2: Hippocampus
The evidence of substantial demyelination of hippocampi beyond atrophic areas constitutes a key contribution of this study and usefully complements previous characterization of hippocampal damage in MS (Rocca et al., 2018). Indeed, demyelination is detected postmortem in 53 to 79% of MS hippocampi (Dutta et al., 2011;Dutta et al., 2013;Geurts et al., 2007)  inconsistently observed in demyelinated hippocampi while synaptic density is systematically decreased (Dutta et al., 2011;Geurts et al., 2007;Papadopoulos et al., 2009). By the same token, chronic inflammation potentially enhances neurogenesis within dentate gyrus (Rocca et al., 2015). Although the functional significance of these cellular changes is still under debate (Pluchino et al., 2008;Zhao, Deng, & Gage, 2008), they may balance neuronal loss, at the structural level (Rocca et al., 2015).

| Pattern 3: Deep gray matter nuclei
Our results show a significant atrophy of DGM in MS and agree with previous reports (Hulst & Geurts, 2011). Thalamic and putamen atrophy relates to significant neuronal and axonal loss. It occurs very early in the disease course and exceeds cortical atrophy (Eshaghi et al., 2018). Due its extensive reciprocal connections with cortical and subcortical structures, thalami are particularly vulnerable to anterograde and retrograde degeneration. This interpretation is supported by the significant inverse relationship between each qMRI parameter value (MT, R1, and R2*) within thalami and the lesion load, which suggests that lesions in connecting WM tracts also alter thalamic microstructure. These neurodegenerative processes likely dominate local inflammatory activity and oxidative injury which were also reported (Haider et al., 2016) but were not sensitively assessed in this study.

| Limitations
This cross-sectional study was run on a relatively small sample size.
To preserve statistical power, the two MS phenotypes were pooled together. We cannot rule out that results are partly driven by the larger proportion of PMS over RRMS patients, although exploratory tests did not show any significant difference between the two patient groups when corrected for multiple comparisons. If considering results significant at p < .001 uncorrected for multiple comparisons, PMS patients showed reduction in MT, R1, and R2* in a number of regions that were already detected in the contrast involving healthy control and the whole MS population. This might indicate that demyelination is even more severe as disease progresses. Because differences in microstructure between RRMS and PMS patients are of paramount importance, they will be assessed in future work, based on larger and independent population samples.
The factorial design used in this study offers the possibility to test separately the effect of disease and the effect of scanner, as well as their interaction. The absence of significant group by scan interaction allows us to disentangle the potential confound introduced by the two different scanners and still reliably discuss the effect of disease on GM microstructure. Furthermore, acquisition protocol has been optimized, pointing out the opportunity for multi-centric studies (Leutritz et al., 2020).
Finally, our results do not confirm previous reports linking thalamic and hippocampal damage to motor performance and cognitive dysfunction in MS patients (Eshaghi et al., 2018;Rocca et al., 2018).
Inferences were conservatively made after correction for multiple comparisons over the whole GM, increasing the risk of Type II error.
In this preliminary study, we indeed considered that conservative inferences had to be preferred to spurious results. Alternatively, it might be the case that microstructural alterations precede the occurrence of clinical symptom: longitudinal studies are needed to answer this question. Moreover, spinal cord lesions were not taken into account although they impact motor performance.

| CONCLUSION
This multiparametric voxel-based approach identifies three different spatially-segregated patterns of GM microstructural/volumetric alterations in MS patients, that might be associated with different neuropathology. The results highlight the usefulness of qMRI parameters and their complementarity with volumetric techniques in assessing GM status in MS.

ACKNOWLEDGMENT
The authors are particularly thankful for the patients and healthy participants who eagerly took part in this study.

CONFLICT INTERESTS
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

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.