Alterations in Resting-State Functional Connectivity and Dynamics in Schizophrenia as a Result of Global not Local Processes

2.1 Patient Sample

The study sample was collected through the Center for Biomedical Research Excellence (COBRE) led by Dr. Vince Calhoun (more information here: http://fcon1000.projects.nitrc.org/indi/retro/cobre.html) and obtained from the SchizConnect database (http://schizconnect.org). This sample has previously been used by our group to explore structural deficits in patients with ScZ (Dimulescu et al. (2021)). From the sample of 43 patients and 43 healthy control participants, we excluded 5 patients due to missing resting-state functional MRI (rs-fMRI) data or artefacts/excessive motion identified during the pre-processing. We thus analyzed a final sample of 43 healthy control subjects and 38 patients with schizophrenia, which we will refer to as the COBRE sample. All patients were receiving antipsychotic medication (see Table 1). Symptom severity in patients was assessed using the Positive and Negative Syndrome Scale (PANSS) (Kay et al. (1989)). Written informed consent was obtained from all participants, and the study was reviewed and approved by the Institutional Review Board of the University of New Mexico.

2.2 Anatomical data

Data collection for the COBRE sample was performed using a Siemens Magnetom Trio 3T MR scanner. Structural images (high resolution T1-weighted) were acquired using a five-echo MPRAGE sequence with the following parameters: repetition time (TR) = 2530 ms; echo time (TE) = 1.64, 3.5, 5.36, 7.22, 9.08 ms; inversion time (TI) = 1200 ms; flip angle (FA) = 7^?; field of view (FOV) = 256 mm × 256 mm; matrix = 256 × 256; slice thickness = 1 mm; 192 sagittal slices. Diffusion tensor imaging (DTI) data were acquired using a single-shot EPI sequence with TR/TE = 9000/84 ms; FA = 90^?; FOV =256 mm × 256 mm; matrix = 128 × 128; slice thickness = 2 mm without gap; 72 axial slices; 30 non-collinear diffusion gradients (b = 800 s/mm²) and 5 non-diffusion-weighted images (b = 0 s/mm²) equally interspersed between the 30 gradient directions. For more information see also Cetin et al. (2014).

For model validation we additionally used a subset of 156 healthy participants from the human connectome project (HCP), which we will refer to as the HCP sample. The diffusion-weighted data were collected with multiband diffusion sequence (HCP version available at http://www.cmrr.umn.edu/multiband). Three different gradient tables are used, each with 90 diffusion weighting directions and six b = 0 acquisitions. More information can be found at https://www.humanconnectome.org/study/hcp-young-adult/document/1200-subjects-data-release.

2.3 Resting-state functional MRI data

COBRE data was acquired using single-shot full k-space echo-planar imaging (EPI) with ramp sampling correction using the intercomissural line (AC-PC) as a reference (TR: 2 s, TE: 29 ms, matrix size: 64x64, 32 slices, voxel size: 3x3x4 mm³). The resting-state scans were acquired in the axial plane with with an ascending slice order (multi slice method; interleaved). For more information see Aine et al. (2017). For the COBRE data set, we preprocessed the rsfMRI data using the FSL FEAT toolbox (Woolrich et al. (2001)). For each data set, we discarded the first five volumes. We corrected head motion using the FSL McFLIRT algorithm and subsequently high-pass filtered the data with a filter cutoff of 100 s. We linearly registered each functional image to the corresponding anatomical image of that subject using FLIRT. We then used the mean volume of the data to create a brain mask using BET. Using the ICA FIX FSL toolbox (Griffanti et al. (2014); Salimi-Khorshidi et al. (2014)), we conducted MELODIC ICA and removed artefactual components (motion, non-neuronal physiological artefacts, scanner artefacts, and other nuisance sources). Finally, we transformed the high-resolution mask volumes from MNI to individual subject functional space and extracted the average BOLD time courses for each cortical region in the AAL2 parcellation scheme using the fslmeants command from Fslutils.

Acquisition details for the functional MRI data from the HCP S1200 release can be found here: https://www.humanconnectome.org/study/hcp-young-adult/document/1200-subjects-data-release. For the HCP data set, we used the data preprocessed according to Glasser et al. (2013) and extracted the average BOLD time courses for each cortical region as described above.

2.4 Measures of connectivity and temporal dynamics

We used the average global brain connectivity (GBC) measure (Cole et al. (2010, 2011)) to assess the changes in connectivity strength. To assess alterations in temporal dynamics we used synchrony and metastability (Deco et al. (2017)). Because of the computational model being restricted to cortical areas, we also restricted our connectivity analysis to cortical areas. However, including subcortical regions did not substantially change the findings (see Supplementary Material).

Specifically, we define the functional connectivity matrix (FC) as the matrix of Pearson correlations of the BOLD signal between two brain areas over the whole time range of acquisition. From the FC matrices we calculate the global brain connectivity (GBC) of a single brain region i as follows (see also Cole et al. (2010, 2011)): where n is the number of regions. The average global GBC can then be defined as the average GBC over all cortical regions i. To calculate the average GBC for a functional subnetwork or generally a set of regions, one simply averages over the regional GBC values for the respective regions.

To assess the temporal dynamics of the functional networks, we used the Kuramoto order parameter as a measure of synchrony and its standard deviation as a measure of metastability, i.e. the variability of the states of phase configurations over time (see for example Deco et al. (2017)). Here the Kuramoto order parameter R(t) is defined as where again n is the number of regions and ϕ_k(t) is the instantaneous phase of the BOLD signal in region k. It measures the global level of synchronization of the BOLD signals from all regions, where a low level close to 0 reflects an almost uniform distribution of the signal phases, and a high value close to 1 reflects near equality of the signal phases. To calculate R, we band-pass filtered the signal in the narrowband 0.04-0.07 Hz (see Deco et al. (2017)) and then extracted the instantaneous phases of the signals at every time step using the Hilbert transform.

2.5 Computational network model

We use a whole-brain network model, where the connectivity, connection strength and delay between network nodes (i.e. brain regions) is derived from brain imaging data. As a model of single-node activity dynamics we employ a mean-field description of a spiking neural network of an excitatory and an inhibitory neural population, where the individual neurons are described by the adaptive exponential integrate-and-fire model (AdEx model; Brette and Gerstner (2005)), developed in our group (Augustin et al. (2017); Cakan and Obermayer (2020)). The following section describes the model in detail.

2.5.1 Single-Node model

A mean-field neural mass model based on a spiking network of coupled excitatory and inhibitory populations, the so-called ALN model (Augustin et al. (2017)), was implemented. The mean-field description offers a drastic speed-up of simulations on the order of about 4 orders of magnitude compared to the spiking model while still retaining its dynamical states and its biophysical parameters. The model has been extensively validated against simulations with the detailed spiking network and overall shows very good agreement (Cakan and Obermayer (2020)).

The mean-field reduction of the spiking neural network utilises the Fokker-Planck approach, i.e. the fact that in the limit of an infinite network size and under the assumption of a sparse, random connectivity, the distribution p(V) of the membrane potentials and the mean firing rate r_a of a population a, can be described by a Fokker-Planck equation (Brunel (2000)). However, to calculate the potential distribution a partial differential equation has to be solved, which is computationally costly. Therefore, the dynamics of a population is captured by a low-dimensional linear-nonlinear cascade model, and can be described by a set of ordinary differential equations (Fourcaud-Trocmé et al. (2003); Ostojic and Brunel (2011)). The mathematical derivation and the underlying assumptions have been detailed in (Augustin et al. (2017)), and we will only provide the final set of model equations in this manuscript.

A single network node in the whole-brain model is represented by the population activity of two interconnected neural populations, an excitatory population E and an inhibitory population I. The dynamics of the membrane currents of a population a ∈ {E, I}, are governed by the following equations: In the above equations μ_a describes the total mean membrane currents, the currents from synaptic activity, the currents from any sources of external input, the external noise input, τ_m the membrane time constant (calculated from the membrane capacitance C and the leak conductance g_L), and τ_s,b the synaptic time constant. Furthermore, is the variance of the membrane currents, and J_ab represent the maximum synaptic current when all synapses from population b to population a are active. The dynamics of the synapses are described by where represents the mean of the fraction of all active synapses, which lies in the range [0, 1], with the extreme cases being no active synapses and no inactive synapses, respectively. Furthermore, is the variance of s_ab.

The timescale τ_a = Φ_T (μ_a, σ_a) of the input-dependent adaptation, the average membrane potential , and the instantaneous population spike rate r_a = Φ_r(μ_a, σ_a) are computed every time step by means of precomputed transfer functions. The mean r_ab and the variance ρ_ab of the effective input rate from population b to population a can be described by and given a certain delay for the spike transmission d_ab. Here c_ab represent the amplitude of the post-synaptic current resulting from one individual spike (for s_ab = 0). Furthermore, K_gl scales the global coupling in the network, and C_ij and D_ij define the connection strengths and the connection delays between regions, estimated from the fibre count and fibre length matrices, respectively. Finally, δ_abE = 1 for a = b = E and 0 otherwise restricting coupling between regions to be exclusively from excitatory to excitatory populations.

The adaptive exponential integrate-and-fire model explicitly accounts for the evolution of a slow adaptation currents that represents both subthreshold and spike-triggered adaptation currents. The subthreshold adaptation current is described by the adaptation conductance α and the spike-triggered adaptation current is denoted by β. In the limit of infinite population sizes, an adiabatic approximation can be employed to describe the mean adaptation current in terms of the mean population firing rate. The mean adaptation current Ī_A can be understood as an inhibitory membrane current whose dynamics are governed by The individual populations a of a single region of the whole-brain network receive an external input current with a given mean and a standard deviation . This background input current can be thought to represent baseline input from extracortical areas in the brain. Additionally, the regions also receive a noise input current modelled as an Ornstein-Uhlenbeck process with a mean of 0 described by Here ξ(t) is a white noise process drawn from a normal distribution with a mean of 0 and a variance of 1. σ_ou determines the fluctuation amplitude of the noise around its mean.

To determine the mean external input to the E (μ_Eext) and I (μ_Iext) populations, the noise strength σ_ou, the subthreshold adaptation conductance α and spike-triggered adaptation increment β parameters for the model in the control condition, we used an evolutionary optimization procedure as described in Cakan et al. (2022). We compared the simulated BOLD FC to the empirical rsfMRI data. We initialized the algorithm with a random population of N_init = 160 individuals and repeated the evolutionary block with N_pop = 80 individuals for 100 generations. Initial parameter values were selected from a uniform distribution across the following intervals for the model parameters: μ_Eext ∈ [0.0, 4.0] mV/ms, μ_Iext ∈ [0.0, 4.0] mV/ms, σ_ou ∈ [0.0, 0.3], a ∈ [0.0, 40.0] nS, and b ∈ [0.0, 40.0] pA. The global coupling strength was set as in Figure 2 of Cakan et al. (2022). All other model parameters were set as given in Table 1 in Cakan et al. (2022) and they are summarised in Table 2.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2. Network parameters.

Overview of the different parameter values for the whole-brain model employed here.

• Download figure
• Open in new tab

Figure 1. Global differences in functional connectivity and temporal dynamics between healthy controls and ScZ patients.

a) Comparison of average GBC per participant for the two groups. Individual dots represent average GBC for one participant. The difference plot on the right shows the difference between the groups in terms of effect size. b) Histogram of region-wise GBC values for the two groups. The histogram displays the region-wise GBC data pooled for all participants in each group. c) Synchrony comparison between the two groups. Each dot represents the mean Kuramoto order parameter (a measure of synchrony) for one participant. The difference plot on the right shows the group difference in terms of effect size. d) Metastability comparison between the two groups. Each dot represents the metastability of one participant. The difference plot on the right shows the group difference in terms of effect size. e) Comparison of global brain connectivity for association areas (Asso. comprising: DMN, Cont, Sal/VAttn) and sensory areas (Sen. comprising: Sommot, Vis, DAttn). f) Comparison of global brain connectivity for the seven functional networks from Yeo et al. (2011): Somato-motor subnetwork (SomMot), Control subnetwork (Cont), Default mode subnetwork (Def), Salience/Ventral attention subnetwork (Sal/VAttn), Dorsal attention subnetwork (DAttn), Limbic subnetwork (Lim), Visual subnetwork (Vis).

