Longitudinal changes of deep gray matter shape in multiple sclerosis.

OBJECTIVE: This study aimed to investigate longitudinal deep gray matter (DGM) shape changes and their relationship with measures of clinical disability and white matter lesion-load in a large multiple sclerosis (MS) cohort.
MATERIALS AND METHODS: A total of 230 MS patients (179 relapsing-remitting, 51 secondary progressive; baseline age 44.5 ± 11.3 years; baseline disease duration 12.99 ± 9.18) underwent annual clinical and MRI examinations over a maximum of 6 years (mean 4.32 ± 2.07 years). The DGM structures were segmented on the T1-weighted images using the "Multiple Automatically Generated Templates" brain algorithm. White matter lesion-load was measured on T2-weighted MRI. Clinical examination included the expanded disability status scale, 9-hole peg test, timed 25-foot walk test, symbol digit modalities test and paced auditory serial addition test. Vertex-wise longitudinal analysis of DGM shapes was performed using linear mixed effect models and evaluated the association between average/temporal changes of DGM shapes with average/temporal changes of clinical measurements, respectively.
RESULTS: A significant shrinkage over time of the bilateral ventrolateral pallidal and the left posterolateral striatal surface was observed, whereas no significant shape changes over time were observed at the bilateral thalamic and right striatal surfaces. Higher average lesion-load was associated with an average inwards displacement of the global thalamic surface with relative sparing on the posterior side (slight left-side predominance), the antero-dorso-lateral striatal surfaces bilaterally (symmetric on both sides) and the antero-lateral pallidal surface (left-side predominance). There was also an association between shrinkage of large lateral DGM surfaces with higher clinical motor and cognitive disease severity. However, there was no correlation between any DGM shape changes over time and measurements of clinical progression or lesion-load changes over time.
CONCLUSIONS: This study showed specific shape change of DGM structures occurring over time in relapse-onset MS. Although these shape changes over time were not associated with disease progression, we demonstrated a link between DGM shape and the patients' average disease severity as well as white matter lesion-load.

Introduction

Multiple sclerosis (MS) is the most common chronic inflammatory, demyelinating and neurodegenerative disease of the CNS in young adults.(Wallin et al., 2019) Among other CNS structures, MS critically affects deep gray matter (DGM) structures such as the thalamus, putamen and globus pallidus, which participate in functionally segregated circuits with motor and non-motor areas of the cerebral cortex and the cerebellum and significantly contribute not only to the control of movement, but also to a variety of cognitive and affective functions.(Bostan et al., 2018, Herrero et al., 2002, Hwang et al., 2017, Lanciego et al., 2012) DGM volume loss occurs in all MS disease types, as well as both in the early and late stages of the disease.(Azevedo et al., 2018, Azevedo et al., 2015, Eshaghi et al., 2018, Magon et al., 2020, Zivadinov et al., 2013) This process most likely reflects the aftermath of both focal demyelination and diffuse neurodegeneration, occurring independently within the DGM or secondary due to retrograde or anterograde myelinoaxonal degradation after white matter injury elsewhere.(Calabrese et al., 2015, Haider et al., 2014, Klaver et al., 2013, Papadopoulou et al., 2019, Vercellino et al., 2009) Moreover, DGM volume loss in particular has been strongly associated with progression of physical and cognitive disability.(Eshaghi et al., 2018, Houtchens et al., 2007, Magon et al., 2014a, Magon et al., 2020, Rocca et al., 2010, Schoonheim et al., 2012, Schoonheim et al., 2015).

So far, DGM volume loss has been evaluated in previous cross-sectional and longitudinal volumetric studies of MS patients.(Eshaghi et al., 2018, Houtchens et al., 2007, Magon et al., 2020) However, analysis of DGM shape changes could provide additional information regarding the exact anatomical location of DGM volume loss and accurately detect the spatiotemporal correlates of tissue injury and disability progression. Several previous studies have shown different trajectories in the progression of atrophy across different CNS regions, suggesting a spatiotemporal heterogeneity of volume loss in MS.(Eshaghi et al., 2018, Tsagkas et al., 2018a, Tsagkas et al., 2021b) In addition, it has been shown that atrophy does not homogeneously progress even within individual CNS structures.(Magon et al., 2020, Tsagkas et al., 2020) Furthermore, previous studies showed that shape analysis is often able to provide spatially localized group differences where volumetric differences are not apparent. This also suggests that shape analyses are sensitive to focal atrophic changes that may precede bulk and detectable volume loss.(Chakravarty et al., 2015, Copersino et al., 2020, Guimond et al., 2021) Thus, shape analysis studies of atrophy progression might be able to capture this within-region heterogeneity and to pinpoint the exact location, where neurodegeneration occurs. This is especially important in DGM structures because of the high functional organization of neuronal populations within the thalamus and basal ganglia (Bostan et al., 2018, Herrero et al., 2002, Nambu, 2011). To that end, a previous cross-sectional study confirmed a correlation between thalamic shape and the expanded disability status scale (EDSS) as well as a relationship between white matter lesion volume and the shape of all DGM structures.(Magon et al., 2014a) However, to our knowledge longitudinal studies examining DGM shape changes in MS patients are currently lacking.

