Functional connectivity dynamics reflect disability and multi-domain clinical impairment in patients with relapsing-remitting multiple sclerosis

Highlights • Patients with early-MS show distinct alterations in functional connectivity dynamics.• Multi-domain clinical outcome scores scale with disability status scores in patients.• Temporal network dynamics vary with disability status in early MS.• Frontoparietal–basal ganglia decoupling is a novel, clinically-relevant signature of MS.


Introduction
Identifying aberrant signatures of brain activity in multiple sclerosis (MS) is a key research target for improving our understanding of the complex interplay between symptom severity and functional network behavior, especially in early MS. However, the characterization of common and clinically relevant functional alterations has been challenging, evidenced by a growing body of literature reporting complex, and occasionally, contradictory patterns of functional brain changes in MS.
Recent studies highlight subcortical brain regions as key players in the functional changes seen in MS: altered functional connectivity (FC) of the thalamus (Lin et al., 2019;Hidalgo de la Cruz et al., 2018;Tona et al., 2014) and basal ganglia (BG) (Finke et al., 2015;Tijhuis et al., 2021) has been consistently linked to cognitive function and fatigue. In contrast, while many studies implicate default mode network (DMN) dysfunction in MS, the directionality and thus the interpretation of this dysfunction remains inconsistent: both increased (Hawellek et al., 2011;Meijer et al., 2017) and opposingly, decreased DMN connectivity  have been identified in the relationship between functional alterations and cognitive impairment in MS.
Potential explanations for such heterogeneity include the temporally dynamic evolution of the disease, the ensuing variance in disease severity and disability, its diverse cognitive and behavioral phenotypes as well as methodological variability across functional MRI (fMRI) studies. Furthermore, much previous research has been limited to a static account of brain activity in which FC is computed on signals over the entire scan duration. However, recent advances in time-resolved approacheswhere FC is computed on a finer temporal scale and clustered into recurrent connectivity stateshave made substantial contributions to identifying functional brain signatures in various states of health  and disease Demirtaş et al., 2016;von Schwanenflug et al., 2022).
Furthermore, evidence from recent investigations into the hierarchical organization of brain dynamics in healthy populations has revealed a principled temporal gradient in which networks operate (Vidaurre et al., 2017). This gradient from primary sensorimotor systems to higher-order cognitive systems is thought to reflect the variance in timescales of information encoding speed and activity across networks (Raut et al., 2020). These differences further highlight the need for continued research into brain dynamics using time-varying techniques, which more sensitively assess moment-to-moment changes that are not detectable with static approaches.
Such time-resolved techniques have recently highlighted reduced network centrality (Eijlers et al., 2019) andFC dynamics (d'Ambrosio et al., 2020) in cognitively impaired MS patients, as well as subtypespecific altered FC and associations with cognitive and motor impairment (Hidalgo de la Cruz et al., 2021). While dynamic whole-brain connectivity studies thus represent a promising avenue to unravel the relationship between functional alterations and clinical outcomes, it remains unclear how altered network dynamics relate to disability status, disease severity, and the multi-domain impairments commonly observed in MS.
Here, we characterize the spatiotemporal patterns of FC alterations in a large sample of patients with early relapsing-remitting MS (RRMS). To this end, we investigate both static and dynamic FC alterations compared to healthy participants and relate these functional dynamics to disability status and clinical impairment across multiple domains including cognitive, behavioral, and affective scores as well as structural brain measures.

Study participants and clinical assessment of patients
Patients were recruited from the outpatient clinic of the Clinical and Experimental Multiple Sclerosis Research Center at Charité -Universitätsmedizin Berlin from 2014 to 2018. At the time of their study visit, all patients fulfilled the McDonald criteria (Thompson et al., 2018) for an RRMS diagnosis (Table 1). Patients did not have any comorbid neurological or psychiatric conditions. Healthy control participants (HC) were recruited using the Charité intranet and were without neurological or psychiatric conditions.
Patients underwent a complete neurological assessment and were evaluated by a board-certified neurologist. Clinical and neuropsychological evaluation occurred on the same day or within one day of MRI scan acquisition. A patient's current level of disability was evaluated using the Expanded Disability Status Scale (EDSS). Motor performance speed was assessed using the Timed 25-Foot Walk Test and the Nine-Hole Peg Test. Visual acuity was assessed using both the Early Treatment Diabetic Retinopathy Study (ETDRS) and Sloan scales at high, medium, and low contrasts.

Neuropsychological assessment
Patients also underwent an extensive battery of neuropsychological tests. Fatigue was assessed using the Fatigue Severity Scale (FSS). Depressive symptoms were assessed using the Beck's Depression Inventory II (BDI-II). It is important to note that, while there was considerable variance in depressive symptoms, average BDI-II raw scores did not suggest manifest clinical depression.
Cognitive performance was assessed using the Brief Repeatable Battery of Neuropsychological Tests (BRB-N), including the following sub-tests: the Selective Reminding Test (SRT) as a measure of verbal learning and memory, the 10/36 Spatial Recall Test (SPAT) as a measure of visuospatial memory, the Symbol Digit Modalities Test (SDMT) as a measure of processing speed and attention, the three-second version of the Paced Auditory Serial Addition Test (PASAT) as a measure of information processing speed and calculation ability, and the Word List Generation Test (WLG) as a measure of verbal fluency.

Group comparisons by disability scores and clinical rationale
Age-and sex-matching of patients and HCs was performed using a custom matching algorithm. Thus, 101 patients and 101 matched HCs were initially included; 92 patients had an EDSS rating from the date of scan acquisition. For the purposes of group comparisons, patients were split into subgroups of mild to moderate disability (EDSS ≥ 2, n = 39) and no disability (EDSS ≤ 1, n = 36) using the upper and lower 30th percentile threshold. Seventeen patients with an EDSS score of 1.5 were thus excluded because their score corresponded exactly to the median of the distribution across patients. Patients with EDSS ≥ 2 were significantly older and had a longer disease duration than patients with EDSS ≤ 1 (Table 1). The final sample used for statistical analysis of group differences thus included 75 patients and 75 HCs. Additional details on the matching procedure and the patient cohort split are provided in supplemental material (SM) sections 1.1 and 1.2.
Patients in the "no disability" subgroup are those with an EDSS rating of 0 or 1. These are patients with either a normal neurological exam (EDSS = 0) or those with no disability but minimal signs of impairment in only one of the functional systems evaluated within the EDSS assessment (EDSS = 1). Overall, these patients are characterized by the absence of clinical disability (Kurtzke, 1983). Patients in the "mild to moderate disability" subgroup are those with an EDSS rating of 2 or greater. Although the range of EDSS ratings in this group was 2 through 5.5, the median EDSS was 2.5 (IQR = 1); no patient had a score of 4.5 or 5, and only one patient had a score of 5.5 ( Supplementary Fig. 2). Thus, this group of patients showed minimal to moderate levels of disability in one or more functional systems incorporated in the EDSS assessment (Kurtzke, 1983).

MRI acquisition
MRI data were acquired at the Berlin Center for Advanced Neuroimaging at Charité-Universitätsmedizin Berlin, Germany, using a 20channel head coil on a 3 T Trim Trio scanner (Siemens, Erlangen, Germany). For each participant, the following sequences were acquired and used in the present study: a 10-minute resting-state-fMRI (rs-fMRI) scan was collected using a repetition time (TR) of 2250 ms (TE = 30 ms, 260 volumes, matrix size = 64 × 64, 37 axial slices, slice thickness = 3.4 mm, voxel size = 3.4 × 3.4 × 3.4 mm 3 ); a T1-weighted structural scan was acquired using a magnetization-prepared rapid gradient echo (MPRAGE) sequence (1 mm 3 isotopic resolution, matrix size = 240 × 240, 176 slices); and a T2-weighted fluid-attenuated inversion recovery (FLAIR) sequence was acquired with a TR of 5000 ms (TE = 388 ms, FOV = 250 × 250 mm 2 , matrix size 250 × 250, 76 slices, slice thickness = 1 mm, 1 mm isotopic resolution). All participants gave written informed consent, and the study was approved by the local ethics committee.

Conventional MRI analysis
T2-hyperintense brain lesion segmentation was performed manually using ITK-SNAP by two MRI technicians with more than 10 years of experience in MS research. The Functional Magnetic Resonance Imaging of the Brain (FMRIB) Software Library (FSL) "cluster" and "fslstats" tools were used to calculate lesion volumes and counts. Lesion load was taken as a patient's total lesion volume in milliliters. Total normalized brain volume was calculated from lesion-filled MPRAGE scans using FSL SIENAX.

Calculation of multi-domain clinical outcome z-scores
To systematically assess the relationship between functional brain measures and variance across the patient cohort in neuropsychological, behavioral, and structural MRI measures, we computed a custom set of domain-specific clinical outcome scores. Individual test scores and structural MRI measures were grouped into the following seven domains: cognition, motor, vision, depression, fatigue, lesion load and brain volume. Table 2 shows a detailed list of measures assigned to each domain. For each domain, a z-score was calculated according to the following procedure: the summary or total score for each test was calculated across its items according to the respective test manual; the test's z-score was then computed across all patient scores; for domains that included multiple tests, a composite z-score was computed by averaging across all tests assigned to that domain. Thus, each patient ultimately obtained a single value for each domain. Please note, to account for directionality differences, test scores for vision, cognition, and brain volume were re-coded by multiplying by − 1. Thus, all domain zscores represent impairment indices, where higher values signify worse performance, more severe symptoms, or higher brain atrophy, respectively.

Preprocessing
All analyses were performed with MATLAB (R2019b) and R (3.6.3). Preprocessing was performed using the CONN toolbox (version 18.b;(Whitfield-Gabrieli and Nieto-Castanon, 2012). The following preprocessing steps were applied using each participants' rs-fMRI and structural scans, in line with the default settings of the CONN toolbox: discarding of first five volumes of the rs-fMRl scan to account for scanner gradient stabilization, realignment to the first volume using b-spline interpolation, slice-timing correction using sinc-interpolation to timeshift and resample functional data slices, direct tissue segmentation and spatial normalization to standard Montreal Neurological Institute space including resampling of functional data to 2 mm isotropic voxels and structural data to 1 mm isotropic voxels, and finally, spatial smoothing by convolving with a Gaussian kernel of 8 mm full width at half maximum.

Group independent component analysis (GICA)
The "spatial GICA" feature of the Group ICA of fMRI toolbox (GIFT v4.0b, http://mialab.mrn.org/software/gift/index.html) was then used to reduce all participants' preprocessed rs-fMRI data into 100 group spatial independent components. To this end, each participant's 4Dfunctional data (255 volumes) were reduced to 150 temporal principal components (PC) that were maximally independent. All participant components were then concatenated, and a second principal component  (Calhoun et al., 2001) for ICA was used to extract 100 group independent components (ICs), which were then checked for stability in ICASSO (Himberg and Icasso, 2003) by performing 20 iterations of the ICA algorithm. Lastly, for each participant, a set of 100 component spatial maps and corresponding component time-courses were extracted with the spatiotemporal regression back-reconstruction algorithm (Du and Fan, 2013) in GIFT. This method has been extensively described in Calhoun et al., 2001).

Component rating and network assignment
Group ICs were manually classified as signal or noise based on criteria recommended in (Griffanti et al., 2017). This was performed by two independent raters (AR, NS), and any disagreement was resolved through rating by a third independent rater (CF). Forty-seven group ICs were classified as "signal", while 63 were discarded as "noise". Signal components were assigned to the following cortical resting-state networks (RSN) according to (Yeo et al., 2011): ventral attention (vATT), dorsal attention (dATT), somatomotor (SMN), DMN, visual (VIS), and frontoparietal (FPN). Upon visual inspection, components overlapping subcortical (SC) and cerebellar (CB) brain regions were assigned to these groups for subsequent analyses. The SC cluster contained only one component which corresponded to the bilateral basal ganglia (BG). Thus, this component will henceforth be referred to as "BG". Additional details on component rating are provided in SM section 1.3.

Functional network analysis
Using the "temporal dFNC" toolbox of GIFT, participants' component time-courses underwent additional steps of nuisance regression including despiking, detrending, regression of realignment parameters and derivatives, and bandpass filtering (0.01-0.15 Hz) -in line with default settings of the toolbox. Static FC was computed using pairwise Pearson's correlations across all time-points. Dynamic FC was performed using sliding window correlations, as described by , with a window size of 22 TR and a slide length of 1TR. Raw correlation coefficients were Fisher-Z transformed. K-means clustering was employed to extract distinct connectivity states using the "cityblock" distance metric, as recommended for high-dimensional data (Aggarwal et al., 2001). Based on convergence between the elbow criterion and Dunn's Index (Dunn, 1974), we identified five states. For reproducibility purposes, the numerical identity of the clusters was recoded according to descending total occupancy across the sample of participants. For each participant, the median over all dFC windows spent in a state was computed for further analyses. Additional details are provided in SM sections 1.4-5.

Dynamic metrics
We calculated the following additional measures of state dynamics: mean dwell time (the mean number of windows spent in a state upon entering it), fraction time (i.e., fractional occupancy; the proportion of scan time spent in a state), transition frequency (the number of switches between each pair of states), and state stickiness (the number of times a participant remained in the same state over two consecutive windows). These dynamic metrics have previously been defined and used to characterize connectivity states in von Schwanenflug et al., 2022).
State-wise average connectivity and modularity were calculated as measures of global connectivity and topology, respectively. Modularity was calculated using the community Louvain algorithm of the Brain Connectivity Toolbox (2019-03-03) (Rubinov and Sporns, 2010). Average connectivity was calculated as the global mean connectivity strength over all windows spent in that state. Additional details are provided in SM section 1.6.

Across-state overall connectivity (ASOC)
To explore the influence of differences in fraction time on a summary measure of global brain connectivity, we computed an additional  "across-state overall connectivity" (ASOC) metric. This was calculated by averaging across a participant's dFC data from all windows, resulting in one grand-average value across all connections and states, as well as 7 intra-network ASOC values (not applicable to BG which had only one component), and 8 inter-network ASOC values. Inter-network ASOC values were computed as the overall connectivity between one RSN and the rest of the brain.

Analysis of group differences
Between-group differences for static and dynamic FC were computed using a "non-parametric" two-sample T-test that estimates the null distribution by permuting group labels (see https://version.aalto.fi/gi tlab/BML/bramila). For static FC and for each dFC state, differences were assessed between three group pairs: HC vs patients with EDSS ≤ 1, HC vs patients with EDSS ≥ 2, and patients with EDSS ≤ 1 vs patients with EDSS ≥ 2. Statistical significance was defined at an alpha level of 0.05 after false discovery rate (FDR) correction for multiple comparisons according to (Benjamini and Hochberg, 1995). For static FC, each participant had a vector of 1081 values (lower triangle of the symmetrical 47-by-47 component matrix). For dynamic FC, each participant had a vector of 1081 values for each of the five connectivity states. We used a mass univariate statistical approach to test the hypothesis that for each of the 1081 FC values, group 1 differed from group 2. For FDR correction in the case of static FC, we considered a family of univariate tests to be within a group pair (e.g., HC vs EDSS ≤ 1) and thus correction was applied over these 1081 tests. In the case of dynamic FC, we considered a family of univariate tests to be within a group pair and within a statethus FDR correction was applied over 1081 tests for each state and each between-group comparison.
Group differences in dynamic metrics were assessed using a withinstate approach. Linear regression was used to remove age-related variance. Using the model residuals, Kruskal-Wallis (KW) omnibus tests on the three groups were performed. If the null hypothesis was rejected, post-hoc pairwise comparisons were performed using Dunn's tests with FDR correction. In the case of temporal dynamic metrics, we considered a family of tests to be within a state and within a measure (e.g., dwell time) -thus, FDR correction was applied over three tests (e.g., state 1 dwell time across HC, patients with EDSS ≤ 1, and patients with EDSS ≥ 2).
To test group differences in ASOC, an identical approach to the dynamic metric analyses was applied. In the case of group effects on internetwork ASOC (e.g., between the FPN and the rest of the brain), an additional analysis was performed between each pair of RSNs to identify which connections were driving the effect.

Correlation analyses
Group differences in domain clinical outcome scores were evaluated using non-parametric permutation tests. FDR correction was applied over seven hypothesis tests, in that patients with mild to moderate disability would differ from patients with no disability in each of the several clinical domains. Spearman's partial correlations were calculated between FC values and clinical scores, controlling for age. In the case of dynamic FC, a total of 37,835 correlations were tested (1081 FC values × 7 domains × 5 states). FDR-correction was applied within a state and within a domain, thus correcting over 1081 tests in each case. For correlations between dynamic metrics, ASOC, and clinical scores, Spearman's correlations were computed using regression model residuals.

Clinical characterization
Patients with mild to moderate disability (EDSS ≥ 2) showed significantly worse outcomes across all clinical domains compared to patients without disability (EDSS ≤ 1), including worse cognitive and motor performance, worse visual acuity, higher depressive and fatigue symptoms, as well as higher lesion load and total brain atrophy (Table 3, Fig. 1). Additionally, rank-based correlation analyses confirmed a significant positive relationship between disability ratings according to EDSS and the distribution of clinical outcome scores across patients (Supplementary Table 1, Supplementary Fig. 3), and this held true across all impairment domains.

Static FC group differences
Group mean sFC and between-group comparisons are shown in Fig. 2. Here, comparisons between patients with mild to moderate disability and matched HCs revealed decreased FC between the DMN and cerebellum in patients. Patients with mild to moderate disability showed a pronounced pattern of inter-network FC alterations compared to matched HCs, predominantly with FC increases between the DMN and FPN to other RSNs. In comparisons between both patient groups, patients with mild to moderate disability showed increased sFC between the FPN, SMN, and DMN and reduced connectivity between VIS-SMN and BG-dATT compared to patients without disability. However, these latter findings were not significant after FDR correction (Table 4). Fig. 3 shows the centroid position of each connectivity state and results of between-group comparisons. State 1 was the most frequently occurring and least globally connected state, while state 5 was the least frequently occurring state with highest global connectivity and lowest modularity. We observed an inverse relationship between the distributions of global average connectivity and modularity across states, as expected from previous network studies (Krohn et al., 2021). Additional state characteristics are shown in SM section 2.1.

Dynamic FC states and group differences in dynamic connectivity strength
State-wise between group comparisons revealed a widespread pattern of FC alterations. Compared to HCs, patients with EDSS ≤ 1 had altered FC in three statespredominantly involving regions of the DMN  (Tables 5-7).

Dynamic metric group differences
No group effects were observed for average connectivity, modularity, or dwell time. In contrast, a significant group effect was present for state 5 stickiness (i.e., the number of times a participant remained in state 5 over two consecutive windows) as well as for transitions between states 3-5 and 4-5. Further analyses showed that this effect was due to differences in how often this state was visited. Accordingly, we observed a significant group effect on state 5 fraction time (i.e., the proportion of total scan duration a participant spent in state 5), such that patients with mild to moderate disability spent significantly more time in state 5 than both patients without disability and HCs (Table 8; Fig. 4).

Across-state overall connectivity
Among all participants, state 5 had the highest average connectivity. We also observed that patients with mild to moderate disability spent more time in this high-connectivity state compared to both patients without disability and HCs. Thus, we hypothesized that patients in the mild to moderate disability group would also display higher global average connectivity across all states. In other words, spending more time in the highest connectivity state would drive a parallel increase in connectivity when considering all states simultaneously. On a global scale, we observed a gradual increase in ASOC from HC to patients with EDSS ≤ 1 to patients with EDSS ≥ 2. However, this was limited to a statistical trend. On the network scale, a significant group effect was found on inter-network connectivity between the FPN and the rest of the brain, such that patients with EDSS ≥ 2 had significantly higher FPN inter-network ASOC than patients with EDSS ≤ 1 and a trend in the same direction compared to HC. We then further investigated which internetwork connections may drive this result. Fig. 5 shows that patients with mild to moderate disability had significantly higher FPN-SMN ASOC than patients with no disability and HCs, as well as significantly lower FPN-BG ASOC compared to patients with no disability and a similar trend compared to HCs (Table 9). Correlation analyses between EDSS scores and FPN-BG overall connectivity strength revealed a significant negative association (rho = -0.3844, p = < 0.001; Supplementary Fig. 5). Thus, an inverse relationship between disability and FPN-BG ASOC was observed on both the group-level and when considering the distribution of EDSS scores across the patient sample. Critically, FPN-BG connectivity was also negatively correlated with fatigue and motor impairment (Fig. 5b). Control analyses using static FC data to estimate overall connectivity corroborated these findings. However, in these static analyses, group effects on overall connectivity between the FPN and the rest of the brain were not significant. Additional details are provided in SM section 2.2. Fig. 6 shows the results of correlation analyses between dFC, dynamic metrics, and clinical outcome scores that persisted when controlling for age. Lesion load was negatively correlated with both dATT-vATT and VIS-SMN connectivity in state 2. In contrast, lesion load was positively correlated with DMN-VIS connectivity in state 2. Depression was positively correlated with both DMN-vATT and dATT-vATT connectivity in state 3. Finally, total brain atrophy was positively correlated with DMN-VIS connectivity in state 3. Correlation analyses between dynamic metrics and clinical scores revealed that depression severity was positively correlated with state 4 average connectivity, and brain atrophy was positively correlated with average connectivity of state 3.

Associations with clinical impairment
Importantly, correlations between static FC strength and clinical scores were not significant. Similarly, no relationship was found between static global average connectivity and any domain clinical score. However, static modularity was positively correlated with brain atrophy (rho = 0.293, p FDR = 0.037).

Fig. 1.
Between-group differences in clinical outcome z-scores across domains. Panels a-g: violin plots show distributions of clinical outcome scores in patients with mild to moderate disability (orange) and those without disability (yellow); Motor, vision, and cognitive symptoms domains are composite scores (see Table 2 for details on individual tests in each domain). Please note, all clinical outcome z-scores are coded such that higher values represent higher impairment (i.e., worse performance, higher symptom severity, or more atrophy). Patients with mild to moderate levels of disability in terms of Expanded disability status scale (EDSS) rating also showed higher scores in all other domains compared to patients with no disability. Panel h: correlation matrix where cells contain Spearman's correlation coefficients between each pair of clinical outcome scores, computed across all patients. Asterisks indicate significance level: * = p FDR < 0.05, ** = p FDR < 0.01, *** = p FDR < 0.001. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.
Whole-brain static functional connectivity matrices and between-group differences. Top: Mean static FC matrices computed by taking the average across all participants within a group. Pearson correlation coefficients were Fisher-Z transformed. Abbreviations on x-and y-axes indicate connections grouped by corresponding resting state networks. Matrix diagonals were zeroed for visualization purposes. Bottom: Between-group comparison results from the relevant permutationbased T-test. Group pairs used for comparison are denoted above each ring graph. Nodes and adjacent numbers identify the corresponding component (see supplementary material seciont 1.3 for detailed component information). Nodes are grouped and color coded according to RSNs. Edge width maps the magnitude and direction of T-statistic. Edges were thresholded at an uncorrected p-value of 0.001. Asterisks indicate tests with a false discovery rate adjusted p-value < 0.05.

Table 4
Between-group differences in component-wise static functional connectivity.

Discussion
In the present study, we show that functional brain dynamics in RRMS exhibit widespread alterations that vary with disability status and reflect clinical outcomes across multiple symptom domainseven in a sample of patients with early-stage MS who show low levels of disability overall. Categorizing patients based on EDSS ratings revealed consistently higher impairment across all clinical domains in patients with mild to moderate disability compared to those without disability, supporting the validity of our approach with respect to clinical outcomes. Despite the overall low disease severity and the early disease stage, we reveal distinct patterns of altered connectivity between patients with without clinical disability and those with mild to moderate disability. With our whole-brain dynamic FC approach, we observed   widespread bidirectional FC alterations across cortical and subcortical components in MS patients. While dynamic analyses revealed an intricate pattern of inter-and intra-network alterations, static analyses only detected inter-network changes. These observations suggest that sFC approaches may obscure transient intra-network connectivity changes that are observed only on the finer temporal scales of time-resolved analyses (i.e., seconds versus minutes). Furthermore, we observed significant relationships between functional connectivity, depressive symptoms, and structural brain measures in patients that remained undetected with a static approach. These findings highlight the advantages of time-resolved methods in studying functional network interactions and emphasize their role in understanding the complex relationships between functional brain dynamics, disability status, and clinical outcomes in MS.
Regarding the implicated brain systems, we found consistent implication of the frontoparietal and default mode networks in componentspecific FC differences between patients with mild to moderate disability and patients without disability in both static and dynamic analyses. In the case of static FC, patients with mild to moderate disability showed predominantly increased FC between FPN, DMN and several other networks compared to patients without disability. On the other hand, dynamic analyses revealed a widespread pattern of connectivity alterations between patient groups, such that FPN connections to other RSNs were predominantly increased and DMN inter-network connectivity strength decreased in patients with mild to moderate disability compared to those without clinical disability. Altered FC in the FPN and DMN have been reported in previous studies of patients with MS. For example, increased static FC of core FPN and DMN regions was associated with loss of cognitive efficiency (Hawellek et al., 2011). Relatedly, another recent study reported increased connectivity between DMN and FPN to the rest of the brain in cognitively impaired patients with MS (Meijer et al., 2017). Static FC results of the present study support these previous findings. However, age-controlled correlations between clinical scores and static FC did not yield any significant associations in this study. In the case of dynamic FC, two previous studies have applied a similar methodological approach to studying whole-brain connectivity states in MS as used here. Cognitively impaired patients with relapsing-remitting MS showed reduced dynamism of resting-state FC compared to cognitively preserved patients, specifically between subcortical and default mode networks (d'Ambrosio et al., 2020). Reduced dFC in sensorimotor, cerebellar and cognitive networks was also identified across the main subtypes of MS, with worsening FC abnormalities correlated to motor and cognitive impairment only apparent in progressive MS (Hidalgo de la Cruz et al., 2021). While these studies also implicated subcortical, default mode, and attention networks, they focused on cognitive impairment in RRMS and relationships between disability and symptom severity across a sample of patients with high clinical variability, respectively. Against this background, a major contribution of the present work is the identification of distinct alterations in FC dynamics between two groups of RRMS patients who are both in early stages of the disease and have overall low levels of disability. Despite a relatively narrow range of EDSS scores in the present population, we show that impairment in cognitive, behavioral, and neuropsychological domains as well as structural brain atrophy is associated with higher disability ratings and relates to several measures of functional brain dynamics, even early on in the disease.
We observed that increased depression symptoms were positively correlated with FC between the vATT-DMN and vATT-dATT in patients with MS. Importantly, the relevant DMN component corresponded to the bilateral hippocampus, which has been repeatedly implicated in depressive symptoms (Sheline et al., 2002) and shown to be affected in MS (Heine et al., 2020;. In a study of hippocampal neuroinflammation, connectivity, and depression symptoms in MS patients, both positive and negative relationships between hippocampal resting-state connectivity and depression scores were reported (Colasanti et al., 2016). While our results relating to hippocampal connectivity complement these findings, we furthermore identify a potential role of ventral and dorsal attention networks in MS-related depressive symptoms. This relationship between depression and FC increase is further supported by the positive association between depressive symptoms and global average connectivity in state 4. Taken together, these findings show that transient increases in FC, both globally and network-specific, are related to depressive symptom severity in patients with MS. The lack of this relationship in static connectivity analyses suggests that it may be specifically altered temporal brain dynamics that play a role in depressive symptoms in MS. Looking to evidence from studies of patients with major depression, alterations in functional dynamics across several RSNs, namely the default mode, frontoparietal, limbic, and salience networks have been consistently related to symptom severity in multiple domains of depression, including negative thinking and rumination, as well as overall disease severity (Kaiser et al., 2016;Marchitelli et al., 2022). Thus, these findings further point to an essential role of temporal dynamics in understanding the link between brain network behavior and clinical presentation.
Fatigue and motor impairment are common features of MS. Several studies have reported inverse relationships between DMN-BG FC alterations and fatigue in MS. With a seed-based sFC approach, a negative relationship between fatigue severity and FC of BG with typical DMN structures has been observed (Finke et al., 2015;Jaeger et al., 2019). More recently, a dynamic approach has similarly reported that lower FC between the BG and DMN was associated with greater fatigue in patients with MS (Tijhuis et al., 2021). Here, we likewise observe the involvement of the BG, but instead implicate its connections to the FPN. We found that decreased overall connectivity between the FPN and BG was systematically associated with higher motor performance impairment, higher fatigue severity scores, and higher EDSS ratings in our patient cohort. While the implication of BG dysfunction and atrophy is well-documented in MS, specific FPN-BG interactions in this population are less understood. In healthy participants, the role of BG in motor control is well-established (Brooks, 1995). Additionally, the FPN is typically regarded as a task-positive network involved in cognitive control and facilitating coordination between other brain networks (Marek and Dosenbach, 2018). With this context, the observed consistent relationship between FPN-BG connectivity and clinical impairment across multiple domains suggests a clinically relevant functional decoupling between the FPN and BG in MS, pointing to an exciting new brain target linking aberrant FC to clinical outcomes.

Limitations
Some limitations of the present work warrant mentioning. We analyzed a cross-sectional sample of data, which places an inherent limitation on investigating causal mechanisms of observed links between FC, disability, and clinical impairment. Additionally, neuropsychological and behavioral data were not sufficiently available for HCs and were thus not included in between-group analyses. This prevented us from making direct comparisons between patient subgroups and HCs when analyzing clinical outcome scores and limited the interpretation of impairment severity against normative values. Additionally, variance in current medication across the patient sample was not considered in statistical analyses. Thus, it is not possible to exclude a potential influence of medication on the observed differences in functional brain dynamics or clinical outcome scores. Finally, it is important to note that Patients with mild to moderate disability have higher across-state overall connectivity (ASOC) between the frontoparietal network (FPN) and the rest of the brain compared to patients with no disability; (Middle) Patients with mild to moderate disability have lower ASOC between the FPN and basal ganglia (BG) compared to patients with no disability, with a similar trend compared to HCs; (Right) Patients with mild to moderate disability have higher ASOC between the SMN) compared to both patients with no demonstrated disability and HC. Panel b: FPN overall connectivity with the BG is negatively correlated with fatigue (left) and motor impairment (right); Asterisks denote significance level: * = p FDR < 0.05.

Table 9
Between-group differences in frontoparietal network across-state overall connectivity. while clustering-based analysis of brain states represents one of the most well-studied approaches to describe functional brain dynamics Calhoun et al., 2001;Calhoun et al., 2014), other approaches such as eigenvector centrality or whole-brain modelling are increasingly explored and might yield additional information in describing functional alterations in early MS. Finally, here we discretized the distribution of EDSS scores across the sample to obtain one group of patients with mild to moderate disability and one group of patients without disability. While further work is necessary to obtain more continuous accounts of the link between FC aberrations and clinical impairment, we also show in rank-based correlation analyses that the association between multi-domain clinical outcomes and EDSS scores holds on an ordinal scale.

Conclusions
These results show that connectivity alterations in early MS are not homogeneous but depend on temporal network dynamics and level of disability, even in a sample of patients with comparably low overall disability status. With multi-domain clinical outcome scores, we show thateven in early stages of the diseaseit is possible to identify significant differences in depressive and fatigue symptom severity, cognitive and motor performance, visual acuity, lesion load, and total brain atrophy between patients with no clinical disability and those with mild to moderate levels of disability according to EDSS. Despite the overall low disability in this population, we reveal distinct connectivity alterations between patient groups that are meaningful for clinical outcomes. Notably, the FPN was consistently implicated in both FC group differences and clinical associations to fatigue and motor performance impairment. Finally, dynamic analyses revealed both FC alterations and behavioral correlations that remained undetected in static analyses, highlighting that temporal brain dynamics are important to help disentangle the relationship between functional alterations, disability, and clinical outcomes in early MS.  . Partial correlations between dynamic functional connectivity, dynamic metrics, and clinical outcome z-scores. Panel a: Whole-brain dynamic functional connectivity values across patients were correlated with domain impairment scores. X-axes of subplots denote the resting-state network assignment, component numbers, and state in which the corresponding functional connection was significantly related to patient impairment indicated by the y-axis labels. Panel b: Dynamic metrics (x-axes) were correlated with domain impairment scores (y-axes). All correlations were performed using Spearman's partial correlations controlling for age, and false discovery rate (FDR) correction was applied within state and domain. Values plotted correspond to variable residuals after rank-based age regression.