Gray matter atrophy patterns in multiple sclerosis: A 10-year source-based morphometry study.

OBJECTIVES: To investigate spatial patterns of gray matter (GM) atrophy and their association with disability progression in patients with early relapsing-remitting multiple sclerosis (MS) in a longitudinal setting.
METHODS: Brain MRI and clinical neurological assessments were obtained in 152 MS patients at baseline and after 10 years of follow-up. Patients were classified into those with confirmed disability progression (CDP) (n = 85) and those without CDP (n = 67) at the end of the study. An optimized, longitudinal source-based morphometry (SBM) pipeline, which utilizes independent component analysis, was used to identify eight spatial patterns of common GM volume co-variation in a data-driven manner. GM volume at baseline and rates of change were compared between patients with CDP and those without CDP.
RESULTS: The identified patterns generally included structurally or functionally related GM regions. No significant differences were detected at baseline GM volume between the sub-groups. Over the follow-up, patients with CDP experienced a significantly greater rate of GM atrophy within two of the eight patterns, after correction for multiple comparisons (corrected p-values of 0.001 and 0.007). The patterns of GM atrophy associated with the development of CDP included areas involved in motor functioning and cognitive domains such as learning and memory.
CONCLUSION: SBM analysis offers a novel way to study the temporal evolution of regional GM atrophy. Over 10 years of follow-up, disability progression in MS is related to GM atrophy in areas associated with motor and cognitive functioning.

Introduction

Multiple sclerosis (MS) is a chronic, inflammatory disease of the central nervous system which was previously considered to primarily affect the white matter (WM). Over the last couple of decades though, an overwhelming amount of research has demonstrated that the gray matter (GM) is also heavily involved. Longitudinal studies have revealed robust associations between GM atrophy and several measures of disease progression such as the advancement of physical disability (Fisher et al., 2008, Jacobsen et al., 2014, Popescu et al., 2013) and cognitive impairment (Bergsland et al., 2016, Pitteri et al., 2017).

There is considerable evidence that GM atrophy begins already quite early in the disease, is clinically relevant, and can be quantified using advanced post-processing pipelines (Bergsland et al., 2012, Filippi et al., 2013, Fisher et al., 2008, Steenwijk et al., 2016, Zivadinov et al., 2013). It is not yet fully clear though whether differences in GM volume are already present in patients who will go on to progress in terms of clinical disability with respect to those who will not. In addition, many questions remain regarding spatial relationships of GM volume loss over time in MS patients. Although a study demonstrated that patterns of cortical atrophy are non-random and are associated with clinical outcomes (Steenwijk et al., 2016), it was limited in its cross-sectional nature. A better understanding of how GM atrophy evolves over time can only be gained with longitudinal studies. Moreover, advanced image processing techniques must be applied which are sensitive enough to capture the heterogeneous nature of GM tissue loss over time.

Against this background, the aims of the present study were two-fold. The first aim was to assess whether differences in GM volume were detectable already at baseline between patients who would subsequently develop disability progression over the course of 10 years of follow-up with respect to those who would not. The second aim was to identify spatial patterns of GM atrophy that are associated with disability progression. Therefore, we utilized an optimized source-based morphometry (SBM) technique, which allows for identifying patterns of covariation across the brain in terms of GM volume (Xu et al., 2009), in a cohort of early relapsing-remitting MS patients.

Materials and methods

Subjects

The current analysis is based on data from 152 patients with early relapsing-remitting MS that were part of the Avonex®-Steroids-Azathioprine (ASA) study (Havrdova et al., 2009), which enrolled a total of 181 subjects. Patients were scanned yearly over the course of ten years. Full details of the study design have been previously reported (Zivadinov et al., 2016).

The original inclusion criteria for the ASA study were clinically definite MS, according to the criteria of Poser et al. (Poser et al., 1983), at least 2 T2 lesions on baseline MRI, the presence of at least two oligoclonal bands in the cerebrospinal fluid, 18–55 years of age, EDSS score ≤ 3.5 on the day of screening, and active disease defined by two relapses in the past year or 3 relapses in the past 2 years. Inclusion criteria for the current sub-study were the availability of MR data at baseline and 10-year follow-up as well as information regarding disability progression.