In a previous longitudinal volumetric study, we identified the ventral lateral nucleus as an important driver of disability progression in MS.(Magon et al., 2020) In contrast to these results, other studies suggested that a gradient of thalamic injury at the ventricular/thalamic borders takes place in pediatric MS patients.(De Meo et al., 2021, Fadda et al., 2019) Therefore, in the present explorative study, we aimed to investigate longitudinal DGM shape changes in MS, in order to pinpoint the exact regions of individual DGM structures, and especially the thalamus, in which atrophy —as a marker of neurodegeneration— progresses faster. Moreover, we aimed at examining whether these shape changes are driven by white matter lesions in order to shed light to the mechanisms leading to DGM atrophy in MS. Last, we hypothesized that specific shape changes of DGM structures —such as the thalamus— relate to disability and disability progression over a 6-year follow-up span.

Materials and methods

Study design and participants

We analyzed clinical and MRI data of relapsing-remitting (RRMS) and secondary progressive MS (SPMS) patients from a large scale cohort study (GeneMSA) at a single center (Multiple Sclerosis Center, University Hospital, Basel, Switzerland), previously reported in former studies.(Magon et al., 2020, Tsagkas et al., 2021b, Tsagkas et al., 2021a, Tsagkas et al., 2020, Tsagkas et al., 2018b, Tsagkas et al., 2018a) Patients were followed over a maximum of 6 years (7 annual time points). In cases, where clinical relapses or treatment with steroids occurred, both clinical assessments and MRI were performed with at least one month distance to the event. The diagnosis of MS was made in accordance with international panel established criteria (McDonald et al., 2001). The local ethics committee approved the study (EKBB-46/04) and all patients provided written informed consent.

Clinical assessment

All patients underwent a standardized neurological examination including the EDSS () performed by trained and certified examiners, 9-hole peg Test (9HPT), Timed 25‐foot walk test (T25fwt) and Paced Auditory Serial Addition Test (PASAT) annually. Patients also underwent an annual Symbol Digit Modalities Test (SDMT) starting at the fourth follow‐up time (third year of monitoring). No parallel test versions were used for the SDMT, whereas two versions of the PASAT were deployed for annual neuropsychological tests. All clinical and neuropsychological metrics were recorded at each follow‐up.

MRI protocol

All annual MRI were performed on one 1.5 Tesla Magnetom Avanto MR scanner (Siemens Healthineers, Erlangen, Germany). The protocol included high-resolution three-dimensional sagittal T1-weighted (T1w) MPRAGE (TR/TI/TE = 2080/1100/3.0 ms, α = 15°, 160 slices, voxel size = 0.98x0.98x1mm3). In addition, a double echo spin echo sequence was performed (TR/TE1/TE2) = 3980/14/108 ms, flip angle = 180°, 40 slices, slice thickness = 3 mm, no gap, in-plane resolution = 1x1 mm2) to obtain proton density weighted (PDw) and T2-weighted (T2w) images.

DGM segmentation and shape analysis

The DGM structures were segmented on the T1w images using the “Multiple Automatically Generated Templates” (MAGeT) brain algorithm. The algorithm uses a single atlas derived from manually segmented serial histological data including the thalamus, striatum and globus pallidus.(Chakravarty et al., 2006) MAGeT provides a vertex-wise analysis of DGM shape with 3016 and 3108 vertexes for left and right thalamus, 6450 and 6178 vertexes for left and right striatum, as well as 1266 and 1138 vertexes for left and right globus pallidus respectively. Thus, this method is able to provide a detailed depictions of the surface of DGM structures.

The T1w images were preprocessed using both the medical imaging netCDF (MINC) toolkit (version 2) and the advanced normalization tools (ANTs) as following: (1) bias field correction using N4-correction algorithm;(Tustison et al., 2010) (2) non-local means denoising;(Manjón et al., 2010) (3) affine registration using a normalized mutual information objective function;(Studholme et al., 2001) and (4) brain extraction using the BEaST algorithm.(Eskildsen et al., 2012) The preprocessed data were then used as input for the MAGeT algorithm.

As a first step, the MAGeT brain algorithm customizes the atlas to a subset of subjects under evaluation that are distributed with respect to age, sex, and disease category to adequately capture neuroanatomical variability within the sample using a region‐of‐interest based nonlinear registration scheme.(Chakravarty et al., 2009, Chakravarty et al., 2008) This segmented subset of subjects acts as a template library for the remainder of the dataset. This has the benefit of averaging different sources of random error prior to the estimation of the final segmentation. Once each subject is matched to each of these templates, there are numerous candidate segmentations that are fused using a voxel‐wise majority voting (i.e., the label occurring most frequently at a specific location is retained;(Magon et al., 2014a)). Specifically, in this study the template library included 21 females and 10 males with a mean age of 47.2 years ± 11.9 (range 19–72), a mean disease duration of 17.2 years ± 11.5 (range 1–47) and a median EDSS score of 3 (range 0–7.5).

