Participants
For this cross-sectional study, we included a cohort of 271 patients with relapsing-remitting or secondary progressive MS according to the 2010 McDonald criteria24 prospectively recruited at the MS Unit in the Hospital Clinic of Barcelona, as well as 54 healthy controls (HC) without any prior or present history of neurological or psychiatric conditions. To explore GM changes at different time-points of the disease, patients were classified into three groups according to their disease duration: less than 5 years from disease onset (MS1); from 5 to 15 years (MS2); and more than 15 years of disease duration (MS3)25.
Physical disability was evaluated in all patients using the Expanded Disability Status Scale (EDSS)24, and in 214 of them, attention and information processing speed was assessed using the Symbol Digit Modalities Test and Paced Auditory Serial Addition Test from the Brief Repeatable Battery of Neuropsychological tests26. We calculated a mean z-score of these two tests (zAttention-Processing speed) using the age-and education-adjusted values of a healthy Spanish cohort27.
The Ethics Committee of the Hospital Clinic of Barcelona approved the study and all the participants signed an informed consent before enrolment onto the study. All procedures were performed according to the principles of the Helsinki Declaration.
MRI acquisition and processing
Structural and diffusion MRI acquisition
MRI was performed on a 3T Magnetom Trio (SIEMENS, Erlangen, Germany) scanner with a 32 channel phased-array head coil. The protocol included structural 3D-Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE), 3D-T2 fluid-attenuated inversion recovery (FLAIR) and diffusion-weighted imaging (DWI) sequences. DWI was acquired with two different sequence parameters. A detailed description of all sequences used in this study is available in the Supplementary material10,28.
Grey matter volume processing
Before computing GM volume, WM lesions were segmented in 3D-MPRAGE space with ITK-SNAP Software (http://www.itksnap.org/pmwiki/pmwiki.php). We used a linear registration of the 3D-FLAIR image to align with 3D-MPRAGE to improve both the identification and delineation of lesions. Subsequently, we applied a WM lesion-filling approach29 and obtained 76 GM regions using the Mindboggle software (https://mindboggle.info), applying the Desikan-Killiany-Tourville cortical labelling atlas (31 cortical labels per hemisphere)30, and the automated subcortical segmentation offered by the FSL-FIRST package for seven subcortical regions in both hemispheres31. The resulting cortical surface parcellation was visually checked and fixed to guarantee the quality of the cortical segmentation labels for further statistical analysis. We excluded 24 patients due to segmentation errors in the cortical surface reconstruction. Finally, we normalised the 76 GM volumes using a volumetric scaling factor provided by SIENAX to reduce the effect of the variability in head-size for quantification32. The flowchart for GM volumetric processing is summarised in Figure 1 (left panel).
Grey matter diffusion processing
The diffusion MRI pre-processing pipeline was described in detail elsewhere23,28 and it involved the following steps: DWI denoising, Gibbs ringing correction, motion-induced distortion correction, a phase unwrapping procedure to correct for geometric distortion and bias field correction33. Subsequently, the MD DTI index was obtained34.
The GM cortical reconstruction process was performed using in-house surface-based approach developed to measure the microstructural changes in neurodegenerative disorders12,13,35. We used Freesurfer tools to calculate the midpoint of each individual’s cortical surface in order to avoid contamination by adjacent WM and cerebrospinal fluid12,13. We then registered the T1-weighted image to the pre-processed diffusion data using a boundary-driven algorithm to realign both the cortical surface mask generated and the 76 GM labels on the DTI maps of each subject36. We selected only those voxels labelled as the midpoint cortical mask to obtain the average of the diffusion metrics for each label. We then eroded each of the 14 deep GM structures using a 3 mm box kernel to minimise the partial volume effects and computed the MD average within these masks. The GM diffusion processing pipeline is summarised in Figure 1 (right panel).
Diffusion and volumetric data were harmonised using the ComBat function in R software37,38 to reduce the variability between sequences.
Statistical analysis
We described the demographic, clinical and neuroimaging data with the mean and standard deviation (SD, for continuous data), or as absolute numbers and proportions (for categorical data), assessing the former’s normal distribution with histograms and Shapiro-Wilk tests. We analysed the demographic and clinical differences among groups with a Chi-squared or Kruskal-Wallis H test, using Dunn’s test when necessary, and we compared the GM diffusion and volumetric data with one-way ANOVA, using Tukey HSD test for the two-group differences (HC vs. MS1, MS2 and MS3). We described the effect size of these differences using Hedges’ g. To ensure that age differences between the HC and MS3 group did not confound the results, we additionally calculated the changes in the MRI metrics between these two groups with a Student t-test or Mann Whitney U Test selecting only the older HC subjects. Finally, in the whole sample of MS patients, we calculated Spearman’s correlations between the GM diffusion or volume and measures of disease severity (number of relapses, normalised lesion volume and EDSS), and the mean z-Attention-Processing speed score.
In all the analyses, we included age as a covariate to control for its potential influence on the results. We considered diffusion and volumetric values less than or greater than the 5th and 95th percentile as outliers, and removed them from the analysis. We used the false discovery rate (FDR) to correct for multiple comparisons, and we set the significance level at a corrected p<0.05. The R Statistical Software (version 3.6.1, www.R-project.org) was used for all statistical analyses.
Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.