Patients were divided into 2 groups based on confirmed disability progression (CDP) status at the end of the 10-year period. The first group contained patients who had CDP, defined as any 1.0-point sustained increase in EDSS score in patients who had a baseline EDSS score of ≥ 1.0 or any ≥ 1.5-point sustained increase in the EDSS score in patients who had a baseline EDSS score of 0. Increases in EDSS were confirmed after 48-weeks. The second group contained patients without disability progression. In addition, 31 healthy controls (HC), group-matched for age and sex, were also imaged at baseline.

The study protocol was approved by the Medical Ethics Committees of the General University in Prague and First Faculty of Medicine, Charles University in Prague; University of Buffalo; and ethics committees in the participating centers. All patients provided their written informed consent.

MRI acquisition and analysis

All MRI assessments were performed using the same Philips Gyroscan 1.5-Tesla scanner (Philips Medical Systems, Best, The Netherlands) that did not undergo major hardware upgrades over a 10-year period (average follow-up after 10.0 ± 0.25 years). Axial brain images were obtained using a 2D T2-weighted fluid-attenuated inversion recovery (T2-FLAIR) with 1.5-mm thickness and axial 3D T1-weighted (3D-T1) images were obtained with 1-mm slice thickness, both without inter-slice gap. Acquisition parameters for the sequences were as follows for FLAIR: echo time (TE) = 140 ms, repetition time (TR) = 11,000 ms, inversion time (TI) = 2600 ms, flip angle (FA) = 90°, field of view (FOV) = 256 mm); 3D-T1: TE = 5 ms, TR = 25 ms, FA = 30°, FOV = 256 mm. Both T2-FLAIR and 3D-T1 acquisitions remained unaltered over the course of the study. Although the study protocol included yearly MRI acquisitions, only the baseline and 10-year data was included for the current study.

T2-lesion volume (LV) was calculated using a reliable semi-automated edge detection contouring/thresholding technique with JIM (version 6) software. To reduce the impact of T1 hypointensities on subsequent analyses (Gelineau-Morel et al., 2012), 3D-T1 images were first lesion filled using co-registered T2-lesion masks. Normalized measures of whole brain and total gray matter volumes were obtained using SIENAX (Smith et al., 2002) while longitudinal changes were derived using SIENA (Smith et al., 2002) and SIENAX multi-time-point (Dwyer et al., 2014), respectively.

An optimized, longitudinal SBM pipeline was implemented via a combination of voxel-based morphometry (VBM) and independent component analysis (ICA) techniques. Advanced Normalization Tools (ANTs) (Avants et al., 2010) was used for all registrations, given that it was ranked number one in terms of accuracy in a comparison of 14 nonlinear, volumetric-based deformation algorithms (Klein et al., 2009). For each patient, non-linear registration of the baseline and 10-year 3D-T1 weighted images (WI) was used to generate an unbiased, subject-specific template. Native space 3D-T1s were then non-linearly registered to the corresponding template and averaged (3D-T1 average). The resulting image was deskulled using the brain extraction tool (BET) (Smith, 2002) and the resulting mask was transformed to the native 3D-T1 WI. Partial volume estimate maps of the cerebrospinal fluid, GM and WM were then obtained (Zhang et al., 2001). Using the 3D-T1 average from all 152 patients, a group template in standard 1-mm3 MNI space was constructed with using the ANTs “antsMultivariateTemplateConstruction2.sh” script with cross-correlation as the similarity metric and symmetric image normalization (SyN) as the transformation model. All registrations were visually checked to ensure that poor spatial normalization would not drive any of the results. A single warp was then generated for each image such that native GM segmentations were aligned into the space of the group template with only one transformation, resampled using trilinear interpolation and downsampled to 2-mm3. A montage of the resulting images is shown in Supplementary Video 1. The spatially normalized GM maps were then modulated by the Jacobian determinant of the corresponding warp and smoothed with an isotropic Gaussian kernel with σ = 3 mm. The resulting 304 images (152 subjects * 2 time points) were next fed into FSL's Multivariate Exploratory Linear Optimized Decomposition into Independent Components (MELODIC) (Beckmann and Smith, 2004) for ICA, which uses the FastICA algorithm, with the number of components (i.e. patterns) to extract set a priori to 8. Each map was thus based on the GM volume co-variation in all of the input images, which has the effect of identifying spatial patterns in terms of areas having similar GM concentrations. All of the resulting maps were visually inspected while component stability was assessed with ICASSO (Himberg et al., 2004) (1000 iterations). Only those components with a stability value of at least 0.95 were subsequently included in the rest of the analysis. ICA spatial maps were converted to Z statistic images via a normalized mixture–model fit and then thresholded at | Z | = 3 for visualization and reporting purposes. As a measure of GM volume for each pattern, we utilized the component loadings obtained from the ICA-derived mixing matrix, as previously described (Xu et al., 2009). SBM was also used to investigate GM volumetric differences between the MS cohort as a whole and HCs at baseline. A separate template was generated using only baseline images only from the MS and HC groups. All baseline images were normalized to this template and an additional SBM analysis was performed using the same procedures described above.