To extract the shape of each DGM structure, an extension of the surface-based methodology was performed, as proposed previously.(Lerch et al., 2008). In a first step, surface-based representations of each DGM structure for each subject were created. Then the 31 nonlinear transformations mapped to the original surface template were concatenated and averaged to increase accuracy and minimize noise errors.(Dorr et al., 2008, Frey et al., 2011) To provide a measurement of local surface displacement, the dot product between the nonlinear transformation and the surface model normals were computed. Inwards and outwards displacements of the structure’s surface were estimated respectively to the reference atlas (Chakravarty et al., 2006) and are differentiated by a positive or negative sign. Hence displacements in the inwards direction indicate that the segmented DGM structure is smaller than the reference atlas, while outwards displacements indicate a larger segmented structure respectively to the template atlas.

White matter lesions were segmented by trained expert observers at the Medical Image Analysis Center (MIAC AG) using Amira (v 3.1.1 by Mercury Computer Systems) and followed standard operating procedures for the analysis of clinical phase II and phase III trials.(Magon et al., 2014b) In this process, lesions were manually marked and subsequently segmented using intensity thresholding with Amira 3.1.1 (Mercury Computer System) with consecutive manual editing. Lesion quantification was performed on PDw with T2w data displayed simultaneously, which allowed for confirmation of the lesion site and extent. All raters underwent standardized training before working on any study, whereas inter- and intra-rater reliability was routinely assessed and verified annually according to the institutional standard operation procedures. Intra-rater and inter-rater reproducibility were high with a mean relative error of <5% und intraclass correlation coefficient >0.99. Finally, each segmentation was validated by a board-certified neuroradiologist as a final quality verification. MAGeT was performed on T1w images after lesion filling using the method previously developed from our group to reduce bias related to tissue misclassification (Magon et al., 2014b) and improve the registration step.(Sdika and Pelletier, 2009).

Ventricle volume segmentation

In MS, as a result from disease-related processes such as both progressive brain volume loss and neuroinflammation, ventricular volume exhibits dynamic changes over time.(Fleischer et al., 2021, Millward et al., 2020, Sinnecker et al., 2020) Since ventricular size influences morphometric DGM measures as well, a correction of shape changes in our statistical analysis using ventricular volumes at each time point was necessary. Therefore, we segmented the brain lateral ventricles using Automatic Lateral Ventricle delineation (ALVIN).(Kempton et al., 2011) ALVIN is a fully automated segmentation-based algorithm to measure the brain lateral ventricle volume normalized to the head size, which has been successfully used in previous longitudinal studies of MS patients.(Sinnecker et al., 2020) This method provides measurement of lateral ventricle volume, which refers to the sum lateral ventricle volume of both sides. To prevent misclassification of periventricular T1w hypointense lesions as part of the ventricles, we used only lesion-filled T1w MPRAGE data as input for ALVIN. After segmentation, the results were reviewed by three experienced readers and segmentation errors were excluded from further analysis.

Statistical analysis

Comparisons of demographic factors, clinical measurements, and number of follow-ups between MS subtypes were made using Welch's and Pearson's chi-squared test with Yate's continuity correction. A logarithmic transformation of the EDSS was performed in order to correct for its nonlinearity in representing physical disability.

Vertex‐wise longitudinal analysis was performed using linear mixed effect models (LMER) in order to explore longitudinal association between the patients' DGM surface displacements as dependent variables, and demographic (including comparisons between disease types), clinical as well white matter lesion load measurements as independent variables. This was done using a random intercept and a random time slope for each subject to allow for within-subject and between-subject variance. In our models, the independent variables were entered blockwise keeping the following sequence: first the lateral ventricle volumes, then demographics and then clinical variables or lesion-load respectively. Each variable was tested first for its correlation to the DGM surface displacement intercept and then for its correlation to the DGM surface displacement slope over time. For instance, for comparisons between MS disease types, the correlation to the DGM surface displacement intercept represents the average between-group DGM shape difference using all available MRI data of each patient; this factor does not take changes of DGM shape occurring over our observational period into consideration. In contrast, in this example, the correlation of MS disease type to the DGM surface displacement slope over time represents the between-group differences in DGM shape changes occurring over our observational period using all available MRI data of each patient. Similarly, both the correlations between average DGM shape measurements and average clinical measurements (e.g. EDSS), as well as between DGM shape changes over time and changes of clinical measurements over time were evaluated.(Diggle et al., 2013) All independent variables without statistical significance with an alpha (p-value) set at 0.05 were excluded from the final model. Furthermore, since 3016 and 3108 vertexes were analyzed for left and right thalamus, 6450 and 6178 vertexes for left and right striatum, as well as 1266 and 1138 vertexes for left and right globus pallidus had to be corrected our results for multiple comparisons. This was done in order to reduce the risk of type I errors. Hence, the results were corrected for multiple comparisons by using the False Discovery Rate approach with a threshold (q-value) of 0.05 representing the rate that features called significant truly reject the null hypothesis. Thus, all of the vertexes reported in our results fulfil the criterion of p < 0.05 as well as the criterion of q < 0.05 after multiple comparisons correction. The strength of each correlation is graphically depicted in our figures using t-value maps. Results are reported in the form of mean t‐value (MTV) ± standard deviation and percentage of significant vertices (PSV) of each cluster in the individual DGM structures in our text and tables. The latter was defined as the number of significant vertices in the respective analysis divided by the total number of vertices of the respective DGM structure (3016 and 3108 for left and right thalamus; 6450 and 6178 for left and right striatum; 1266 and 1138 for left and right globus pallidus respectively). Due to the explorative nature of our study with no prespecified key hypothesis during data collection, we did not correct our results for the conduction of multiple tests (e.g. in the between-group comparisons).(Bender and Lange, 2001, Goeman and Solari, 2011).