• Download figure
• Open in new tab

Figure 2. Computational model

a) Modelling approach combining a model for the regional dynamics with anatomical input that defines the structural network. b) Average FC matrix for the COBRE sample c) Average FC matrix for the HCP sample d) Model FC matrix e) comparison of the correlation of empSC to empFC (blue) and the correlation of simFC and empFC (yellow) for the COBRE (left) and the for the HCP (right) data sets.

3.1 Demographic and clinical characteristics

The control and the patient group did not differ significantly in terms of age and gender (see Table 1). Patients also did not show a change in symptomatology or type/dose of antipsychotic medication during the three months before the assessment (for more details see Aine et al. (2017)).

3.2 Global differences in connectivity strength and temporal dynamics between ScZ patients and healthy controls

Global GBC was significantly reduced in patients with ScZ (effect size g = −0.65; see Figure 1 a) and Table 3). Comparing both groups a substantial shift from high GBC towards medium to low GBC values occurs in ScZ patients (Figure 1 b) and Table 3). Synchrony, as measured by the Kuramoto order parameter was lower in the patient group (effect size g = −0.44; see Figures 1 c) and Table 3). However, variability in synchrony, measured by metastability, did not significantly differ between groups (Figure 1 d) and Table 3).

Reductions of functional connectivity strengths affected all seven subnetworks (effect sizes ranging from g = −0.57 to g = −0.83; see Table 3), with the dorsal-attention, the somato-motor and the visual subnetworks showing the strongest effects (Figure 1 f)).

