Global Rearrangement of Degree Centrality Reflects Cognitive Impairment and Fatigue in Multiple Sclerosis

2.1 Study Design and Participant Selection

The secondary analysis was performed on cross-sectional imaging and behavioral data of 65 PwMS and 65 HC matched for age and sex, with participants recruited from multiple sclerosis centers across Czechia. The same cohort has been partially analyzed and published using other methods.^6,26 Original inclusion criteria were: positive diagnosis of multiple sclerosis;²⁷ spastic paraparesis as a prominent clinical feature; stable clinical status for at least three months preceding the study (determined by neurologist); physical ability to undergo clinical testing – consistent with the EDSS score ≤ 7.5. Subjects were excluded in case of missing imaging data or imaging artifacts and conditionally excluded if exceeding the motion outlier criteria (see Section 2.5 Data pre-processing and analysis below). To address potential attrition bias, a sensitivity analysis was performed in a sample including motion outlier subjects.

2.2 Clinical Assessment and Questionnaires

Clinical parameters (regressands) comprised an objective assessment of cognitive performance using a Symbol Digit Modalities Test (SDMT)²⁸ and possible confounding factors such as fatigue (Fatigue Scale for Motor and Cognitive Functions, FSMC),²⁹ global disability (Expanded Disability Status Scales, EDSS),¹² and motor performance (Timed Up and Go Test, TUG).³⁰

2.3 Standard Protocol Approvals, Registrations, and Patient Consents