All statistical analyses of the shape were performed in R (3.5.1 ) using the RMINC package ().

Results

Participants

Due to segmentation failures, 222 MRI datasets (including a total of three datasets, in which ALVIN segmentations of ventricular volumes failed) were excluded from further statistical analysis. A total of 230 patients (179 RRMS, 51 SPMS) and 1158 MRI datasets were entered in the statistical analysis. Patients were followed-up over a maximum of 6 years (mean 4.32 ± 2.07 years). Baseline demographics, clinical and MRI characteristics are described in Table 1.

Table 1

Baseline demographics, clinical and MRI characteristics of all MS patients.

Subjects	All Patients	RRMS	SPMS
Number of MS patients	230	179	51
Age at baseline (years, mean ± sd)	44.5 ± 11.3	41.36 ± 10.22	55.29 ± 7.62
Sex (female/male)	159/71	132/47	27/24
Disease duration at baseline (years, mean ± sd)	12.99 ± 9.18	11.33 ± 8.29	18.96 ± 9.74
EDSS at baseline (range: 0–7.5, mean ± sd)	3.05 ± 1.63	2.57 ± 1.36	4.74 ± 1.38
9HPT – dominant hand at baseline (seconds, mean ± sd)	21.84 ± 6.50	20.45 ± 4.96	26.75 ± 8.64
9HPT – non dominant hand at baseline (seconds, mean ± sd)	24.10 ± 10.12	22.40 ± 7.81	30.04 ± 14.32
T25FWT at baseline (seconds, mean ± sd)	7.50 ± 10.76	5.67 ± 7.09	14.15 ± 17.39
PASAT at baseline (mean ± sd)	43.78 ± 11.65	45.04 ± 11.24	39.28 ± 12.10
SDMT at baseline (mean ± sd)	46.86 ± 13.42	48.42 ± 13.94	41.30 ± 9
T2w lesion volume at baseline (mm3, mean ± sd)	6153.07 ± 6970	5884.79 ± 6818.79	7113.51 ± 7479.45

Abbreviation: sd. = standard deviation.

DGM shape changes over time

DGM shape changes over time in our MS cohort are displayed in Fig. 1. A significant inwards displacement over time at the ventrolateral pallidal surface bilaterally was observed (MTV/PSV left: −3.17 ± 0.34/27.9%; right: −2.96 ± 0.28/23.4%), whereas a significant outwards displacement over time was seen at the anterodorsal pallidal surface bilaterally with a left predominance (MTV/PSV left: 3.00 ± 0.36/13.2%; right: 2.70 ± 0.17/5.5%). Moreover, a significant inwards displacement over time at the left posterolateral striatal surface was observed (MTV/PSV −3.30 ± 0.46/5.0%), whereas a significant outwards displacement over time was seen at the medial striatal surface (MTV/PSV −3.47 ± 0.44/8.6%). No shape changes over time were observed at the bilateral thalamic and right striatal surfaces.

Fig. 1

Shape changes over time in DGM structures. As shown by the t-values extracted from our linear mixed effect models, the gradient from yellow to red indicates a lower to higher positive correlation (outward displacement over time) and the gradient from pink to blue indicates a lower to higher negative correlation (inward displacement over time) respectively. The upper and lower borders of each color bar were defined by the extracted from the t-values in each model, after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Comparisons of DGM shape between MS phenotypes

The association between DGM shape and sex, age, and disease duration are displayed in the Appendix, Supplementary Table 1 and Supplementary Figs. 1-3. Compared to RRMS, SPMS patients had an average inwards shape displacement of a small cluster in the lateral left striatum surface, the right lateral and left antero-dorsal pallidal surface, as well as an average outwards shape displacement of the left ventral striatal surface and ventromedial pallidal surface bilaterally. DGM shape changes over time did not significantly differ between RRMS and SPMS. Comparisons between MS disease types are shown in Table 2, Azevedo et al., 2015 and Fig. 2.

Table 2

Association of DGM shape with disease type (SPMS compared to RRMS).

Variable	Surface Displacement	Left Thalamus	Right Thalamus	Left Striatum	Right Striatum	Left Globus Pallidus	Right Globus Pallidus
DiseaseType(SPMS)	Inwards	n.s.	n.s.	−4.48 ± 0.50(0.17%)	n.s.	−3.03 ± 0.24(14.8%)	−3.33 ± 0.40(24.7%)
	Outwards	n.s.	n.s.	5.10 ± 1.03(7.4%)	n.s.	3.79 ± 0.77(18.7%)	3.61 ± 0.72(26.7%)
Disease Type (SPMS) *Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.