For comparison purposes, a VBM analysis was conducted on the same data that was processed with MELODIC for SBM. The “randomise” tool was used with threshold-free cluster enhancement (Smith and Nichols, 2009) with 5000 permutations while age, sex and head size were included as variables of no interest. Baseline differences in GM concentrations were assessed between 1) all MS patients versus HCs and 2) patients without CDP compared to those with CDP. The two patient groups were also compared in terms of GM atrophy over 10 years of follow-up. Results were considered significant at p < 0.05, family-wise error corrected.

Mapping of structures for both SBM and VBM analyses was done with the aid of the MNI Structural and Harvard-Oxford Subcortical/Cortical atlases provided with FSL.

Statistical analysis

All statistical analyses were carried out using SPSS 22.0 (IBM Corp., Armonk, NY) and R (3.2.1, R Foundation for Statistical Computing, Vienna, Austria). Differences at baseline between the groups in terms of demographic, clinical, and MR imaging characteristics were assessed using Student's t-test, Fisher's exact test, and Mann-Whitney U test, as appropriate. Analysis of Covariance (ANCOVA) was used to compare whole brain atrophy, total GM atrophy and T2-LV changes over the follow-up in MS patients while adjusting for age and sex. At baseline, ANCOVA models were used to compare SBM-derived loading factor differences between HCs and MS patients while adjusting for age, sex, and head size. For the MS patients only, repeated measures ANCOVA models were used to compare changes in ICA-derived loading factors between baseline and 10-year follow-up while adjusting for age, sex, and head size. The group by time interaction term was included in the model and used to assess when the change in loading factor was significantly different for the ICA-derived patterns. Profile plots were also examined in order to evaluate the direction of the effect of group on GM atrophy. Multiple comparisons were controlled for using the Benjamini-Hochberg false discovery rate method (Benjamini et al., 2001) and a corrected p-value (i.e. q-value) of < 0.05 was considered significant.

Results

Clinical and demographic measures at baseline and follow-up

Table 1 shows baseline clinical and demographic between MS patients without CDP (n = 67) and with CDP (n = 85). The CDP group was somewhat older on average (about 2.7 years, p = 0.035) but were otherwise similar in terms of clinical and demographical characteristics. As expected, the CDP group presented with a significantly higher EDSS score after 10 years of follow-up (p < 0.0001). Compared to the MS cohort as a whole, the HC group was about 1.5 years older, although this was not significant (p = 0.316) nor was the proportion of females (p = 0.644).

Table 1

Demographic and clinical characteristics of multiple sclerosis patients at baseline and change over follow-up.

	Without CDP group(n = 67)	CDP group(n = 85)	p-Value
Female, n (%)	51 (76.1)	69 (81.1)	0.548a
Age, mean (SD)	29.0 (7.3)	31.7 (8.1)	0.035b
Disease duration in years at baseline, mean (SD)	4.8 (3.4)	5.8 (5.8)	0.154b
Baseline EDSS, median (range)	2 (0–3.5)	2 (0–3.5)	0.287c
Follow-up EDSS, median (range)	2 (0–3.5)	4 (1.5–6.5)	< 0.0001c

Legend: CDP – confirmed disability progression; EDSS – Expanded Disability Status Scale.

Fisher's exact test.

Student's t-test.

Mann-Whitney U test.

MRI characteristics at baseline and follow-up

Table 2 shows baseline and follow-up changes of T2-LV, whole brain and total GM volumes in the MS cohort. No significant differences were detected between the groups at baseline, but the CDP group had significantly greater whole brain and total GM volume loss over 10 years of follow-up. At baseline, MS patients presented with significantly smaller whole brain and GM volumes compared to the HC group (p < 0.0001).

Table 2

MRI volumetry of multiple sclerosis patients at baseline and changes over follow-up.

	Without CDP group (n = 67)	With CDP group (n = 85)	HC group (n = 31)	Without CDP vs. CDP p-value	MS patients vs. HC p-value
