Gray matter integrity predicts white matter network reorganization in multiple sclerosis

Abstract Multiple sclerosis (MS) is a chronic inflammatory and neurodegenerative disease leading to gray matter atrophy and brain network reconfiguration as a response to increasing tissue damage. We evaluated whether white matter network reconfiguration appears subsequently to gray matter damage, or whether the gray matter degenerates following alterations in white matter networks. MRI data from 83 patients with clinically isolated syndrome and early relapsing–remitting MS were acquired at two time points with a follow‐up after 1 year. White matter network integrity was assessed based on probabilistic tractography performed on diffusion‐weighted data using graph theoretical analyses. We evaluated gray matter integrity by computing cortical thickness and deep gray matter volume in 94 regions at both time points. The thickness of middle temporal cortex and the volume of deep gray matter regions including thalamus, caudate, putamen, and brain stem showed significant atrophy between baseline and follow‐up. White matter network dynamics, as defined by modularity and distance measure changes over time, were predicted by deep gray matter volume of the atrophying anatomical structures. Initial white matter network properties, on the other hand, did not predict atrophy. Furthermore, gray matter integrity at baseline significantly predicted physical disability at 1‐year follow‐up. In a sub‐analysis, deep gray matter volume was significantly related to cognitive performance at baseline. Hence, we postulate that atrophy of deep gray matter structures drives the adaptation of white matter networks. Moreover, deep gray matter volumes are highly predictive for disability progression and cognitive performance.


| INTRODUCTION
The conception of the human brain as a complex network includes elements that integrate neural information (gray matter), as well as their connecting pathways (white matter). In multiple sclerosis (MS), a chronic inflammatory autoimmune disease, this network is reconfigured due to demyelinating and inflammatory processes, which strongly alter neural processing (Fleischer, Radetz, et al., 2019). Ultimately, consequences include physical disability and cognitive impairment in domains such as processing speed, episodic memory, and executive functions (Johnen et al., 2017).
Advanced imaging techniques allow investigation of both the integrity of white and gray matter in order to evaluate disease progression (Fleischer et al., 2017;Fleischer, Radetz, et al., 2019;Muthuraman et al., 2016). By means of magnetic resonance imaging (MRI), a decrease of cortical thickness has been observed in the posterior cortices, temporal lobe, and precentral gyrus with progressing MS (Eshaghi et al., 2018;Steenwijk et al., 2016). Subcortical volume mainly decreased in deep gray matter structures including thalamus, caudate, and brain stem, as well as cerebellar regions (Deppe et al., 2016;Eshaghi et al., 2018) in early stages of MS. Atrophy of deep gray matter nuclei is strongly related to cognitive impairment in MS (Damjanovic et al., 2017). In particular, thalamic atrophy is predictive of cognitive decline (Papathanasiou et al., 2015;Štecková et al., 2014) and correlates with physical disability (Deppe et al., 2016). Gray matter atrophy, measured with MRI, reflects neuronal and axonal loss (Popescu et al., 2015). However, although addressed by several studies (Bergsland et al., 2015;Bodini et al., 2009;Henry et al., 2009), it is not yet clear whether it is primarily driven by gray matter neurodegeneration, or a consequence of white matter damage (Steenwijk et al., 2015).
Diffusion-weighted imaging allows to evaluate white matter microstructural integrity, as more intact myelination leads to more pronounced directional water flow in nerve tracts (Basser, Mattiello, & LeBihan, 1994), which can be quantified using MRI. Connections between brain regions are inferable with methods such as probabilistic tractography (Basser, Pajevic, Pierpaoli, Duda, & Aldroubi, 2000). The impact of pathological processes on both white, as well as gray matter network configuration and efficiency in MS can then be assessed using graph theory (Fleischer, Radetz, et al., 2019).
Applying this approach, patients with relapsing-remitting MS (RRMS) could be distinguished from those with clinically isolated syndrome (CIS) based on stronger decomposition of their network into distinct modules (Muthuraman et al., 2016). Lower global and local network efficiency within white matter networks of RRMS patients was significantly associated with clinical disability, disease duration, and white matter lesion load (Shu et al., 2011). Further, lower global and local efficiency, as well as clustering, were related to worse attention and information processing capacities (Shu et al., 2016). In a further study, a decrease in clustering and global efficiency was observed in patients with RRMS compared to healthy controls, while the distance measure path length was increased (Llufriu et al., 2017). As a first step to investigate network alterations in relation to disease progression, groups of patients with different disease durations were investigated (Fleischer et al., 2017). Modularity and clustering coefficient appeared to increase in the first year after disease onset and subsequently decrease, suggesting an initial compensatory response that cannot be obtained with disease progression (Fleischer et al., 2017).
In a recent approach, information from gray matter and white matter networks was combined in order to evaluate which MRI parameters best explain cognitive performance (Charalambous et al., 2019). While controlling for age, gender, lesion load, and brain volume, the authors found that global efficiency significantly explained additional variance of processing speed and clinical disability scores.
Firstly, we aimed to identify regions of pronounced gray matter atrophy over the one-year study period. We were then interested in how neurodegenerative processes in these regions relate to white matter network properties of modularization, characteristics relevant for local information processing, and overall network distance change (i.e., emerging disconnection) over time. These measures have previously proven to be highly relevant markers of network functioning in MS and other neurodegenerative disorders (Fleischer et al., 2017;Koirala et al., 2018).
However, it is still an unanswered question whether gray matter damage or white matter network reorganization occurs first. Hence, we examined whether gray matter integrity is predictive of white matter network reconfigurations, or whether white matter network parameters can predict gray matter integrity variations, that is, gray matter atrophy. Lastly, we assessed which of the gray and white matter integrity parameters best predict physical disability and cognitive performance.