Results represent the correlation strength shown in the form of t‐values extracted from our linear mixed effect models after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. Lower (or more negative) t‐values reflect a stronger negative correlation in inwards surface displacements, whereas higher (or more positive) t-values reflect a stronger positive correlation in outwards surface displacements. All results are reported in the form: mean t‐value ± SD, (percentage of significant vertices in the respective DGM structure).

Abbreviation: n.s. = not significant.

Fig. 2

Association of average shape of DGM structures associated with SPMS, as compared to RRMS. As shown by the t-values extracted from our linear mixed effect models, the gradient from yellow to red indicates a lower to higher positive correlation (outward displacement in SPMS) and the gradient from pink to blue indicates a lower to higher negative correlation (inward displacement in SPMS) respectively. The upper and lower borders of each color bar were defined by the extracted from the t-values in each model, after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Association of DGM shape with T2w lesion-load

The association between DGM shape and T2w lesion load in our MS cohort are displayed in Table 2 and Fig. 3. Differences of these associations between RRMS and SPMS groups are shown in Supplementary Table 2. Rough anatomical localizations of these associations are reported in Table 4.

Fig. 3

Association between the average shape of DMG structures and lesion load. As shown by the t-values extracted from our linear mixed effect models, the gradient from yellow to red indicates a lower to higher positive correlation (outward displacement with higher lesion load) and the gradient from pink to blue indicates a lower to higher negative correlation (inward displacement with higher lesion load) respectively. The upper and lower borders of each color bar were defined by the extracted from the t-values in each model, after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 4

Association between the shape of DGM structures and all clinical outcomes.

Variable	Surface Displacement	Left Thalamus	Right Thalamus	Left Striatum	Right Striatum	Left Globus Pallidus	Right Globus Pallidus
logEDSS	Inwards	−3.89 ± 1.04(49.6 %)	−3.89 ± 0.83(40.4 %)	−4.07 ± 1.02(34.7 %)	−3.60 ± 0.71(28.9 %)	−3.62 ± 0.81(23.1 %)	−3.71 ± 0.72(21.5 %)
	Outwards	3.33 ± 0.77(5.90 %)	2.78 ± 0.38(1.06 %)	3.25 ± 0.71(17.3 %)	3.27 ± 0.59(17.0 %)	3.60 ± 0.63(23.0 %)	3.54 ± 0.63(21.2 %)
logEDSS*Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
T25fwt	Inwards	−3.63 ± 0.91(33.3 %)	−3.61 ± 0.90(35.2 %)	−3.72 ± 0.86(24.6 %)	−4.40 ± 1.17(25.6 %)	−3.34 ± 0.41(27.7 %)	−3.69 ± 1.25(21.7 %)
	Outwards	3.28 ± 0.58(24.3 %)	3.63 ± 0.78(25.1 %)	3.57 ± 0.69(22.0 %)	4.21 ± 1.17(28.8 %)	3.60 ± 0.61(32.4 %)	3.73 ± 1.16(26.6 %)
T25fwt*Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
SDMT	Inwards	(only SPMS)−3.73 ± 0.23(2.18 %)	n.s.	−3.19 ± 0.43(14.6 %)	−3.33 ± 0.50(14.4 %)	−3.58 ± 0.57(21.0 %)	−3.53 ± 0.53(24.1 %)
	Outwards	(only SPMS)3.51 ± 0.13(0.76 %)	n.s.	3.43 ± 0.48(32.4 %)	3.32 ± 0.42(31.6 %)	2.93 ± 0.26(26.6 %)	3.03 ± 0.36(17.9 %)
SDMT*Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
9HPT	Inwards	−3.51 ± 0.85(38.8 %)	−2.95 ± 0.38(36.1 %)	−4.18 ± 1.41(27.6 %)	−3.07 ± 0.22(8.38 %)	−4.23 ± 1.10(32.2 %)	n.s.
	Outwards	2.52 ± 0.17(1.06 %)	n.s.	3.70 ± 0.79(15.2 %)	3.12 ± 0.26(4.80 %)	3.57 ± 0.62(19.1 %)	n.s.
9HPT*Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
PASAT	Inwards	n.s.	(only SPMS)−3.51 ± 0.50(8.59 %)	−4.20 ± 1.03(23.7 %)	−3.60 ± 0.65(9.22 %)	−4.67 ± 1.06(32.1 %)	−3.28 ± 0.39(14.4 %)
	Outwards	n.s.	(only SPMS)4.61 ± 1.19(26.0 %)	4.27 ± 1.10(36.4 %)	3.88 ± 0.81(17.2 %)	4.0 ± 0.82(36.6 %)	2.73 ± 0.02(0.17 %)
PASAT*Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.

Results represent the correlation strength shown in the form of t‐values extracted from our linear mixed effect models after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. Lower (or more negative) t‐values reflect a stronger negative correlation in inwards surface displacements, whereas higher (or more positive) t-values reflect a stronger positive correlation in outwards surface displacements. All results are reported in the form: mean t‐value ± SD, (percentage of significant vertices in the respective DGM structure).

Abbreviation: n.s. = not significant.