Baseline T2 lesion volume, mean (SD)	6.1 (9.9)	9.1 (10.5)	–	0.080	–
Baseline whole brain volume, mean (SD)	1515.8 (76.3)	1502.4 (81.1)	1592.0 (87.5)	0.305	< 0.0001
Baseline gray matter volume, mean (SD)	813.2 (48.7)	805.2 (54.1)	847.1 (50.9)	0.349	< 0.0001
T2 lesion volume absolute change, mean (SD)	2.9 (5.5)	4.2 (7.5)	–	0.273	–
Percent brain volume change, mean (SD)	− 5.3 (3.0)	− 7.5 (3.8)	–	< 0.0001	–
Percent gray matter volume change, mean (SD)	− 5.8 (2.8)	− 7.1 (3.7)	–	0.007	–

Legend: CDP – confirmed disability progression.

Baseline p-values were calculated using Student's t test while changes over follow-up were derived using analysis of covariance, correcting for age and sex. Volumes are presented in mL.

Gray matter patterns at baseline and changes over follow-up

All patterns were deemed stable (minimum stability of 0.99) and thus were included in the statistical analyses. In general, the SBM-identified components corresponded to structurally or functionally connected GM regions. For the SBM analysis utilizing only MS patient data from both time points, details of the individual patterns are reported in Table 3 and shown in Fig. 1. No significant differences were detected at baseline between MS patients with CDP compared to those without after 10 years of follow-up. For the SBM analysis of baseline data from MS patients and HCs, details are reported in Table 4 and shown in Fig. 2. With respect to HCs at baseline, the mean loading factors from patterns 1 and 7 were significantly different after correction for multiple comparisons (q < 0.0001), reflecting widespread atrophy in MS patients throughout the deep GM and cortex. In the repeated measures analysis comparing loading factors between the two time points, the group by time interaction term was significant in patterns 1 and 7 after correcting for multiple comparisons. In general, the effects were seen bilaterally for all structures except where explicitly noted. Pattern 1 (q = 0.001) included the anterior and posterior cingulate gyri, precentral gyrus, supplementary motor cortex, insula, thalamus, caudate, nucleus accumbens, pallidum, and cerebellum. Pattern 7 (q = 0.007) included parahippocampal gyrus, lingual gyrus, hippocampus, nucleus accumbens, and putamen. Including baseline T2-LV or change over the follow-up as an additional covariate did not change any of the outcomes nor did controlling for treatment switch or relapse rate (results not shown).

Table 3

Independent component analysis-derived gray matter patterns identified in patients with multiple sclerosis at baseline and follow-up.

Pattern	Anatomic regions	Cluster extent	Peak Z statistic	X	Y	Z
1	L./R. anterior cingulate gyrusL./R. posterior cingulate gyrusL./R. precentral gyrusL./R. supplementary motor area	2990	8.08	0	− 16	46
	L./R. thalamusL./R. caudateL./R. nucleus accumbensR. pallidum	396	5.48	2	2	10
	Cerebellum	233	4.16	38	− 62	− 22
	L/R. insula	149	4.52	44	− 6	− 2
	R. thalamus	149	4.45	16	− 34	4
2	Cerebellum	6346	15.20	26	− 72	− 52
	L. thalamus	139	4.64	− 18	− 34	0
	R. thalamus	108	4.49	18	− 34	4
	L. superior parietal lobule	18	3.34	− 32	− 52	46
	L. precuneus	16	3.15	6	− 50	56
3	L. lateral occipital cortex	4930	8.04	− 36	− 86	20
	R. lateral occipital cortex	4319	6.39	32	− 84	18
	L./R. precuneusL./R. posterior cingulate gyrus	775	4.7	2	− 62	46
4	Cerebellum	9511	12.40	24	− 78	− 28
	R. nucleus accumbensR. pallidumR. putamenR. caudate	813	4.91	10	6	− 6
	L. nucleus accumbensL. pallidumL. putamen	776	5.24	− 12	8	− 8
	R. superior frontal gyrusR. supplementary motor area	143	3.96	20	− 6	72
	R. angular gyrusR. supramarginal gyrus	132	4.5	56	− 52	50
	L. supramarginal gyrus	39	3.74	− 60	− 44	44
5	CerebellumL./R. hippocampusL./R. thalamusL./R. lingual gyrus	3880	5.47	− 18	− 38	− 2
	L./R. precuneusL./R. postcentral gyrusL./R. superior parietal lobule	2893	5.61	4	− 48	72