We further tested whether the GBC differences we found were specific to association areas as indicated by a previous study (Yang et al. (2016)). We grouped the default mode subnetwork, the control subnetwork and the ventral attention subnetwork together as the association areas and the somatomotor subnetwork, the visual subnetwork and the dorsal attention subnetwork as the sensory areas. We found reduced GBC for ScZ patients in both groupings, with the sensory areas showing an even stronger effect than the association areas (effect sizes g = −0.78 for sensory areas versus g = −0.61 for association areas, see Figure 1 e) and Table 3).

3.3 Mechanisms underlying connectivity and dynamics alterations

4.0.1 Global changes in connectivity and temporal dynamics

Evidence for large-scale dysconnectivity in functional networks has been accumulated over the last years in ScZ (Liang et al. (2006); Bluhm et al. (2007); Bassett et al. (2008); Liu et al. (2008); Bullmore and Sporns (2009)). However, it is still unclear, how these changes relate to changes on the microscopic level. To address this gap, we analysed resting-state fMRI data from healthy participants and patients with chronic ScZ. We identified a global reduction in functional connectivity that affected both sensory and association areas equally and that was present for all functional subnetworks together with a moderate decrease of temporal synchrony. Using a biophysical network model, we found that a decrease in global coupling or an increase in global noise levels could explain the connectivity reduction and the increase in synchrony best, whereas local changes to the glutamatergic or GABAergic system did not produce changes matching our experimental findings. However, both changes also yielded an increase in metastability in our model, which we did not find in the experimental data.