Higher average lesion-load was associated with an average inwards displacement of i) the global thalamic surface with relative sparing on the posterior side (slight left-side predominance), ii) the antero-dorso-lateral striatal surfaces bilaterally (symmetric on both sides) and iii) the antero-lateral pallidal surface (left-side predominance). Higher average lesion-load was also associated with an average outwards displacement of i) small clusters of the left ventral thalamic and right posterior thalamic surfaces, ii) the bilateral ventromedial striatal surfaces (left-side predominance) and iii) the medial pallidal surfaces. Lesion-load changes over time were not associated with DGM shape changes over time.

Association of DGM shape with clinical outcomes

The association between DGM shape and clinical outcomes in our MS cohort are displayed in Table 4 and in Fig. 3, Fig. 4. Differences of these associations between RRMS and SPMS groups are shown in Supplementary Table 2. Rough anatomical localizations of these associations are reported in Supplementary Table 3.Fig 5.

Fig. 4

Association of average shape of DGM structures associated with motor clinical outcomes. As shown by the t-values extracted from our linear mixed effect models, the gradient from yellow to red indicates a lower to higher positive correlation (outward displacement with higher EDSS, 9HPT and T25fwt scores) and the gradient from pink to blue indicates a lower to higher negative correlation (inward displacement with higher EDSS, 9HPT and T25fwt scores) respectively. The upper and lower borders of each color bar were defined by the extracted from the t-values in each model, after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. Upper row: association between the average shape of DGM structures and the EDSS. Middle row: association between the average shape of DGM structures and the T25fwt. Lower row: association between the average shape of DGM structures and the 9HPT. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

Association of average shape of DGM structures associated with cognitive clinical outcomes. As shown by the t-values extracted from our linear mixed effect models, the gradient from yellow to red indicates a lower to higher positive correlation (inward displacement with lower SDMT and PASAT scores) and the gradient from pink to blue indicates a lower to higher negative correlation (outward displacement with lower SDMT and PASAT scores) respectively. The upper and lower borders of each color bar were defined by the extracted from the t-values in each model, after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. Upper row: association between the average shape of DGM structures and the SDMT. Lower row: association between the average shape of DGM structures and the PASAT. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Higher average EDSS was associated with an average inwards displacement of i) the antero-lateral thalamic surface bilaterally with left predominance, ii) dorsolateral striatal surface bilaterally, and iii) the ventro-antero-lateral pallidal surface bilaterally. In addition, higher average EDSS was associated with an average outwards displacement of i) small clusters of the posteriolateral thalamic surface bilaterally, ii) the ventromedial striatal surface bilaterally and iii) dorso-postero-medial pallidal surface bilaterally. EDSS changes over time were not associated with any DGM shape changes over time.

Higher average T25fwt was associated with an average inwards displacement of i) the antero-lateral thalamic surface bilaterally, ii) left anterolateral and right anterior striatal surface and iii) anterolateral pallidal surface bilaterally with left predominance. In addition, higher average T25fwt was associated with an average inwards displacement of i) the posteromedial thalamic surface bilaterally, ii) ventromedial striatal surface bilaterally and iii) medial pallidal surface bilaterally. T25fwt changes over time were not associated with any DGM shape changes over time.

Higher average 9HPT was associated with an average inwards displacement of i) the anterolateral thalamic surface bilaterally, ii) lateral striatal surface bilaterally and iii) left ventrolateral pallidal surface. In addition, higher average 9HPT was associated with an average outwards displacement of i) a small cluster of the left posterolateral thalamic surface, ii) medial striatal surface bilaterally and iii) left dorsomedial pallidal surface. 9HPT changes over time were not associated with any DGM shape changes over time.

Lower average SDMT and PASAT scores were associated with an average inwards displacement of i) the right anterodorsal (only PASAT) and a small cluster of the left medial thalamic surface, ii) bilateral dorsolateral striatal surface, and iii) bilateral lateral pallidal surface (only left for PASAT). In addition, lower average SDMT and PASAT scores were associated with an outwards displacement of i) small clusters in the right ventral (only PASAT) and left lateral (only SDMT) surface, ii) medial striatal surface bilaterally with left predominance and iii) lateral pallidal surface bilaterally (only left for PASAT). SDMT and PASAT changes over time were not associated with any DGM shape changes over time.

Discussion

In this study we investigated the shape of DGM structures in a large longitudinal relapse-onset cohort of MS patients followed-up over 6 years and explored their association with several demographic factors, white matter lesion-load and clinical disability assessments. We observed significant shape changes over time in the globus pallidus and left striatum, which were not influenced by the patients’ sex, age, disease duration and clinical disease phenotype. Similarly, DGM shape changes over time were generally not associated with the increase of MS-related disability over time. However, we found strong correlations between the average shapes of all DGM structures and the severity of motor and cognitive disability.