The secondary data analysis was pre-registered at osf.io (https://osf.io/v8ejw). The original study was approved by ethics committee of the Ethics Committee of the Faculty Hospital Královské Vinohrady, approval No. EK-VP/22/0/2014. All patients gave their written informed consent to participate in the study.

2.4 MRI Data Acquisition

Imaging was performed using a 3T magnetic resonance scanner (Siemens Trio Tim, Erlangen, Germany) equipped with a 12-channel phased-array head coil. The MRI protocol included blood oxygenation level-dependent (BOLD) resting-state fMRI (repetition time [TR] = 2,500 ms, echo time [TE] = 30 ms, flip angle = 70°, voxel size = 3×3×3 mm, 44 contiguous slices, 240 volumes, acquisition time = 10 min), as well as high-resolution 1-mm T1-weighted and T2-weighted imaging, and 2-mm diffusion-weighted imaging (DWI). Detailed acquisition parameters have been published elsewhere.^6,26

2.5 Data Pre-processing and Analysis

The initial pre-processing consisted of motion correction, correction of susceptibility-induced distortions and normalization to standard space using Advanced Normalization Tools (ANTs, v2.3.5.dev212-g44225).³¹ Brain coverage was evaluated using Mask_explorer.³² Next, motion outlier detection, anatomical component-based denoising procedure, and band-pass filtering within the frequency range 0.008 Hz – 0.09 Hz were applied in CONN toolbox v. 21a.³³ Participants with maximum volume-to-volume displacement exceeding 2 mm and/or mean volume-to-volume displacement exceeding 2 standard deviations above the sample mean (i.e., >0.26 mm) were marked as outliers (n = 13). In the main analysis, all remaining steps were performed after excluding the outliers. However, additional sensitivity analysis consisting of the same steps was conducted in a parallel pipeline in the sample with outliers. Please also refer to Supplementary Methods for further details on pre-processing.

Following the pre-processing, subject-specific functional connectivity matrices containing Fisher-transformed Pearson’s r coefficients were computed in CONN for 4,632 large voxels created using 6-mm resampling of a common gray matter mask, see Supplementary Methods for more details. Degree centrality was then computed for each subject using a brain connectivity toolbox (BCT, https://sites.google.com/a/brain-connectivity-toolbox.net/bct/) with 10% link density.²⁴ Regional degree centrality was extracted by averaging nodal degree from the DMN (4 ROIs),^18,19,21 basal ganglia (6 ROIs),^19,21 thalamus, hippocampus, and cerebellum (5 ROIs),²¹ and from the multimodal ROIs explicitly participating in SDMT: superior parietal lobule (2 ROIs), dorsolateral prefrontal cortex (2 ROIs), and anterior cingulate cortex (ACC; 1 ROI),¹⁵ see Supplementary Fig. S1 and Supplementary Table 1 for detailed ROI description. Finally, k_D was calculated using custom Matlab scripts implementing a modified approach according to Mansour et al.²⁴ In brief, k_D was obtained for each individual using a linear regression (y = k_D *x + b), where y = individual nodal degree – mean control group nodal degree, x = mean control group nodal degree, and b = intercept of the regression, see Supplementary Methods for more details.

The pre-processing of T1-weighted and DWI data, calculation of gray matter volume (GMV) as a measure of cortical atrophy and extraction of fractional anisotropy (FA) as a measure of white matter integrity, as well as ROI definition for GMV and FA is described elsewhere.⁶ Finally, the lesion load (LL) was calculated using the Lesion Segmentation Tool (LST, www.statistical-modelling.de/lst.html with lesion prediction algorithm.³⁴

2.6 Statistical Analysis

Initially, normality was assessed for all continuous variables using Kolmogorov-Smirnov test. Non-parametric tests were applied in case normality was violated. Additionally, regressors considerably deviating from normal distribution (LL) were log-transformed prior to any subsequent analysis to meet regression model assumptions. Demographic variables were compared between groups using Fisher’s exact test and Mann-Whitney U test. Pairwise deletion was applied in case of missing clinical data. A summary of all variables and outcome measures is provided in Supplementary Table 2.

The hypotheses were evaluated using following tests: (1) Group differences in k_D(primary outcome measure) between PwMS and HCs were assessed using a Mann-Whitney U test. (2) Differentiation between PwMS and HC based on k_D was assessed using (2a) area under the receiver operating characteristic (ROC) curve (AUC) and compared with differentiation based on regional degree centrality in each ROI using pair-wise asymptotic significance. For comparison with structural imaging parameters (2b), multiple logistic regression χ² statistics was used to test whether model with k_D, GMV, FA, and log(LL) significantly improved differentiation between PwMS and HCs compared to a nested model without k_D as a regressor. The correlations between (3) k_D and SDMT, (4a) regional degree centrality and SDMT, as well as auxiliary correlations between k_Dand (5) additional clinical variables (EDSS, FSMC, TUG) or (7) structural imaging parameters (GMV, FA, LL) were evaluated using Spearman’s rank correlation coefficient. For correlations of regional degree centrality (4a), Bonferroni-Holm correction for multiple comparisons across 18 ROIs was applied (α = 0.0028). To test whether k_D improves regression of clinical scores on top of structural imaging, ordinal regression models including k_D, GMV, FA, log(LL), age, sex, and disease duration as regressors of (4b) SDMT, (6a) EDSS, (6b) FSMC, and (6c) TUG were compared to nested models without k_D by means of likelihood ratio test. (8) Finally, differences in all imaging parameters (regional degree centrality, global GMV, global FA, and lesion load) were assessed using multiple Mann-Whitney U tests with Bonferroni-Holm correction for multiple comparisons. All tests were performed using SPSS v29.0.1.1 (IBM, Armonk, NY, USA). The p < 0.05 was considered significant. One-tailed tests were used where superiority was assumed by the hypotheses or statistics with one-tailed distribution was employed (2a, 2b, 4b, and 6), with two-tailed tests applied otherwise. Details on figure preparation are provided in Supplementary Methods.

2.7 Post-hoc Analyses

In order to visualize local contributions to significant global correlations between k_D and clinical scores, a post-hoc voxel-wise analysis was performed in randomise, part of FSL v. 6.0.3.³⁵ First, nodal degree centrality was back-projected to the original 6-mm voxels. Next, general linear model with two-sample t-test (group differences) and regression contrast was employed to evaluate the correlation using non-parametric threshold-free cluster enhancement (TFCE) correction for multiple comparisons with 10,000 permutations and family-wise error corrected alpha = 0.05.

In case of significant correlations with EDSS, FSMC or TUG, correlation with regional degree in the 18 pre-defined ROIs was additionally assessed using Spearman rank correlation to allow further interpretation.

2.8 Power Analysis

No k_Ddata in PwMS for a power analysis were available. Based on a published difference in k_Dbetween patients and HC of |k_D| = 0.21 in a different patient cohort,²⁴ analysis results in a minimum sample of 12 participants to achieve power to detect a significant correlation of 90%. With the existing data set (n = 64) differences of down to |k_D| = 0.09 can be identified with the same power of 90%.

3.1 Study Sample

Out of the original sample of 65 PwMS and 65 HC, one PwMS and one HC were excluded due to missing data (incomplete field of view and susceptibility artifact), another PwMS was excluded due to a suspected vascular lesion. In the remaining sample, 13 subjects with excessive motion levels were identified (7 PwMS and 6 HCs), see Supplementary Fig. S2 for inclusion/exclusion diagram. Here, only results in 56 PwMS and 58 HCs after excluding motion outliers are reported (“final” sample), whereas results in the sample with outliers are provided in Supplementary Results. While both analyses yielded mostly similar results, two differences are explicitly stated below. Demographic details of the “final” sample and summary statistics for clinical parameters are provided in Table 1.

3.2 Group Differences and Differentiation between PwMS and HCs

PwMS showed significantly lower k_D compared to HCs (PwMS: median = −0.298, inter-quartile range [IQR] = 0.549; HCs: median = −0.058, IQR = 0.542; p = 0.009, Mann-Whitney U test), see Fig. 1. Underlying raw degree centrality data are summarized in Supplementary Fig. S3.

• Download figure
• Open in new tab

Fig. 1. Group differences in degree rank order disruption index (k_D).

Raincloud plots illustrating individual k_D values and distribution. Patients with multiple sclerosis are shown in green, healthy controls in red. Mann-Whitney U test p value is provided in annotation.

The ROC analysis for differentiation between PwMS and HCs yielded significant above-chance area under curve (AUC) for k_D (AUC = 0.642, p = 0.007, two-tailed asymptotic significance for null hypothesis AUC = 0.5), the left lateral parietal portion of the DMN (DMN-LLP; AUC = 0.671, p = 0.001), and the ACC (AUC = 0.619; p = 0.026), see Fig. 2 and Supplementary Table 3. In pair-wise comparisons, AUC for k_D was significantly higher than AUC for 6 ROIs and did not significantly differ from the remaining ROIs (Supplementary Table 3).

• Download figure
• Open in new tab

Fig. 2. Receiver operating characteristic (ROC) analysis.

ROC curves for differentiation between patients with multiple sclerosis (MS) and healthy controls using degree rank order disruption index (k_D; area under curve [AUC] = 0.642; p = 0.007, two-tailed asymptotic significance for null hypothesis AUC = 0.5, uncorrected); the left lateral parietal portion of the default mode network (DMN-LLP; AUC = 0.671; p = 0.001), and the anterior cingulate cortex (ACC; AUC = 0.619; p = 0.026).

We observed no significant improvement in a multiple logistic regression model differentiating between PwMS and HCs after adding k_D as an additional regressor on top of GMV, FA, log(LL) (χ² step = 0.007, p = 0.934).

3.3 Correlation with Cognitive Processing Speed

We detected a significant correlation between k_D and SDMT (Spearman’s ρ = 0.32, p = 0.019, n = 55, Fig. 3). For the regional degree centrality, significant correlation with SDMT was observed in the medial prefrontal part of the DMN, yielding slightly higher effect size than the k_D, see Table 2. An ordinal regression model including GMV, FA, log(LL), age, sex, and years since diagnosis as regressors of SDMT score was significantly improved after adding k_D (χ² step = 4.49, p = 0.034, likelihood ratio test, see Table 3). In contrast, neither significant correlation with SDMT nor improvement of regression model for SDMT were observed in analysis with motion outliers (see Supplementary Results).

• Download figure
• Open in new tab

Fig. 3. Correlation between k_D and clinical scores.

Scatter plots illustrating relationship between the degree rank order disruption index (k_D) and cognitive processing speed (Symbol Digit Modalities Test, SDMT), global disability (Expanded Disability Status Scales, EDSS), fatigue (Fatigue Scale for Motor and Cognitive Functions, FSMC), and motor performance (Timed Up and Go Test, TUG). Spearman’s rank correlation coefficient ( ρ), two-tailed uncorrected significance, and number of valid observations are provided.

3.4 Correlation with Fatigue, Global Disability, and Motor Performance

We detected a significant correlation between k_D and FSMC (Spearman’s ρ = −0.34, p = 0.010, n = 56), but not for EDSS (Spearman’s ρ = −0.06, p = 0.674, n = 56) or TUG (Spearman’s ρ = −0.18, p = 0.211, n = 52), see Fig 4. In ordinal regression, k_D significantly improved the model fit for fatigue (FSMC) when added on top of GMV, FA, log(LL), age, sex, and years since diagnosis, but not for EDSS or TUG (Table 3).

• Download figure
• Open in new tab

Fig. 4. Voxel-wise correlation between degree and clinical scores.

Color overlays illustrating spatial distribution of correlation of the voxel-wise degree centrality with (top) cognitive processing speed (Symbol Digit Modalities Test, SDMT) and (bottom) fatigue (Fatigue Scale for Motor and Cognitive Functions, FSMC). Red-yellow and blue-lightblue overlays indicate positive and negative correlation, respectively. Thresholded using non-parametric threshold-free cluster enhancement (10,000 permutations, family-wise error-corrected p = 0.05). In the top panel, yellow overlay shows overlap with the medial prefrontal cortex default mode network (MPFC-DMN) region of interest.

3.5 Relationship between k_D and Structural Imaging Biomarkers

We observed no significant correlation between k_D and structural imaging parameters, i.e., GMV, LL, and global FA, see Supplementary Fig. S4. In analysis with motion outliers, however, k_D was significantly correlated with LL (see Supplementary Results). All structural imaging parameters significantly differed between PwMS and HCs, see Supplementary Table 4.

3.6 Post-hoc Analyses

Voxel-wise group differences are illustrated in supplementary Fig. S3. Voxel-wise regression analysis for degree centrality as response variable and SDMT as explanatory variable yielded significant cluster with negative effect mainly in MPFC and ACC, and to a lesser degree in the subcallosal cortex and the right nucleus accumbens, overlapping in part with the DMN-MPFC ROI (Fig. 4 and Fig. S3 for unthresholded data). The regression for FSMC yielded a positive effect in the predominantly right cerebellum, right caudate nucleus, right inferior frontal gyrus and a negative effect in the left temporo-occipital fusiform and lateral occipital cortex (Fig. 4 and Fig. S3 for unthresholded data).

The ROI analysis for FSMC yielded significant correlation in the left and right caudate nuclei (left: Spearman’s ρ = 0.41, p = 0.002; right: Spearman’s ρ = 0.40, p = 0.002; n = 56, Bonferroni-Holm-corrected across the 18 ROIs, i.e., alpha = 0.0028). See Table 2 for complete results.

4.1 Differentiation between PwMS and HCs

Our primary observation of decreased k_D in PwMS captures the global character of changes in nodal centrality (both degree and eigenvector) that have recently been reported on the local and network-wide level in MS.^{9,18,20,21,36} Lower k_D in PwMS suggests less centralized and more diffusely distributed rsFC, which is analogous to a previously described disruption of the rich-club topology of the brain network in MS.³⁷ However, our multiple regression analysis indicated that k_Ddid not improve differentiation between PwMS and HCs when added on top of structural imaging parameters. Taken together with the considerable overlap between k_Ddistributions in PwMS and HC (Fig. 1), it can be inferred that degree reordering is predominantly driven by factors other than presence of MS. In fact, changes in centrality of the network hubs have been consistently observed to depend on the cognitive status of PwMS.^19–22 As such, we conclude that changes in k_D are more likely related to the resulting neurological deficits than to the underlying MS-related pathology per se.

4.2 Correlation with Cognitive Processing Speed

Our next main analysis demonstrated that k_D was a significant explanatory variable for cognitive processing speed when compared to global structural imaging parameters. While the negative correlation between k_D and SDMT is a novel finding, it is in line with previous evidence for weak positive association between SDMT score and rsFC of the peripheral nodes outside the rich club.³⁷ The correlation between k_D and SDMT yielded an effect size similar to some previously reported structural imaging biomarkers of cognitive processing speed, including fractional anisotropy in the superior longitudinal fascicle¹⁶ or gray matter atrophy-based brain age gap,³⁸ but was lower than overall effect size in a recent meta-analysis of multimodal structural MRI data.³⁹

In our dataset, a superior correlation with SDMT was achieved using regional degree centrality in the MPFC hub of the DMN, yielding an effect size comparable with structural imaging biomarkers.^39,40 Correspondingly, cognitive impairment in MS has been shown to be associated with increased degree or eigenvector centrality in the DMN^19–23 and decreased centrality in the visual^20,22,23 and sensorimotor network,²⁰ while eigenvector and degree centrality show high agreement even within the same group.²¹ Nevertheless, our results add novel evidence for the strength of the relationship between degree centrality in the DMN and cognitive processing speed since no such correlation in MS has been reported before.

4.3 Role of the DMN in Pathophysiology of Cognitive Deterioration in MS

The ROI and post-hoc voxel-wise analyses indicated that the MPFC and ACC (i.e., anterior DMN) provided the highest correlation between degree and cognitive processing speed. The role of DMN in CI remains controversial.¹ Although rsFC abnormalities in DMN have been considered to be non-specific, i.e., linked to the dysfunction of the entire brain network,¹ an excessively central and less dynamic DMN has also been proposed to directly hinder externally oriented cognitive processing by superfluous introspective thoughts.^20,41 In contrast, cognitive rehabilitation has been associated with increased activation and rsFC of the DMN, suggesting its possible participation in compensatory mechanisms.⁴² Additionally, some rsFC studies using the Paced Auditory Serial Addition Test (PASAT), i.e., another frequently employed cognitive test in MS research, also point to dysfunction of DMN and subcallosal cortex.^36,40 Hence, graph-theoretical measures extracted from anterior portion of the DMN are potential future candidates for even more accurate biomarkers of cognitive processing speed than k_D, although studies using more elaborate cognitive testing are warranted to assess their specificity with respect to other cognitive domains.

4.4 Association with Fatigue

Our exploratory analyses indicated, that apart from cognitive processing speed, k_Dwas a significant explanatory variable for the self-reported global fatigue score (FSMC). While cognitive performance and fatigue have been shown to be associated with similar rsFC dysfunctions, such as increased rsFC in posterior DMN and reduced rsFC in the anterior DMN,⁴³ our post-hoc analyses on ROI and voxel-wise level suggested potential differentiation between CI and global fatigue: Whereas SDMT correlated with average degree in DMN-MPFC, fatigue scores were more strongly associated with degree in the caudate nuclei, cerebellum, and fusiform cortex (Table 2 and Fig. 4).

Although cytokine and neuroendocrine abnormalities have been discussed as potential causes of fatigue in MS,^44,45 it has also been associated with damage to cortico-subcortical pathways and with particular involvement of prefrontal cortex.⁴⁵ On the computational (metacognitive) level, it has been proposed to result from mismatch between predicted and measured output from cognitive and sensorimotor networks.^1,44 Our results fit in by emphasizing the role of basal ganglia⁴⁶ and cerebellum, which is involved in maintaining internal forward models and error monitoring.⁴⁷ Future dedicated studies should evaluate the specificity and stability of the here identified biomarkers and their accuracy with respect to the motor and cognitive sub-domains of FSMC.

4.5 Relationship with Global Disability and Motor Performance

Despite some evidence for increased eigenvector centrality in the basal ganglia,¹⁹ DMN,^18,19 and lower degree in executive control and visual networks³⁶ in patients with higher EDSS, we found no association between k_D and global disability (EDSS) and gait (TUG). In contrast, previous machine learning study in the same dataset achieved a medium correlation (not exceeding r = 0.5) with EDSS or TUG using rsFC, white matter integrity and gray matter atrophy features either separately or combined.⁶ We therefore speculate that motor deficits are better explained by local changes in rsFC than by overall network disarray. This is supported by fact that PwMS in our paraparetic cohort had significantly lower degree in supplementary motor area and paracentral lobule than HCs as shown by our post-hoc analysis (Supplementary Fig. S4). This may be, at least in part, due to strong influence of strategic lesions in pyramidal or cerebellar pathways (including the spinal cord) on motor performance including gait,^48,49 which dominates the EDSS score.¹²

4.6 Correlation with Structural Imaging Parameters

Our results showed that k_D was largely independent on global structural imaging parameters, such as brain atrophy, global white matter integrity or lesion load, supporting its applicability as an add-on biomarker for some of the clinical scores. In previous studies, however, degree centrality in the thalamus and cerebellum correlated with the lesion load,¹⁸ while degree centrality in the insular cortex, cerebellum and precuneus (DMN) correlated with brain and gray matter volume.¹⁸ Correlation between LL and degree centrality of the executive and salience networks was also described.³⁶ We speculate that the lack of correlation between structural imaging and k_D might reflect its closer relationship with actual functional reserve that may vary even in patients with similar degree of structural damage.

4.7 Limitations and Future Directions

While the main strength of the study was rigorous pre-registration of all main analyses, there are also several limitations related to the fact that this study was carried out as a secondary analysis: SDMT reflects mainly cognitive processing speed and involves visual processing. Further assessments across multiple cognitive domains and sensory modalities, as well as consideration of depression as a potential confounding factor are warranted. Furthermore, our cross-sectional data cannot differentiate between state (e.g., affected by the current level of exhaustion) and trait (i.e., effect of structural damage and subsequent (mal-)adaptive plasticity) as well as premorbid patient characteristics. Normative data, assessment of test-retest reliability and longitudinal evaluation are thus needed. Future directions also include assessment of time-resolved k_D, as temporal dynamics of centrality seem to provide additional insights into network dysfunction underlying CI.²⁰ In case of fatigue, more fine-grained distinction into motor and cognitive fatigue would also further aid the interpretation. Finally, our data cannot be currently generalized to all individuals with MS (see our inclusion criteria), hence, a cross-validation of our results in an independent dataset is necessary before translating the results into clinical practice.