Our findings of reduced global brain connectivity are in line with previous research. For example, Lynall et al. (2010) and Bassett et al. (2012) both found significantly reduced global integration in patients with schizophrenia. However, we did not find stronger connectivity disturbances in association areas compared to sensory areas, as previously reported (Yang et al. (2016)).

Our analysis of the temporal dynamics of the activity, i.e. synchrony and metastability, revealed a decrease in synchrony but no change in metastability. Our finding of unchanged metastability is in line with previous findings of Lee et al. (2018) on the same dataset but in contrast to very recent work from Hancock et al. (2023a), proposing metastability as a candidate biomarker for schizophrenia. However, we have to note that Hancock et al. (2023a) introduced a new measure of metastability with increased sensitivity to detect the differences between healthy controls and ScZ patients. This new measure of metastability did not rely on predefined brain parcellations but rather flexibly defined recurring spatio-temporal modes, so-called ‘communities’ where single brain regions may be grouped into more than one community. As this approach was not applicable to our computational network model we did not employ it in our analysis. Overall, several different metastability measures have been proposed and have been applied in different contexts in neuroscience (Hancock et al. (2023b)).

4.0.2 Mechanistic explanations of global changes in ScZ

Reduced global coupling and increased global noise levels are in line with earlier modelling studies. For example, several studies, using both simple phase oscillator models and dynamic mean-field models, have shown that a decrease of global coupling compared to the best model fit to human resting-state data led to a decrease in connectivity and a more random, less integrated graph structure (Cabral et al. (2012a, 2013, 2012b)). Similar to the model presented here, the operating point is chosen close to a bifurcation point from a silent down state to a limit-cycle which produces oscillating activity. In this regime, both functional connectivity and temporal dynamics best match empirical data. Therefore, the reduced coupling or the increased global noise disturbs this specific state and thus reduces global connectivity, synchrony and more complex network properties.

