Explaining the heterogeneity of functional connectivity findings in multiple sclerosis: An empirically informed modeling study

Abstract To understand the heterogeneity of functional connectivity results reported in the literature, we analyzed the separate effects of grey and white matter damage on functional connectivity and networks in multiple sclerosis. For this, we employed a biophysical thalamo‐cortical model consisting of interconnected cortical and thalamic neuronal populations, informed and amended by empirical diffusion MRI tractography data, to simulate functional data that mimic neurophysiological signals. Grey matter degeneration was simulated by decreasing within population connections and white matter degeneration by lowering between population connections, based on lesion predilection sites in multiple sclerosis. For all simulations, functional connectivity and functional network organization are quantified by phase synchronization and network integration, respectively. Modeling results showed that both cortical and thalamic grey matter damage induced a global increase in functional connectivity, whereas white matter damage induced an initially increased connectivity followed by a global decrease. Both white and especially grey matter damage, however, induced a decrease in network integration. These empirically informed simulations show that specific topology and timing of structural damage are nontrivial aspects in explaining functional abnormalities in MS. Insufficient attention to these aspects likely explains contradictory findings in multiple sclerosis functional imaging studies so far.

and decreased connectivity can arise in the same set of patients in different disease stages, and how both can relate to worse cognitive outcomes at the same time. We postulate that the reason for such contradictory results can only be elucidated in a completely controlled system. Discrepancy between results on connectivity in MS may be driven by several factors (Stam 2014). Heterogeneity is one factor that could play a major role, both heterogeneity in demographic factors and disease-specific factors. For example, studies differ with respect to proportion of disease types in the population of interest, disease duration, female versus male ratios, age, physical status, but also in terms of MRI characteristics, for example, the amount of diffuse white matter damage or grey matter atrophy. Other more technical factors may also have contributed to divergent results on connectivity and activity in MS, for instance, suboptimal processing of fMRI data and suboptimal statistics using well-known functional imaging toolboxes (Eklund, Nichols, & Knutsson, 2016).
However, another more fundamental hurdle is a lack of understanding how functional connectivity changes in MS relate to structural abnormalities, such as thalamic or cortical atrophy, assumed to reflect loss of neuronal connections and neurons (Popescu et al., 2015). We hypothesize that understanding the relationship between structure and function in MS could elucidate how discrepant connectivity results may emerge in a disease with ubiquitous heterogeneity. Disentangling the separate effects of white matter (WM) and grey matter (GM) pathology on connectivity and networks from empirical data is currently challenging due to several factors such as lack of longitudinal data, limited sample sizes, heterogeneity in patient groups, use of different MRI scanners with different field strength and hardware, and differences in imaging pipelines. Biophysical models, mimicking electrophysiological signals, have been used successfully to understand pathophysiological processes in several other neurological diseases (Breakspear et al., 2006;Modolo et al., 2010;Moran et al., 2011;van Dellen et al., 2013).
These models allow for realistic longitudinal data simulations that may help to understand the effects of demyelinating and neurodegenerative MS pathology on functional connectivity and functional network organization. In other words, at this stage, these models are both an adequate and a necessary alternative for longitudinal empirical data and allow analyzing the influence of structural damage (in its various forms) on the functional network in a completely controlled and systematic way.
In this study, we employ a biophysical model, informed by empirical data, to investigate the separate effects of white matter, cortical, and thalamic degeneration on functional connectivity and network organization in MS. We hypothesize that the stage and amount of grey and white matter MS pathology can elucidate the contradictory results that were previously explained as "functional reorganization." 2 | M E TH ODS 2.1 | Cortico-thalamic mean field model Given the strong clinical relevance of cortico-thalamic involvement for cognition in MS (Houtchens et al., 2007;Schoonheim et al., 2015a), we employ a well-known cortico-thalamic mean field model in this study (Abeysuriya, Rennie, & Robinson, 2014;Robinson, Loxley, O'connor, & Rennie, 2001;Robinson, Rennie, Rowe, & O'Connor, 2004). This model is informed by empirical data and optimized to produce realistic power spectra mimicking electroencephalogram/magnetoencephalography recordings. It is also one of the few large-scale biophysical models that include the thalamus. In short, the model consists of units, where each unit comprises two cortical (excitatory and inhibitory) and two thalamic (relay nuclei and a reticular nucleus) populations ( Figure 1a). We denote these groups as a 2 e; i; r; s f g , where e; i; r; s denote excitatory, inhibitory, reticular, and relay, respectively.
For each group, a, and region j, the mean membrane potential is denoted by V a;j and the mean firing rate by Q a;j , which were interrelated by a sigmoid function: Here, Q max refers to the maximum firing rate in Hz, u is the mean firing threshold in mV, and r is the standard deviation of this threshold.
The mean membrane potential V a;j itself fluctuates under the influence of incoming firing input from other groups within the same population Q a;j À , a 2 e; i; r; s f g Þ and from other excitatory regions belonging to other populations, / e;n , outside its own region. The mean membrane potential V a;j can be considered as a low-pass filter. For every group, its dynamics can be described by Here a and b are constants, which are independent of time and correspond to the synaptic rise and decay rates in s 21 ; v aa denotes the synaptic densities between the population types and e corresponds to the global structural coupling strength between populations. The last term in Equation 2c corresponds to noise, which is defined as  (Gong et al., 2009), which is of the same size as previously described empirical networks in MS (Tewarie et al., 2013(Tewarie et al., , 2014a(Tewarie et al., , 2014bTewarie, van Dellen, Hillebrand, & Stam, 2015). Therefore, external firing input from other populations is mediated by the presence of an structural connection A jn , where A denotes a 78 3 78 adjacency matrix of a literature based structural network (Gong et al., 2009). Input from the thalamus to the cortex and vice versa is delayed by s ct , and cortico-cortical excitatory input between regions is also mediated with a delay s jn , depending on the Euclidian distance in the AAL atlas between region j and n. As incoming excitatory firing input from other populations / e;n is propagated over long range white matter tracts, its initial activity Q e;n is damped by the following expression The constant g refers to the cortical damping rate. Equation 3 usually has an extra term on the left hand side, which contains the differential Laplace operator r 2 . Similar to Robinson, Rennie, and Rowe (2002), we ignore the spatial dynamics by setting r 2 50 and therefore all observables in the model become independent of position in space.
Note that there is still spatial dependence based on the incoming connections over the white matter tracts in the network (Equation 2a). Values for all constants can be found in Supporting Information, Table S1 and are the same parameters as used to simulate resting-state activity in previous studies (Hindriks & van Putten 2012Robinson et al., 2002), with the only exception that we increased v ei to account for the extra excitatory input from the network. Simulations based on the cortico-thalamic mean-field model were performed using the Euler-Maruyama method with an integration time step of 1 3 10 24 . The first 2 s (20,000 samples) of the simulated data were discarded to exclude any non-oscillatory data. Only the time series of the 78 excitatory cortical populations / e;n t ð Þ were used as the model output to mimic MEG signals (Robinson et al., 2002), as it is believed that measurable and empirical neurophysiological signals outside the head mainly originate from excitatory pyramidal neurons (Hämäläinen, Hari, Ilmoniemi, Knuutila, & Lounasmaa, 1993). Note that the thalamic populations are required to generate realistic alpha band oscillations in a thalamocortical loop (Hindriks & van Putten, 2013); however, they are not used as model output to mimic MEG signals since they are usually difficult to detect using MEG. Thus, signals from excitatory cortical populations were used and filtered in the alpha band (8-13 Hz), which is the dominant frequency band of the model, but also the frequency band which most often shows cognitively relevant altered connectivity and disrupted network organization in MS (Cover et al., 2006;Schoonheim et al., 2013;Tewarie et al., 2013Tewarie et al., , 2014aTewarie et al., , 2014b. Analysis of the signals from the excitatory populations enables us to study a realistic system of coupled brain regions in terms of their activities and functional connectivity between them.

| Outcome measures
We define four outcome measures: mean activity, mean functional connectivity, network diameter, and leaf fraction. connectivity is defined as phase synchronization measured by the phase locking value (PLV) (Lachaux, Rodriguez, Martinerie, & Varela, 1999).
The PLV assumes that there is high connectivity when the phase difference of two signals is stable over a certain amount of time, and was used to estimate alpha band connectivity between all possible region pairs to obtain a weighted functional connectivity matrix. We subsequently average across regions to obtain a single average PLV value for each simulation.
In addition, we examined network organization corresponding to the functional connectivity matrix obtained in the previous paragraph.
The difference between functional connectivity and functional network organization is that functional connectivity refers to the strength of connections between brain regions, whereas functional network organization (topology) corresponds to the pattern of connections between regions in the brain. Specifically, we assess network organization by

| Experiments
Given the clinical and cognitive impact of white matter, cortical, and thalamic degeneration, we focus on the effect of these types of damage on connectivity and networks. As the initial parameters of the model were obtained after parameter estimation based on empirical data from healthy controls ), we consider these parameter settings as our initial working point. For each experiment described below, we simulate network activity followed by computing the outcome measures. All simulations were repeated 50 times, and the mean and standard deviations across 50 realizations for each iteration were visualized in a graph. We follow a three-step approach.

| Thalamic degeneration
Thalamic degeneration was modeled by linearly decreasing the strength of within thalamic unit connections. The reductions were informed by thalamic volumes measured in an empirical study in longstanding MS (n 5 202) (Daams et al., 2016). The maximum reduction was equal to the relative reduction (%) in thalamic volumes in MS patients compared to healthy controls (100% 2 (18.50/20.78 3 100%) 5 11%, see table 1 in a previous study) (Daams et al., 2016).
As within-unit connections in both the cortex and the thalamus correspond to neuronal connection densities , alteration of these parameters in experiments 2 and 3 correspond more to neuronal/axonal/synaptic loss (neurodegeneration/atrophy) within the cortex than demyelination, which is not modeled in this study. The inclusion of neuronal/axonal loss can be justified given findings from a recent study demonstrating that MRI-measured atrophy in MS is driven by neuronal/axonal loss (Popescu et al., 2015). Modeling of grey matter demyelination would require within-unit delays, which is computationally not straightforward and requires future studies. In addition, the stronger association between neurodegeneration and disability justifies focusing on that aspect of the disease in the current study rather than demyelination (Geurts & Barkhof 2008).

| White matter damage
The effect of white matter damage on functional activity, connectivity, and network organization is shown in Figure 2. Increasing white matter damage induces a monotonic decrease in mean neuronal activity, while functional connectivity changes are characterized by an inverted Ucurve, meaning that white matter damage initially causes an increase in connectivity, followed by a subsequent decrease. This inverted U-curve is explained by damage to long range tracts that initially leads to higher local connectivity (i.e., connectivity of local groups is not perturbed by distant groups), followed by a global collapse of the network associated with lower connectivity. The effect of white matter damage on network integration and segregation remains small during the initial stages, but subsequently features a sudden and steep change during the transitional period between increased and decreased connectivity. The increased MST diameter (Figure 2c) can be interpreted as a decrease in network integration, whereas the decrease in MST leaf fraction ( Figure   2d) can be interpreted as an increased network segregation. Figure 2c could give the impression that the initial stage of stable MST diameter values (during 0-10 iterations), is actually characterized by a decrease in MST diameter. However, note that there is still a large overlap between the standard deviations during this regime, and therefore this slight decrease cannot be considered as significant.

| Grey matter damage: Cortical damage
Simulations for the effect of cortical damage on functional activity, connectivity, and network organization are shown in Figure 3 (blue traces). Increasing cortical degeneration induces initially an increase in cortical neuronal activity (Figure 3a), followed by a plateau. Increasing cortical damage also brings about an increase in cortical functional connectivity ( Figure 3b), but note that functional connectivity hardly changes in the initial phases of cortical damage. However, the effect of cortical matter damage on activity and connectivity is very distinct compared to the effect of white matter damage on these functional outcome measures. This discrepancy is however not very evident with respect to functional network organization. For cortical damage we can observe that MST diameter increases with increasing damage, that is, the networks become less integrated. The difference with white matter damage is, however, that this increase follows very quickly after the initial stages of damage, in contrast to the effect of white matter damage on the MST diameter, which initially stays quite stable. MST leaf fraction shows a rapid decrease after a short stable phase. Again the sign of the change is similar to that of white matter damage, being again the difference that network segregation is less robust against cortical damage than white matter damage.

| Grey matter damage: Thalamic damage
Simulations for the effect of thalamic damage on functional activity, connectivity and networks are shown in Figure 3 (red traces). Similar to cortical damage, increasing thalamic degeneration induces an increase in neuronal activity, however, less pronounced than for cortical damage. Change in neuronal activity in the cortex during the course of thalamic damage follows a nonlinear effect with two maxima over time.
Thalamic degeneration also brings about an increase in mean cortical functional connectivity, which occurs earlier during the course than for cortical damage, but thalamic damage has a less strong effect on cortical functional connectivity than cortical damage. However, thalamic damage seems to induce a comparable effect on network integration (MST diameter) and segregation (MST leaf fraction) to cortical damage, with the difference of a steeper slope for thalamic damage and a lack of an initial stable phase.

| D I SCUSSION
In this study, we used biophysical modeling to investigate the separate effects of white matter, cortical and thalamic degeneration on functional connectivity/activity and network organization in MS. We observed a nonlinear effect of white matter damage on functional connectivity, in the form of an initially increased, but subsequently decreased overall brain connectivity. Grey matter damage only resulted in increased connectivity, especially for simulated cortical degeneration.
At the same time, both white and grey matter degeneration led to decreased network integration and segregation, albeit occurring at a faster pace for grey matter degeneration. Taken together, these findings could explain the heterogeneity of connectivity results in the literature.
Our simulations demonstrated that white matter damage can induce both an increase and decrease in functional connectivity depending on the amount of white matter damage (in terms of axonal/ neuronal damage). This inverted U-curve curve of functional connectivity during the course of increasing white matter damage may explain its inconsistent relationships with clinical outcome measures (Rocca et al., 2010(Rocca et al., , 2015bSbardella et al., 2016;Tona et al., 2014). This would especially be the case for cross-sectional studies, where patients with different amounts of white matter damage could be found on different points along the curve in Figure 2b. At the same time, functional network organization seems to be fairly robust against white matter damage up to a critical point after which the functional network undergoes a rapid disintegration.
In contrast to the white matter, damage of the gray matter induced functional network disintegration at an earlier stage, that is, less cortical damage is required to lead to disintegration of functional networks. Apparently, functional network organization is less robust for grey matter damage than for white matter damage. This could explain why cognitive and clinical measures correlate more strongly to grey matter than white matter alterations in MS. The discrepancy between damage within the nodes (i.e., atrophy) and damage between the nodes (i.e., structural connections) on network organization seems to be a generic principle in a diverse range of networks, which is related to the concept of degeneracy (Fornito, Zalesky, & Breakspear, 2015): due to the fact that the brain has many more structural connections than nodes (e.g., nodes are usually connected to several other nodes), network function can be retained in the presence of damage between the nodes by using one of the many alternative connection routes. The pool of alternatives for damage within the nodes (i.e., in the gray matter) is much smaller (i.e., there is less redundancy for grey matter than for white matter connections), which could explain its sensitivity for network disruption.  (Minagar et al., 2013), while cortical atrophy is thought to be more prominent and developing in patterns in later stages (Fisniku et al., 2008). Whether the supposed sequential occurrence of thalamic and cortical atrophy is driven by one of the changes remains to be investigated. In contrast to functional connectivity, the direction in which network integration and segregation changes, seems to be unambiguous for white and grey matter damage, indicating that it may be more fruitful to analyze network topology measures in empirical studies as a correlate for cognitive and clinical outcomes.
Some methodological aspects of our study warrant discussion. First, our results apply to resting-state functional networks and not task-based connectivity findings, which requires a different way of modeling in terms of external inputs that drive the dynamics. Second, all large-scale computational models lack regional specificity along white matter tracts.
Thus we cannot distinguish between lesional damage and diffuse white matter damage along a white matter tract. Finally, we did not investigate the effects of heterogeneity of conduction velocities along white matter tracts to simulate diffuse white matter demyelination, and did not use thalamo-cortical connectivity as an outcome measure, warranting future studies. Last, the empirical data of the different types of structural damage that has been used to inform the simulations differs in disease duration. This will have influenced the end points (last iteration) in the graphs (Figures 2 and 3), but not the shape of the curves.
In conclusion, we have demonstrated that joint white and grey matter damage can induce concurrent increases and decreases in functional connectivity, warranting caution when interpreting such changes as a compensatory mechanism. In addition, these concurrent increases and decreases in connectivity could also explain the reported heteroge-