6	L./R. caudateL./R. nucleus accumbensL./R. putamenL./R. pallidumL./R. thalamus	2876	5.16	− 4	10	0
	R. insula	1018	5.14	42	14	− 12
	L. insula	872	5.13	− 42	2	− 6
	L./R. frontal pole	659	4.87	− 2	64	4
	L./R. frontal poleL./R. anterior cingulate gyrus	239	4.85	30	60	6
	L./R. frontal pole	124	3.84	− 22	46	44
7	CerebellumL./R. hippocampusL./R. parahippocampal gyrus L./R. lingual gyrus	2757	5.44	− 8	− 56	− 8
	R. nucleus accumbensR. putamen	101	3.95	10	8	− 12
	L. nucleus accumbensL. putamen	100	3.96	− 10	8	− 14
8	L./R. intracalcarine cortexL./R. supracalcarine cortexL./R. lingual gyrusL./.R. precuneusL./R. cuneus	2593	5.51	18	− 64	8

Legend: L. – Left; R. – Right.

Cluster extent refers to the number of contiguous 2mm3 voxels. X, Y and Z refer to standard Montreal Neurological Institute coordinates.

Fig. 1

Source-based morphometry results showing the eight patterns of gray matter volume that were identified in all 152 multiple sclerosis patients at baseline and ten years of follow-up. Z-scores, thresholded at | 3 |, are shown in red-yellow for each component, with yellow being associated with larger values. Although no differences were identified at baseline, patients with confirmed disability progression (CDP) after ten years of follow-up had a significantly greater rate of atrophy in Patterns 1 and 7 compared to those without CDP, after correction for multiple comparisons. The X, Y and Z values refer to the slice shown in Montreal Neurological Institute coordinates.

Table 4

Independent component analysis-derived gray matter patterns identified in patients with multiple sclerosis and healthy controls at baseline.

Pattern	Anatomic regions	Cluster extent	Peak Z statistic	X	Y	Z
1	L./R. precuneus	419	4.59	0	− 52	62
	L. central opercular cortexL. insula	385	4.39	− 44	− 10	6
	R. insula	373	4.09	44	16	− 8
	R. central opercular cortex	297	5.40	44	− 10	8
	L./R. superior frontal gyrus	255	4.66	0	40	46
	L./R. caudateL./R. thalamus	209	5.67	− 10	4	6
	R. postcentral gyrusR. precentral gyrus	188	4.04	60	− 8	42
	R. precentral gyrus	132	4.15	46	8	34
	L. superior parietal lobule	96	4.01	− 34	− 52	44
2	Cerebellum	3971	7.37	20	− 72	− 36
	L. caudate	326	4.54	− 12	14	− 2
	R. caudate	145	3.82	14	10	10
	R. insula	137	3.75	38	− 18	6
	L. insula	122	3.65	− 36	− 14	8
3	R. lateral occipital cortexR. occipital poleR. angular gyrus	3802	5.53	46	− 74	26
	L. lateral occipital cortexL. occipital pole	432	4.86	− 24	− 88	32
	L. postcentral gyrusL. precentral gyrusL. superior temporal gyrus	264	4.43	− 64	− 10	16
4	L. temporal poleL. frontal orbital cortex	1649	6.52	− 32	14	− 26
	R. temporal poleR. frontal orbital cortex	1340	6.08	36	12	− 24
	R. cuneusR. occipital poleR. lateral occipital cortex	978	5.55	12	− 82	24
	L. cuneusL. lateral occipital cortex	554	5.08	− 10	80	24
	L/R. subcallosal cortex	242	3.89	0	10	− 6
	L/R. anterior cingulate gyrus	212	3.83	0	38	0
	R. frontal pole	78	3.38	34	52	22
	L. frontal pole	74	3.32	− 12	58	20
5	L./R. lateral occipital cortexL./R. precuneusL./R. cuneus	6948	6.97	26	− 74	52
	L./R. lingual gyrusCerebellum	2413	5.71	− 10	− 82	− 16
6	Cerebellum	2412	6.05	8	− 46	− 18
	L./R. supplementary motor area	331	4.97	0	− 10	56
	R. putamen	26	3.19	18	6	− 12
	L. putamen	24	3.18	− 20	10	− 12