In MS, the occurrence of volume loss over time in DGM structures – the thalami and the thalamic nuclei in particular – is well established.(Azevedo et al., 2018, Azevedo et al., 2015, Eshaghi et al., 2018, Magon et al., 2020, Zivadinov et al., 2013) However, in our longitudinal analysis, no significant thalamic shape changes were found over time. This finding does not contradict the presence of thalamic atrophy progression in our patients, and is most likely linked to the high between-patient spatial variability of longitudinal shape changes taking place in this DGM structure. Nevertheless, our analysis evaluated the outer thalamic borders, whereas a previous volumetric study of thalamic nuclei of the same cohort indirectly evaluated the intra-thalamic borders of these structures as well.(Magon et al., 2020) This might partly explain the discrepancy between the two approaches and point to the advantages of volumetry for evaluation of atrophy in the thalamic nuclei, since atrophy might also occur within the outer thalamic borders and cannot be investigated with the shape analysis applied in this study. A recent study of pediatric MS patients showed an “in-gradient” of thalamic atrophy progression on the ventricular side,(Fadda et al., 2019) suggesting that the presence of a possible soluble neurotoxic factor has a deleterious effect on this structure. In our study, this finding was not reproduced. In contrast, we observed a significant shrinkage over time of the posterolateral putamen. This putamen territory receives large sensorimotor cortico-striatal projections,(Bostan et al., 2018, Parent and Hazrati, 1995) belongs to a bihemispheric network including the cerebellum (Viñas-Guasch and Wu, 2017) and is involved in sensory, motor and language functions. In parallel, we observed an expansion over time in multiple clusters at the left medial putamen surface, in areas where striatal efferents arise to the internal and external globus pallidus and to the substantia nigra, as part of the cortico-basal ganglia-thalamo-cortical circuits.(Bostan et al., 2018, Lanciego et al., 2012) Furthermore, we found a shrinkage in the ventrolateral pallidum, including regions known for their involvement in motor and motivational-behavioral functions, as well as an expansion in the anterodorsal or “associative” pallidum.(Saga et al., 2017b, Saga et al., 2017a) These findings provide evidence for an MS-related structural reorganization of the basal ganglia with parallel shrinkage and expansion within the same structures over time.

Interestingly, larger lesion-volumes in our patients were associated with an almost global bilateral thalamic shrinkage. The thalamus is not only an important relay station of information from the periphery and subcortical structures to the cortex, but is also involved in the integration of information originating across multiple functional brain networks.(Guillery and Sherman, 2002, Hwang et al., 2017) For this reason, axonal transection of the afferent and efferent thalamic fibers is thought to at least contribute -at least partially- to neurodegeneration occurring in this DGM structure and a cross-sectional and longitudinal association of thalamic atrophy and white matter lesion load has been reported in several previous studies.(Fabiano et al., 2003, Magon et al., 2020, Sepulcre et al., 2006) These results are also in line with a recent study identifying a lesion-led radiologic MS disease phenotype,(Eshaghi et al., 2021) which was also highly associated with pronounced and early DGM atrophy, most probably driven by the thalamus.(Eshaghi et al., 2018) In addition to the thalamus, higher lesion-loads were associated with a shrinkage in large DGM areas receiving efferent fibers from the cortex or from more laterally located DGM structures, such as the anterolateral striatum and the external globus pallidus. At the same time, an expansion of the medial surface of these structures was observed, maybe suggesting compensatory changes within the striatopallidal pathway.(Lanciego et al., 2012).

Moreover, our analyses revealed an association between shrinkage of large lateral DGM surfaces as well as an expansion in smaller medial surface clusters with higher clinical disease severity and larger white matter lesion-loads. Nevertheless, a slight dissociation was apparent between motor related scores, such as the EDSS, T25fwt and 9HPT, and assessments of the cognitive function such as the SDMT and PASAT. Motor-related disability seemed to be more driven by atrophic changes in the anterior parts of DGM structures, whereas cognitive-related disability was driven by more diffuse and posterior located DGM atrophy. This difference is most probably explained with the fact that separate DGM areas are involved with motor and cognitive functions, with motor-related regions mostly located anteriorly.(Helie et al., 2013, Turner and Desmurget, 2010) On the other hand, the significant overlap of DGM areas related with both motor and cognitive areas is expected, since similar DGM regions are associated in motor and cognitive agility or speed through a function of dopamine.(Hanakawa et al., 2017).

Despite the longitudinal nature of our study with annual MRI and clinical measurements over a time span of 6 years, no associations between shape changes over time and clinical worsening or lesion-load increase were found in our analysis. However, a longitudinal correlation between DGM volumes and clinical outcomes has been well established in multiple previous studies, (Eshaghi et al., 2018, Houtchens et al., 2007, Magon et al., 2014a, Rocca et al., 2010, Schoonheim et al., 2015, Schoonheim et al., 2012) including a recent longitudinal volumetric DGM analysis of the same MS cohort.(Magon et al., 2020) DGM shape changes over time essentially reflect subregional volumetric changes of these structures, which should be associated with clinical worsening in a similar way to atrophy. Nevertheless, to that end, a number of methodological aspects have to be taken into account. First, both significant DGM atrophy and global CNS volume take place in MS.(De Stefano et al., 2016) The latter leads to an expansion of the ventricles,(Sinnecker et al., 2020) which possibly changes the spatial relationship between DGM structures and other natural landmarks of the brain. In addition, recent studies also proposed neuroinflammation as a possible driver of dynamic ventricular volume changes over time.(Fleischer et al., 2021, Millward et al., 2020) Hence, adjusting for this effect was necessary in our statistical models. However, it may well be that this statistical manipulation made the detection of subtle changes in the DGM surfaces difficult. On the other hand, it may well be that shape analyses do not have the same sensitivity for the slow degenerative processes occurring in MS, compared to volumetric assessments of DGM structures, which quantify atrophy in a much more straightforward way. Finally, in MS, a between-patient heterogeneity with regard to clinical and radiologic features is well known;(Thompson et al., 2018) hence, it cannot be excluded, that the large spatial variance of longitudinal DGM shape changes hampered the detection of specific surface clusters in our analysis.