| Participants
We analyzed data of 83 patients (58 female, mean age = 35.33 years, SD = 10.84 years) with early RRMS (n = 61) and CIS (n = 22), who were diagnosed according to the revised McDonald criteria 2010 (Polman et al., 2011) from the Clinic of Neurology. CIS and RRMS were classed as early if the disease duration (defined as time between the manifestation of a patient's first symptoms and the date of the first MRI examination) was ≤3 years. Table 1 provides descriptive demographic details of all patients and an overview of the distribution of the patients' ongoing therapies at the time of MRI examinations.
The following exclusion criteria were applied for all patients: any preexisting medical condition known to be associated with brain pathology; pregnancy; previous or current addiction to substances; relapses or systemic therapy with steroids (intravenous, intrathecal, or oral) within the month before the MRI examination; a history of additional neurological or psychiatric disorders. All patients were measured at two time points (mean = 1.12, SD = 0.48 years apart) in the same MRI scanner. At the time of both MRI examinations, all patients were physically and neurologically examined. The same specially trained and certified neurologist (JK) assessed Expanded Disability Status Scale (EDSS) scores in all patients (Kurtzke, 1983). From our cohort, 15 patients experienced one relapse between the two measurements, and one patient experienced two. Relapses occurred at least 1 month, and maximally 11 months after the first MRI measurement (mean = 4.5, SD = 3.5 months). Relative to the follow-up measurement, relapses occurred at least 1 month, and maximally 18 months before (mean = 10.8, SD = 5 months).
The study was approved by the interdisciplinary ethics committee of the University of Münster and the Physicians' Chamber of Westphalia-Lippe (Ärztekammer Westfalen-Lippe, 2010-378-b-S, 2017-754-f-S). Patients were informed of the study content in both oral and written form. All subjects gave written informed consent before participating in this study. (time interval between neuropsychological assessment and MRI was mean = 2.4, SD = 3.4 months). The SDMT is a short screening test for information processing speed in which numbers have to be assigned to abstract symbols as fast as possible. Visual attention and task switching abilities were assessed by the TMT, in which patients are instructed to connect a set of dots with numbers and letters according to specific rules. For all analyses including cognitive variables, we used percentile rank scores obtained by comparing patient's raw scores with testspecific normative data stratified for age and education. Average raw scores and SDs of all three assessments are depicted in Table 1.

| MRI
Patients were scanned in a 3T MAGNETOM Prisma fit MRI scanner using a 20-channel (matrix) head coil (Siemens Healthcare, Erlangen, Germany).
Employing the same MRI parameters and protocol for each patient at both time points, we obtained a 3D isotropic T1-weighted magnetization pre-

| Quantification of lesion volume
Lesion volume was estimated using the lesion growth algorithm (Schmidt et al., 2012) as implemented in the lesion segmentation toolbox (LST) version 2.0.15 (https://www.applied-statistics.de/lst.html) of the statistical parametric mapping (Penny & Henson, 2006) software (SPM12). T1-weighted volumes were segmented into cerebrospinal fluid, gray matter, and white matter maps. TSE-FLAIR images were coregistered to the T1-weighted volumes and information about the tissue classes used to compute lesion belief maps. These maps were thresholded with an initial threshold of κ = 0.1, which was chosen as the optimal value after visual inspection. A binary lesion map was obtained by growing lesions along voxels that appeared hyperintense in the TSE-FLAIR images.

| Data preprocessing, network reconstruction, and network parameter computation
The cortical surface was automatically reconstructed from T1-weighted images using the longitudinal processing stream (Reuter, Schmansky,  (Reuter & Fischl, 2011) using robust, inverse consistent registration (Reuter, Rosas, & Fischl, 2010). Reliability and statistical power is significantly increased by initializing with common information from this template, while performing processing steps including skull stripping, Talairach transforms, atlas registration, spherical surface maps, and parcellations (Reuter et al., 2012). DWI were preprocessed using FSL version 5.0.9 (https://fsl.fmrib.ox.ac.uk/fsl/), including motion artifact and eddy current correction (diffusion toolbox of FSL). Crossing fiber distributions and the probability of major and secondary fiber directions were computed. In order to track fibers through regions of crossing or complexity, a multi-fiber model was fit to the diffusion data at each voxel (Behrens et al., 2003). To estimate the probability distribution of connections from each voxel, 5,000 streamlines were drawn. Seed masks were individually defined based on the 68 cortical (Desikan et al., 2006) and 26 deep gray matter regions (Fischl et al., 2002) reconstructed with FreeSurfer. Tracts that passed through at least one other seed mask were retained.
The probabilistic tractography matrix was then generated for each subject. Here, each cell represented the number of streamlines passing through any two regions, normalized for the total number of streamlines drawn from a region. We computed graph theoretical measures of network modularity, distance, and local information transfer with the Brain Connectivity Toolbox (Rubinov & Sporns, 2010; https://sites.google.com/site/bctnet/), as these measures have been shown to be clinically relevant (Fleischer et al., 2017;Koirala et al., 2018). In the graph theoretical framework, brain regions are considered as nodes, and their connections as edges (Fornito, Zalesky, & Bullmore, 2016). Modularity is a measure of the degree, to which the network is subdivided into strongly interconnected nodes with few connections to other groups of nodes (Newman, 2006). The eccentricity of a node is the maximum of its finite distances to all other vertices (Sporns, 2003), and we computed the average of all eccentricity values within the network. Local efficiency, computed with Dijkstra's algorithm, is the average of the inverse distance matrix in subclusters and reflects to what extent neighboring regions of a node remain connected after this node is taken out of the network (Latora & Marchiori, 2001). An overview of the processing pipeline is depicted in Figure 1. FreeSurfer's Qdec (Hagler Jr, Saygin, & Sereno, 2006), and average thicknesses of vertices belonging to significant clusters were extracted for each subject.
For deep gray matter regions, we extracted average volume values as computed using FreeSurfer at both time points. We identified regions of interest based on the regional deep gray matter atrophy rate ((GM t2 − GM t1 )/t diff ), with GM t1 and GM t2 as the volume values at the time point of baseline (t 1 ) and follow-up (t 2 ) MRI examination, respectively, and t diff as the time interval between any two longitudinal measurements. Subsequently, we applied a false discovery rate (FDR)-correction to all 26 comparisons computed with t tests against zero atrophy using MATLAB (MathWorks, Natick, MA) and included significantly atrophied regions in our further analyses. To further resolve, which brain stem regions account for possible observed effects, we applied FreeSurfer's brain stem parcellation algorithm (Iglesias et al., 2015) and

| Analysis 3: Structural integrity and physical disability
We computed an ordinal regression to predict physical disability (EDSS score) at t 2 . In one model, we entered the network parameters modularity, eccentricity, and local efficiency at t 1 as variables of interest additionally to the control variables in a stepwise manner. In a second model, we entered deep gray matter volume and cortical thickness of atrophying regions at t 1 additionally to control variables again in a stepwise manner.

| Analysis 4: Structural integrity and cognitive performance
The relation between gray matter integrity, white matter network configuration, and cognitive performance was investigated using stepwise linear regression analysis. We inserted deep gray matter volume and cortical thickness at t 1 of regions of interest, as well as eccentricity, modularity and local efficiency at t 1 as variables of interest in three separate models each for SDMT, TMT-A, and TMT-B as dependent variable.
For all regression analyses, we included age, gender, lesion, and total brain volume as control variables. Variables with p < .01 were included in the final models.

| Analysis 1: Gray matter atrophy
We identified regions of interest by assessing cortical and deep gray matter atrophy over the 1-year measurement period. Left and right middle temporal (MT) cortices were significantly atrophied (both p < .05, Monte Carlo Simulation; Figure 2). Atrophied deep gray matter included bilateral thalamus, caudate, brain stem, and right putamen (p < .05, FDR corrected; Figure 3). After further parcellating the brain stem, we detected that the pons mainly accounted for the observed atrophy effects (p < .05, FDR corrected, while p > .05 for medulla, superior cerebellar peduncle, ventral diencephalon, and midbrain).

| Analysis 2: Gray matter integrity and white matter networks
We modeled white matter network reconfigurations regarding modularity, eccentricity, and local efficiency controlling for age, gender, lesion volume, and total brain volume at t 1 . Average volume of bilateral thalamus, caudate, brain stem, and right putamen as well as bilateral MT cortex thickness at t 1 as atrophying regions of interest were added in a stepwise linear regression analysis (Table 2). While lesion volume alone accounted for 23% of the variance of eccentricity rate (p < .001), 32% of the variance was explained by additionally including the volume of F I G U R E 1 Overview of the processing pipeline. Based on diffusion-weighted images, probabilistic tractography was conducted. This yielded a white matter network connectivity matrix, from which we computed network measures of modularity, distance, and local information transfer. Reconstructed from T1-weighted images, cortical thickness and deep gray matter volume of regions of interest was estimated. Finally, we set gray matter integrity and white matter network measures into relation with parameters of physical disability and cognitive impairment. Abbreviation: ROI, region of interest atrophying deep gray matter at t 1 (p = .002, Table 2). We performed the same analysis with modularity rate as dependent variable. The best model involved the volume of deep gray matter regions of interest at t 1 as predictor (p < .005; Table 2). Cortical thickness did not explain variance of eccentricity or modularity rate. Regarding local efficiency rate, neither deep gray matter volume, nor cortical thickness accounted for variance in addition to the covariates.
We next investigated the reverse direction, that is, whether atrophy can be predicted by network parameter values at t 1 . No additional variance could be explained by network parameters at t 1 either for deep gray matter, or for cortical atrophy (all p > .01).

| Analysis 3: Structural integrity and physical disability
Predicting physical disability, the strongest model of the ordinal regression analysis included deep gray matter volume at t 1 in addition to control variables, explaining 17% of the variance of EDSS score at t 2 (p = .003, Table 3). No other variable of interest improved the model significantly.
We identified models that predicted each eccentricity and modularity rate. In addition to the control variable lesion volume, volumes of deep gray matter regions of interest at baseline explained further variance in the change of eccentricity, where higher lesion volume and lower deep gray matter volume were associated with an increase in eccentricity over time. Lower deep gray matter volume in regions of interest was also predictive of an increase in modularity and of physical disability at 1-year follow-up. In a subsample of the patients, we related MRI variables to neuropsychological parameters. Lower lesion volume was associated with better performance in SDMT.
Female gender and higher deep gray matter volume in regions of interest were related to higher TMT-B ranks.

| Gray matter atrophy
Previous works have detected atrophy in deep gray matter regions of patients with MS, including the thalamus, caudate, putamen, and brain stem. These structures have a wide range of functions, including sensory and motor processing, which become restricted in MS due to the accumulating inflammatory and neurodegenerative processes (Deppe et al., 2016;Eshaghi et al., 2018;Steenwijk et al., 2016). We could show that the pons was the structure driving the atrophy effects in the brain stem. The importance of pathological processes in the pons has been outlined previously, as lesions in the pons are predictive of physical disability (Bakshi, Benedict, Bermel, & Jacobs, 2001). Cortically, we observed a significant reduction in thickness of MT areas.
This was related to cognitive impairment in a previous study, where the temporal lobe was found to be thinner among all cortical regions in cognitively impaired MS patients than in patients with preserved cognitive function (Tillema et al., 2016). However, here, we did not find an association with neuropsychological scores. This could be due to the small sample size of patients who were neuropsychologically examined or the type of neuropsychological assessments used.

| Gray matter integrity and white matter networks
Bringing gray matter integrity and white matter network characteristics together, we found that lower deep gray matter volume was predictive of an increase in the distance measure eccentricity. Higher eccentricity in a network indicates a disruption of information flow, which is rerouted via longer pathways. Deterioration of deep gray matter could cause such a rerouting of information, maintaining information flow, but with reduced efficiency. Equally, higher lesion load causes an increase in eccentricity. Measures of path distance have been shown previously to be increased in patients compared to healthy controls (Llufriu et al., 2017), probably a direct consequence of focal tissue damage. We further showed that increasing modularity can be predicted by lower deep gray matter volume. Previous work from our group revealed that patient networks initially decompose stronger into distinct modules compared to healthy brain networks (Fleischer et al., 2017). This finding raises the question how far variations in modularity and, more generally, functional reorganization of white matter networks are maladaptive or adaptive (Penner & Aktas, 2017;Rocca & Filippi, 2017). A modularity increase could result from a maladaptive breakdown in long-range connections, dividing the network in a higher number of smaller tightly interconnected units, or from an adaptive compensatory response with stronger intramodular connections (Fleischer et al., 2017). While it was possible to predict changes in white matter network properties based on gray matter atrophy, this was not the case for the reverse direction, such that atrophy was neither predictable in cortical, nor in deep gray matter structures. We hence assume that white matter networks adapt as a response to gray matter depletion.

| Structural integrity and physical disability progression
Lower volume of atrophying deep gray matter regions at baseline served as the best predictor for physical disability at 1-year follow-up.
We, therefore, found that lower deep gray matter volume not only causes adaptations in white matter networks, it is also related to physical disability progression. Conversely, network parameters did not explain any additional variance. This again supports the view that white matter networks initially are less affected by the disease than gray matter integrity. For investigations of patients with longer disease duration, we would expect a relation between white matter network parameters and physical disability as well.

| Structural integrity and cognitive performance
We additionally investigated the relation of gray matter integrity and white matter network properties with cognitive performance in a part of the patient sample at baseline. Lower lesion volume was related to an increase in information processing speed, as quantified by SDMT.
Male patients scored lower in the TMT-B, reflecting executive functioning. Gender of healthy subjects does not correlate with TMT scores (Tombaugh, 2004). However, since males are less frequently, but more severely affected by the disease (Greer & McCombe, 2011), an adverse effect of gender on both physical disability and cognitive scores was not unexpected. To further explore this, we additionally compared median EDSS scores between males and females within the subgroup, and observed higher disability in male patients. Our results are in line with prior studies that assessed cognition in MS patients and observed stronger decline in males (Lin et al., 2016;Schoonheim et al., 2012), which moreover was related to EDSS scores (Savettieri et al., 2004). In our cohort, better performance in the TMT-B was further associated to higher deep gray matter volume in our regions of interest, again in accordance with previous research (Damjanovic et al., 2017;Modica et al., 2016;Schoonheim et al., 2012). The less complex TMT-A is potentially too insensitive to reflect initial variations in gray matter integrity or network configurations but might be more suitable to distinguish patients of different disease stages.
For cognitive performance, we would also expect a relation to white matter network measures after longer disease progression. However, in this retrospective study, only scores from the first measurement time point were available; hence, this needs to be investigated in future work.

| Limitations and prospect
Our patient cohort was rather homogeneously mild impaired regarding physical disability and cognitive impairment.  (Fleischer, Koirala, et al., 2019). Such longitudinal studies including sensitive markers of pathology are rare, however, are increasingly needed for early predictions of the clinical outcome. A complicating factor in MS research is the heterogeneity of patient samples. The disease and its manifestation in the central nervous system is very heterogeneous, and additional variables such as age, gender, treatment, relapse rate, and numerous other parameters are challenging to account for in analyses. In our study, we only considered early stage MS patients; however, it would be interesting to investigate other disease types. For example, more severe cognitive impairment has been observed in MS patients with the primary progressive type (Johnen et al., 2017), while only a few studies have analyzed underlying network reorganization processes . Considering the variety of confounding variables, analyses would be most informative if they were based on large patient cohorts, for example, from multicenter studies, allowing several sub-analyses by controlling different influences.

| Conclusion
Our results suggest (I) that neurodegeneration predominantly appears in deep gray matter regions after disease onset, and (II) that white matter network reconfiguration occurs as a response to deep gray matter pathology in patients with CIS and early RRMS. Physical disability progression and cognitive deterioration appear to depend on deep gray matter atrophy and the capacity of the network to adaptively reconfigure and maintain information processing efficiency.

ACKNOWLEDGMENTS
We thank all patients who participated in the study and Rosalind Gilchrist for proofreading the manuscript. Parts of this research were conducted using the supercomputer Mogon and advisory services offered by Johannes Gutenberg University Mainz (hpc.uni-mainz.de), which is a member of the AHRP and the Gauss Alliance e.V. We gratefully acknowledge the computing time granted on Mogon.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.