7	L. cuneusL. lingual gyrusL. posterior cingulate gyrus	416	5.34	− 16	− 54	4
	L. middle temporal gyrusL. inferior temporal gyrus	300	4.33	− 64	− 12	− 22
	L./R. occipital pole	236	5.34	0	− 92	28
	L. inferior temporal gyrus	207	4.91	− 44	− 28	− 24
	L. hippocampus	141	3.85	− 20	− 14	− 20
	L. hippocampusL. thalamus	102	4.09	− 20	− 38	0
	R. inferior frontal gyrus	99	3.72	52	18	26
	L./R. anterior cingulate gyrusL./R. supplemental motor area	91	3.41	− 2	− 6	46
	R. caudate	72	3.66	10	2	12
	L. caudate	55	3.34	− 8	16	4
	L./R. occipital pole	52	4.20	− 2	− 92	26
8	L./R. intracalcarine cortexL.R. supracalcarine cortexL./R. precuneus	3987	6.21	− 22	− 66	6
	L./R. occipital pole	201	4.36	16	− 100	− 4
	R. frontal pole	87	3.82	20	62	− 2

Legend: L. – Left; R. – Right.

Cluster extent refers to the number of contiguous 2mm3 voxels. X, Y and Z refer to standard Montreal Neurological Institute coordinates.

Fig. 2

Source-based morphometry results showing the eight patterns of gray matter (GM) volume that were identified in all 152 multiple sclerosis patients and 31 healthy controls at baseline. Z-scores, thresholded at | 3 |, are shown in red-yellow for each component, with yellow being associated with larger values. Compared to healthy controls, multiple sclerosis patients had significantly reduced GM volume in Patterns 1 and 7 after correction for multiple comparisons. The X, Y and Z values refer to the slice shown in Montreal Neurological Institute coordinates.

VBM results at baseline and changes over follow-up

No significant differences were detected between patients without CDP and those with CDP at baseline in the VBM analysis. With respect to HCs at baseline, the MS cohort as a whole presented with decreased GM volume within the deep GM as well as in the cortex, albeit to a smaller degree (Supplemental Table 1 and Supplementary Fig. 1). Results of the longitudinal VBM analysis are shown in detail in Supplemental Table 2 and Supplementary Fig. 2. In brief, patients with CDP showed a greater rate of GM atrophy in most of the deep GM structures, cerebellum as well as several cortical regions. Although there was overlap with several of the regions identified in the SBM analysis, the overall pattern was more restricted.

Supplementary Fig. 1

Voxel-based morphometry results comparing gray matter differences between multiple sclerosis patients and healthy controls at baseline. Areas of significantly decreased gray matter volume are shown in red-yellow, with yellow being associated with smaller p-values. Effects related to age, sex, and head size were regressed out from the analyses and results are shown with threshold-free cluster enhancement at p < 0.05, family-wise error corrected. The X, Y and Z values refer to the slice shown in Montreal Neurological Institute coordinates.

Supplementary Fig. 2

Voxel-based morphometry results comparing gray matter atrophy in patients with confirmed disability progression (CDP) and those without after ten years of follow-up. Areas with a significantly greater rate of gray matter loss are shown in red-yellow, with yellow being associated with smaller p-values. Effects related to age, sex, and head size were regressed out from the analyses and results are shown with threshold-free cluster enhancement at p < 0.05, family-wise error corrected. The X, Y and Z values refer to the slice shown in Montreal Neurological Institute coordinates.

Discussion

In this study, we used a data-driven approach for identifying common spatial patterns of GM density in a cohort of relapsing-remitting MS patients. Within the identified patterns, we found that GM atrophy progresses in a manner which follows brain regions that are both functionally and/or structurally related to each other. These findings extend those from a previous report involving the same cohort where it was shown that patients with CDP demonstrated greater whole brain and total GM tissue loss compared to those without CDP (Zivadinov et al., 2016). The SBM approach allowed us to identify spatial patterns of GM volume common to the cohort of MS patients in this study. Although SBM has been suggested to be more sensitive in detecting GM atrophy than other methods such as VBM (Ciarochi et al., 2016, Kasparek et al., 2010), we did not detect differences at baseline between patients split on CDP status at the end of the study. Nonetheless, the CDP patient group experienced a significantly greater rate of GM atrophy in two of the identified patterns. Interestingly, several of the corresponding regions were significantly reduced in volume for the MS cohort as a whole when compared to the HC group at baseline (e.g. insula, hippocampus, thalamus, supplemental motor area, anterior and posterior cingulate gyri). Taken together, these results suggest that already atrophic structures may be at risk for an even greater rate of tissue loss in patients who progress in terms of clinical disability.