One limitation of this study is related to the lack of a control group of healthy subjects, which would have possibly allowed for a more clear differentiation between age-related and MS-specific DGM shape changes. In addition, DGM structures with an established central role for cognitive symptoms in MS(Ciolac et al., 2021, Rocca et al., 2018) were not evaluated in our work. Future studies should consider including these structures, especially when cognitive disability in MS is being the main study focus. Another potential limitation of the present study is the use of a single MRI contrast (T1w) for DGM segmentation. To that end, there is evidence showing that additional T2-weighted MR-contrast or diffusion tensor imaging could contribute to this task.(Traynor et al., 2011) However, although T2-weighted images were acquired in our cohort, they were exclusively 2D with a lower resolution compared to the T1w data, insufficient for the methodology used in this work. Furthermore, diffusion tensor imaging data was not available. Lesions in our study were delineated on PDw images with the according slices of T2w images displayed at the same time, which allowed for confirmation of the lesion site and extent. As the patient recruitment of this longitudinal study first started in June of 2004, T2w 3D techniques like 3D FLAIR were not available at this time. With such 3D techniques, a number of novel methods for fully automatic lesion segmentation have been recently developed.(Andermatt et al., 2018, Egger et al., 2017, Kuijf et al., 2019, La Rosa et al., 2020, Zeng et al., 2020) However, T2w lesion load assessed manually on 2D PDw/T2w is an established imaging approach for lesion segmentation in MS. (Barkhof et al., 2007, Kappos et al., 2015, Kappos et al., 2010, Offenbacher et al., 1993, Radue et al., 2012, Rocca et al., 2013, Sprenger et al., 2020) Further, this study evaluated 1.5T MR-images; future studies could utilize MR-systems of higher magnetic field (e.g. 3T or 7T), which could potentially allow for precise delineation of individual DGM subnuclei (e.g. caudate and putamen). In this study, the potential effect of disease modifying drugs (DMD) on inflammation, and ultimately on DGM shape variation was not included in our analysis. In our cohort 68% of patients were treated with DMD at baseline including primarily first line injectables (63%). While injectables also show an effect on brain atrophy we believe that this effect is rather negligible.(Favaretto et al., 2018) Nevertheless, even a potential therapeutic effect of DMD to DGM shape should not have influenced the association between atrophy and clinical disability. In our explorative statistical analysis, we did not correct our results for the conduction of multiple tests (e.g. in the between-group comparisons).(Bender and Lange, 2001, Goeman and Solari, 2011) Despite the fact, that this practice in still under ongoing debate, this could have led to an inflation of probability of type I error because of the increase in the number of comparisons, in a part of our analysis, specifically the comparisons between the RRMS and SPMS patients. Hence, these results should be considered with caution. Finally, the variability across MRI scans due to methodological or physiological factors also poses a possible limitation. Nevertheless, using the MAGeT segmentation algorithm, we generated a template library using MRI-sessions from all timepoints and patients with different demographic and clinical characteristics in our dataset. In that way, we believe that our template library was representative of our sample and that this approach mitigated possible sources of variability across the input data set.

Conclusions

In conclusion, this study showed specific shape change of DGM structures occurring over time in relapse-onset MS. Although these shape changes were not associated with disease progression, we demonstrated a link between DGM shape and patients’ average disease severity as well as white matter T2w lesion load.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 3

Association between the shape of DGM structures and T2w lesion load.

Variable	Surface Displacement	Left Thalamus	Right Thalamus	Left Striatum	Right Striatum	Left Globus Pallidus	Right Globus Pallidus
T2w lesion load	Inwards	−4.90 ± 1.55(77.4%)	−4.91 ± 1.84(85.7%)	−5.89 ± 2.10(51.4%)	−5.21 ± 1.63(44.5%)	−3.50 ± 0.87(43.9%)	−5.11 ± 1.32(44.5%)
	Outwards	2.62 ± 0.37(0.76%)	2.35 ± 0.25(0.35%)	4.91 ± 1.62(29.2%)	3.86 ± 1.12(26.2%)	4.93 ± 1.63(30.0%)	6.09 ± 1.63(36.0%)
T2w lesion load *Time	Inwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.
	Outwards	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.

Results represent the correlation strength shown in the form of t‐values extracted from our linear mixed effect models after correction with the false discovery rate approach for multiple comparisons set at q < 0.05. Lower (or more negative) t‐values reflect a stronger negative correlation in inwards surface displacements, whereas higher (or more positive) t-values reflect a stronger positive correlation in outwards surface displacements. All results are reported in the form: mean t‐value ± SD, (percentage of significant vertices in the respective DGM structure).

Abbreviation: n.s. = not significant.