Previous work on the effects of changes to the glutamatergic and GABAergic system has demonstrated profound alterations on the cortical microcircuit level. For example, numerous computational studies have shown that ScZ-associated changes on the microcircuit level can lead to substantial reductions in gamma power in auditory steady-state response tasks (Metzner et al. (2016); Metzner and Steuber (2021); Metzner et al. (2019); Vierling-Claassen et al. (2008)). Since local gamma oscillations have been hypothesized to at least partially determine the large-scale functional connectivity and temporal dynamics of resting-state activity (Cabral et al. (2014, 2022)), it seems surprising that changes to either of the systems did not produce changes in global brain connectivity in our model. One reason for the lack of impact of the changes might be that we applied them homogeneously. In the work presented here, we only varied glutamatergic or GABAergic strength globally, i.e. without any spatial heterogeneity. Therefore, it seems plausible that these changes disturbed the local, regional nodes all in a similar fashion and thereby did not substantially alter their interrelation, thus not changing global brain connectivity. Indeed, several studies have demonstrated that heterogeneous models of cortex, which explicitly incorporate regional differences in dynamics, match experimental resting-state functional connectivity more accurately (Demirtaş et al. (2019); Kong et al. (2021)). Importantly, these regional differences in dynamics covary with expression profiles for markers of glutamatergic and GABAergic neurotransmission and E-I balance (Burt et al. (2018); Demirtaş et al. (2019)). Therefore, a more detailed, heterogeneous model might be able to shed more light on the effect of E-I balance changes associated with ScZ on large-scale functional networks.

4.0.3 Limitations

The computational model that we have employed in this study, while generally showing a very good fit to the experimental data, is not fully biophysically realistic. Moreover, the model used an average connectome and was not able to provide subject-specific, individual results for each participant. Furthermore, the anatomic parcellation (AAL2 Rolls et al. (2015)) is relatively coarse-grained with a number of 80 cortical regions.

The ALN model that was used to simulate regional activity has been demonstrated to approximate cortical resting-state activity (Cakan and Obermayer (2020); Cakan et al. (2022)). However, it is restricted to the cortex. Including subcortical regions such as the thalamus into whole-brain models is still in its infancy and rarely goes beyond coupling a single cortical and thalamic region (e.g. Jajcay et al. (2022), but see Griffiths et al. (2020)).

The ALN model also presents a simplification of the regional circuitry as it approximates and neglects both the variability of cell types, especially the diversity of inhibitory interneurons, and the laminar structure of the cortex. Therefore, the inclusion of more detailed models of regional activity, both in terms of cell type diversity and of laminar structure and connectivity, seems likely to further our understanding of ScZ dysconnectivity and its underlying mechanisms.

Lastly, the regional ALN model we used had the same parameters regardless of the cortical region it represented, i.e. we implemented a homogeneous model in that respect. As already discussed above, cortical regions are known to differ in various important aspects, whose incorporation are likely to provide additional insight into the pathophysiology of schizophrenia.