The present study is one of the first to utilize SBM in a longitudinal setting. We implemented an optimized pipeline which included the generation of individual subject templates that were subsequently used to create a high-quality group template from all 152 subjects. We then identified 8 patterns of GM volume co-variation from all of the available imaging data. In general, the identified patterns consisted of bilateral GM regions which are known to be either structurally and/or functionally related. For example, Pattern 1 consisted of areas known to be involved with motor control, including deep GM structures (thalamus, caudate and nucleus accumbens), cortical regions (precentral gyrus, supplementary motor area and cingulate cortex) and the cerebellum. Given that EDSS is heavily weighted towards motor dysfunction, it is perhaps relatively unsurprising that the CDP group of patients presents with significantly increased atrophy within these areas over the course of 10 years of follow-up. Interestingly, Pattern 7 consisted of GM areas associated with both learning and memory storage/retrieval (cerebellum, hippocampus, nucleus accumbens, putamen, parahippocampal gyrus, and lingual gyrus). The CDP group developed significantly more GM atrophy in these areas over the follow-up as well. Unfortunately, detailed neuropsychological assessments were not available in this cohort. However, clinical disability, as measured with EDSS, has been linked to cognitive decline in a number of studies (Caneda and Vecino, 2016, Deloire et al., 2010, Moccia et al., 2016, Uher et al., 2017), lending further credence to our results.

Univariate imaging methods, such as VBM, have enjoyed widespread adoption within the MS neuroimaging literature (Eshaghi et al., 2014, Hofstetter et al., 2014). Although VBM is a conceptually simple approach to analyzing group data, our results suggest that it may lack a certain degree of sensitivity to GM atrophy, especially in a longitudinal setting. Similar discrepancies between SBM and VBM have been reported in cross-sectional studies of other neurological diseases such as schizophrenia (Kasparek et al., 2010) and Huntington's disease (Ciarochi et al., 2016). In the current study, Pattern 1, which underwent significantly greater atrophy in the CDP group, contained, in part, the anterior and posterior cingulate gyri, which were not identified by the VBM analysis. Damage to the cingulate cortex has been implicated in a number of a pathological features of MS, including depression (Bonavita et al., 2017), cognitive impairment (Steenwijk et al., 2016), and clinical disability as measured by EDSS (Rudko et al., 2016). Thus, that VBM fails to show any significant differences in the cingulate cortex appears to be a false negative rather than a false positive in the SBM analysis. In fact, it has been proposed that multivariate techniques such as SBM may provide advantages not offered by univariate methods (Kasparek et al., 2010). For example, the multiple comparison problem inherent in univariate methods is reduced given that statistical analysis is done on a summary measure as opposed to across all voxels. Moreover, the multivariate nature of the analysis inherently allows for the identification of underlying relationships between regions with similar characteristics (e.g. GM volume). This can potentially aid in teasing out the multi-faceted effects of GM atrophy. Indeed, we were able to identify a pattern which corresponded primarily to areas involved in motor control. A univariate method, such as VBM, is inherently limited in this regard as voxels are tested individually and then significant clusters are determined after the fact. In addition, SBM can potentially provide additional insight into what drives the total amount of observed GM atrophy in a way that is not possible with VBM. The ICA-derived loading factors in an SBM analysis relate to the GM volume within the corresponding pattern. Thus, one could theoretically investigate whether certain patterns of GM volume contribute more heavily than others to global GM atrophy.

Although we have shown that GM atrophy progresses in a manner which follows structurally and functionally related areas, we cannot conclusively state why this is so. It is conceivable that related GM areas may atrophy in relative unison due to mechanisms such as retrograde and anterograde degeneration following axonal transection (Dalton et al., 2004, Sepulcre et al., 2009). This does not appear to be the only explanation though given that the two groups had similar T2-LVs at baseline and experienced similar absolute changes over the follow-up. In addition, the overall results did not change when T2-LV was added as an additional covariate in the statistical analysis. Thus, it seems likely that other pathophysiological processes are also at work, including primary mechanisms such as meningeal inflammation (Zivadinov et al., 2017), which has been shown to be associated with cortical atrophy. Other hypotheses include selective vulnerability of certain types of neurons to particular aspects of the disease (Geurts and Barkhof, 2008). Nevertheless, the exact mechanisms remain to be elucidated over the coming years.

Our study is not without limitations. First, patients were divided into two groups based solely on whether or not they had CDP after 10 years of follow-up. Utilizing other criteria, such as cognitive impairment status, for defining the sub-groups may have shown increased atrophy in some of the patterns that were not significantly different when comparing patients with CDP to those without. Moreover, we only included subjects who had available MRI as well as clinical data at baseline and 10 years. However, the original protocol included yearly MRI examinations and > 85% of the enrolled subjects had MRI data available any time point over the duration of the entire study (Zivadinov et al., 2016). Thus, it would be highly desirable to exploit all of the available imaging data in order to gain a more precise understanding of the manner in which GM network atrophy evolves over time. One approach would be to utilize the “dual regression” (Filippini et al., 2009) technique which has been widely utilized in the resting state fMRI literature. Implementing such an analysis, however, would require careful consideration with respect to how to handle missing data. We will address this aspect in a future manuscript. Regarding spatial normalization, like VBM, SBM results are influenced by the registration quality. A surface-based strategy would likely have yielded better inter-subject alignment (Steenwijk et al., 2016) than could be achieved using volumetric normalization, as was done in the current study. This stems from the fact that the geometry of the GM/WM boundary can be used to drive the registration with a surface-based technique, rendering alignment quality essentially invariant to WM atrophy. In the current study, it might be the case that the observed GM patterns are not entirely independent from WM volume loss. Nonetheless, a state-of-the-art normalization algorithm was used and all registrations were deemed sufficient during quality control. In addition, some discussion on the implementation of the SBM pipeline is warranted. The ICA algorithm used likely has an impact on the identified components. Although a number of studies have utilized the Infomax algorithm (Kasparek et al., 2010, Steenwijk et al., 2016, Xu et al., 2009), FastICA (as used in the current study) has also been used (Coppen et al., 2016, Hafkemeijer et al., 2014). Moreover, we chose, a priori, to extract 8 components based on other SBM studies which used a similar number (Coppen et al., 2016, Hafkemeijer et al., 2014, Steenwijk et al., 2016). Although automated methods exist for selecting the number of components, they are not consistently used in SBM studies nor it is not clear that this approach results in the optimal partitioning of the data. In addition, more accurate segmentation of the deep GM structures can be obtained by using other model-based approaches (Patenaude et al., 2011). However, the integration of such techniques into the current SBM pipeline would have its own challenges, especially since partial volume estimates at structural borders are generally not generated. Finally, the lack of healthy control data at 10-year follow-up prevented us from assessing whether particular GM networks are more prone than others to atrophy even in patients that do not go on to develop disability progression, as measured by EDSS. Such information would have been of particular interest given that the patient group without CDP experienced an annualized rate of whole brain atrophy that exceeded the proposed pathological cut-off of − 0.40% (De Stefano et al., 2016).

In conclusion, the SBM approach employed in the current study offers a means to identify patterns of GM tissue loss in a manner not afforded by traditional region of interest or univariate methods. The results from this analysis represent an additional step in characterizing the role that GM atrophy plays in the progression of clinical disability over mid-term follow-up in early relapsing-remitting MS patients.

The following are the supplementary data related to this article.

Supplementary tables

Supplementary Video 1

The spatially normalized gray matter segmentations for all 152 MS patients at baseline and 10-year follow-up. Images are shown prior to Jacobian modulation.

Disclosures

N. Bergsland reports no disclosures.

D. Horakova received compensation for travel, speaker honoraria and consultant fees from Biogen Idec, Novartis, Merck Serono, Bayer Shering, and Teva, as well as support for research activities from Biogen Idec.

M.G. Dwyer has received consultant fees from Claret Medical and EMD Serono and research grant support from Novartis.

T. Uher received financial support for conference travel and honoraria from Biogen Idec.

M. Vaneckova received financial support for research activities from Biogen Idec.

M. Tyblova received financial support for research activities from Biogen Idec.

Z. Seidl received financial support for research activities from Biogen Idec.

J. Krasensky received financial support for research activities from Biogen Idec.

E Havrdova received speaker honoraria and consultant fees from Biogen Idec, Merck Serono, Novartis, Genzyme and Teva, as well as support for research activities from Biogen Idec and Merck Serono.

R Zivadinov received personal compensation from EMD Serono, Genzyme-Sanofi, Novartis, Claret-Medical, Celgene for speaking and consultant fees. He received financial support for research activities from Claret Medical, Genzyme-Sanofi, QuintilesIMS Health, Intekrin-Coherus, Novartis and Intekrin-Coherus.