Exploring global and local processes underlying alterations in resting-state functional connectivity and dynamics in schizophrenia

Introduction We examined changes in large-scale functional connectivity and temporal dynamics and their underlying mechanisms in schizophrenia (ScZ) through measurements of resting-state functional magnetic resonance imaging (rs-fMRI) data and computational modelling. Methods The rs-fMRI measurements from patients with chronic ScZ (n=38) and matched healthy controls (n=43), were obtained through the public schizConnect repository. Computational models were constructed based on diffusion-weighted MRI scans and fit to the experimental rs-fMRI data. Results We found decreased large-scale functional connectivity across sensory and association areas and for all functional subnetworks for the ScZ group. Additionally global synchrony was reduced in patients while metastability was unaltered. Perturbations of the computational model revealed that decreased global coupling and increased background noise levels both explained the experimentally found deficits better than local changes to the GABAergic or glutamatergic system. Discussion The current study suggests that large-scale alterations in ScZ are more likely the result of global rather than local network changes.


Introduction
ScZ is a severe mental disorder with a high burden of disease [Lopez and Murray (1); Charlson et al. (2)].However, the underlying mechanisms remain elusive.While no single brain area accounting for the heterogeneous symptom profiles has been identified, the notion that ScZ can be understood in terms of a general dysconnectivity has emerged [Friston et al. (3); Friston (4), Bullmore et al. (5); Pettersson-Yeo et al. (6)].
Experimental evidence for the dysconnection hypothesis comes from neuroimaging studies.Analyses of resting-state fMRI connectivity have shown widespread changes of functional connectivity.However, there is still a debate whether correlations of neural activity between regions are decreased [Liang et al. (7); Bluhm et al. (8)] or increased in ScZ [Zhou et al. (9)].There is also growing evidence for possible longitudinal changes of functional connectivity over the course of the disorder.Anticevic et al. (10) demonstrated that prefrontal cortical connectivity is increased in early-course ScZ while the opposite pattern was observed in chronic ScZ patients.Going beyond pairwise correlations between brain regions, graph theoretic measurements have identified reductions in integration, hierarchy, clustering, efficiency and small-worldness [Bassett et al. (11); Liu et al. (12); Bullmore and Sporns (13); Lynall et al. (14)].
Yet, the origin of functional dysconnectivity patterns in ScZ is still unclear.One hypothesis is that cellular and synaptic changes associated with ScZ disrupt local processing and thus impact on large-scale connectivity.Indeed changes at the microcircuit level have been identified in ScZ.Excitatory and inhibitory neurotransmission is disturbed, for example a reduced excitatory drive onto GABAergic inhibitory neurons [Chung et al. (15,16) and a decreased inhibitory output (Hashimoto et al. (17); Morris et al. (18); Moyer et al. (19)].Changes to the glutamatergic system, such as increased recurrent excitation, have been suggested to lead to deficits in large-scale connectivity with a gradient along the cortical hierarchy [Yang et al. (20)].
Computational models of large-scale brain circuits can be used to investigate dynamical circuit mechanisms linking local ScZassociated alterations to global changes in the functional organisation of the brain.Leveraging such computational models, studies have shown that decreases in global inter-regional connectivity strengths can lead to wide-spread functional disruptions [Cabral et al. (21)], increased global signal variance [Yang et al. (22)] and altered topological characteristics of functional brain networks (Cabral et al. (23,24) resembling ScZ.However, except for Yang et al. (22), these studies only investigated a global scaling of the inter-regional connectivity.Yang et al. (22) manipulated local and global neuronal coupling and demonstrated that both could increase signal variance as seen in ScZ but did not explore their potentially differential effects on large-scale functional connectivity.Thus, so far the effect of ScZ-associated local changes to glutamatergic and GABAergic neurotransmission and the effect of increased background noise on large-scale functional connectivity has not been explored.
To address this question, we quantified functional connectivity differences in a data set of healthy controls and chronic ScZ patients.We then implemented local microcircuit and global network parameter changes in a computational model of largescale cortical dynamics and compare the resulting connectivity changes to the experimental data.Furthermore, we also explored the temporal dynamics of the resting-state brain and characterised potential deficits in large-scale synchrony and metastability in ScZ patients and compared them to the different computational models, thus identifying mechanistic links underlying these changes.

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. (25)].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. (26)].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.
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-youngadult/document/1200-subjects-data-release.

Resting-state functional MRI data
COBRE data was acquired using single-shot full k-space echoplanar 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 3 ).The restingstate scans were acquired in the axial plane with with an ascending slice order (multi slice method; interleaved).For more information see Aine et al. (28).For the COBRE data set, we preprocessed the rsfMRI data using the FSL FEAT toolbox [Woolrich et al. (29)].For each data set, we discarded the first five volumes.We analyzed the relative mean framewise displacement as the root mean square (RMS) of the translation parameters.We found an average RMS of 0.15(± 0.09) for the control group and an RMS of 0.20 (± 0.10) for the patient group (t=2.0884,p=0.04).These results are in line with previous studies indicating that ScZ patients have higher framewise motion displacement than healthy controls [Guo et al. (30)].We thus 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. (31); Salimi-Khorshidi et al. (32)], 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. (33) and extracted the average BOLD time courses for each cortical region as described above.

Measures of connectivity and temporal dynamics
We used the average global brain connectivity (GBC) measure (Cole et al. (34,35) to assess the changes in connectivity strength.To assess alterations in temporal dynamics we used synchrony and metastability [Deco et al. (36)].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. (34,35): 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. (36)].Here the Kuramoto order parameter R(t) is defined as: where again n is the number of regions and f 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.07Hz[see Deco et al. (36)] and then extracted the instantaneous phases of the signals at every time step using the Hilbert transform.

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 (Figure 2).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 (37)], developed in our group [Augustin et al. (38); Cakan and Obermayer (39)].The following section describes the model in detail.

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. (38)], 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 (39)].
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 (40)].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. (41); Ostojic and Brunel (42)].The mathematical derivation and the underlying assumptions have been detailed in [Augustin et al. (38)], 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: m syn a = J aE s aE (t) + J aI s aI (t) In the above equations µ a describes the total mean membrane currents, m syn a the currents from synaptic activity, m ext a the currents from any sources of external input, m ou a the external noise input, t m the membrane time constant (calculated from the membrane capacitance C and the leak conductance g L ), and t s,b the synaptic time constant.Furthermore, s 2 a 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 s ab 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, s 2 s ab is the variance of s ab .The timescale t a = F T (µ a , s a ) of the input-dependent adaptation, the average membrane potential and the instantaneous population spike rate r a = F r (µ a , s a ) are computed every time step by means of precomputed transfer functions.The mean r ab and the variance r ab of the effective input rate from population b to population a can be described by: 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, d 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 a and the spike-triggered adaptation current is denoted by b.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 I A can be understood as an inhibitory membrane current whose dynamics are governed by: The individual populations a of a single region of the wholebrain network receive an external input current with a given mean m ext a and a standard deviation s ext a (t).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 m ou a (t) modelled as an Ornstein-Uhlenbeck process with a mean of 0 described by Here x(t) is a white noise process drawn from a normal distribution with a mean of 0 and a variance of 1. s 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 s ou , the subthreshold adaptation conductance a and spike-triggered adaptation increment b parameters for the model in the control condition, we used an evolutionary optimization procedure as described in Cakan et al. (43).We compared the simulated BOLD FC to the empirical rsfMRI data.We initialized the algorithm with a random population of

BOLD model
In order to compare the model output, i.e. the neural activity of the regions, to the BOLD signal of the rs-fMRI data, the firing rates of the excitatory population of each region had to be converted into model BOLD signal timecourses.Here, we used the well-established Balloon-Windkessel model [Friston et

Network connectivity
Structural images were preprocessed employing a semiautomatic pipeline implemented in the FSL toolbox (www.fmrib.ox.ac.uk/fsl,FMRIB, Oxford).For the anatomical T1weighted images we used the brain extraction toolbox (BET) in FSL to remove non-brain tissue and to generate the brain masks.After manual quality checks, 80 cortical regions were defined according to the automatic anatomical labelling (AAL2) atlas [Rolls et al. (47)].For the diffusion-weighted images, we performed a brain extraction as well and corrected the images for head motion and eddy current distortions afterwards.Probabilistic fibre tracking, using the Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques (BEDPOSTX) and the PROBTRACKX algorithms implemented in FSL [Behrens et al. (48)], was then used with 5,000 random seeds per voxel to extract individual connectomes.Since the tractography does not yield directionality information and the connectome matrices are non-symmetric, we explicitly enforced symmetry by averaging the entries from region i to region j and from region j to region i for all pairs of regions.Furthermore, we normalised each connectome by dividing each matrix entry by the maximum matrix entry, thus ensuring compatibility between participants.The resulting connectome then determines the relative coupling strength between regions in the above described computational whole-brain model.The fibre tracking also yielded matrix fibre lengths for each participant, which, when multiplied with the signal speed, determines the delay of signal propagation between any two regions in the model.

Modelling ScZ-associated changes
We implemented four different sets of parameter changes that are thought to represent the following four ScZ-associated alterations: 1) local GABAergic inhibition, 2) local glutamatergic  [Hashimoto et al. (17)], in cortical regions in ScZ.We implemented these changes as a reduction of the inhibitory weights J EI and J II in the ALN model of the regional dynamics.We varied the strength of the inhibition onto the excitatory population J EI and onto the inhibitory population J II simultaneously in the range from 100% to 60% in steps of 5%.
Next, we systematically reduced the glutamatergic, excitatory drive onto inhibitory neurons in our model of regional activity.These changes reflected the reduced and more varied colocalization of glutamatergic pre-and postsynaptic markers on PV interneurons in dorsolateral prefrontal cortex (DLPFC) (Chung et al. (15,16).Specifically, we reduced the excitatory weight onto inhibitory neurons J IE in the ALN model in a range from 100% to 60% in steps of 5%.
Global dysconnectivity might also be explained by a simple reduction of the global connectivity strength.Therefore, to test whether the differences we found experimentally could alternatively be explained by an overall network decoupling, we reduced the global coupling strength K gl in the range from 100% to 60% in steps of 5%.
Finally, the global alterations of functional connectivity might also be the result of an increase in background noise disrupting functional connectivity in the network [Winterer et  Consequently, we increased the global background noise s ou in a range from 100% to 140% in steps to 5%, to test whether a global increase in noise level can account for the connectivity differences found in the experimental data.

Simulation details
Simulations were implemented using the neurolib Python framework [Cakan et al. (52)].The differential equations of the model were solved numerically using an Euler forward scheme with a time step of 0.1ms.For all described simulations the duration was 70s and we discarded the transient response in the first 5 s before calculating any of the above described measures.To assess the robustness of our results, we created 40 virtual subjects by changing the seed for the random number generator underlying the Ornstein-Uhlenbeck noise process.These 40 virtual subjects were then kept fixed for all implemented changes allowing for a direct comparison to the default, 'healthy' condition.

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. (28)].

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 1A and Table 3).Comparing both groups a substantial shift from high GBC towards medium to low GBC values occurs in ScZ patients (Figure 1B and Table 3).Synchrony, as measured by the Kuramoto order parameter was lower in the patient group (effect size g = −0.44;see Figure 1C and Table 3).However, variability in synchrony, measured by metastability, did not significantly differ between groups (Figure 1D and Table 3).
Reductions of functional connectivity strengths affected all seven subnetworks (effect sizes ranging from g = −0.57to g = −0.83;see Table 3), with the dorsal-attention, the somato-motor a n d t h e v i s u a l s u b n e t w o r k s s h o w i n g t h e s t r o n g e s t effects (Figure 1F).
We further tested whether the GBC differences we found were specific to association areas as indicated by a previous study [Yang et al. (20)].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.61for association areas, see Figure 1E and Table 3).
Since our sample showed a significant difference in head motion between controls and patients, we investigated whether the changes in connectivity and dynamics were still present when applying a very strict threshold for head motion.Specifically, we had 4 control participants and 5 patients with a framewise displacement > 0.3.We repeated the analysis of FC and temporal dynamics after removing these subjects.Details of this analysis can be found in the Supplementary Table S3.Removal of the participants did not alter the results significantly with one exception.For the synchrony measure the mean difference and effect size both slightly decreased and did not reach statistical significance anymore.

Mechanisms underlying connectivity and dynamics alterations 3.3.1 Control model
We derived a model of healthy large-scale cortical activity that matched the behaviour of the control group data from the COBRE study well in terms of functional connectivity (Figures 2B, D).The correlation between simulated FC (simFC) and empirical FC (empFC) (r = 0.33 ± 0.09; Figure 2E) was higher than the correlation between empirical structural (empSC) and empFC (r = 0.19 ± 0.07; Figure 2E).
To further assert that the default model captures the restingstate functional connectivity of healthy subjects well, we also validated the model behaviour against a set of 156 subjects from the HCP S1200 release.Here, we also found a good fit for functional connectivity (Figures 2C, D).
Overall, the model functional connectivity correlated well with the empirical functional connectivity of individual HCP subjects (r = 0.43 ± 0.08; see Figure 2E).This correlation was again substantially higher than the correlation of structural connectivity and empirical functional connectivity (r = 0.20 ± 0.08; see Figure 2E).

Modelling ScZ-associated alterations
We systematically performed perturbations to four key aspects of the model that have been associated with schizophrenia: 1) local GABAergic inhibition, 2) local glutamatergic excitation of inhibitory cells, 3) global interregional coupling, and 4) global noise levels.
We found that changing the inhibitory weights (model perturbation 1) did not alter the global GBC and the GBCs for sensory and association areas significantly.Furthermore, the changes did not alter the synchrony and the metastability (see Supplementary Table S3).As for the local changes to the inhibitory system, changes to the glutamatergic excitatory drive to the inhibitory population (model perturbation 2) did not result in significant changes in GBC on all levels, synchrony and metastability (see Supplementary Table S4).A reduction of global coupling (model perturbation 3) resulted in a strong decrease in global brain connectivity as well as connectivity within the sensory and association systems (Table 4).Additionally, synchrony decreased strongly and metastability increased for larger reductions (Table 4).
An increase in noise levels (model perturbation 4) yielded a strong decrease in global brain connectivity as well as connectivity within the sensory and association systems, even stronger than for the global coupling changes (Table 5).Additionally, synchrony decreased strongly and metastability increased for larger noise strengths (Table 5).12); Bullmore and Sporns (13)].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. (14) and Bassett et al. (54) both found significantly reduced global integration in patients with schizophrenia.However, we did not find stronger Comparison of average global GBC, average GBC in sensory areas, average GBC in association areas, average synchrony, and average metastability with increased noise (from 105% to 140% in steps of 5%).Shown are the mean differences, i.e. the mean of the default condition minus the respective increased noise condition and in brackets the effect size (Hedge's g).The mean in each condition is calculated over the 40 virtual subjects.Significant differences, i.e. a permutation p value of< 0.001, are highlighted in bold.Permutation tests were performed using 5,000 permutations of labels.Comparison of average global GBC, average GBC in sensory areas, average GBC in association areas, average synchrony and average metastability for reduced global coupling (from 95% to 60% in steps of 5%).Shown are the mean differences, i.e. the mean of the default condition minus the respective reduced global coupling condition and in brackets the effect size (Hedge's g).The mean in each condition is calculated over the 40 virtual subjects.Significant differences, i.e. a permutation p value of< 0.001, are highlighted in bold.Permutation tests were performed using 5,000 permutations of labels.
connectivity disturbances in association areas compared to sensory areas, as previously reported [Yang et al. (20)].
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. (55) on the same dataset but in contrast to very recent work from Hancock et al. (56), proposing metastability as a candidate biomarker for schizophrenia.However, we have to note that Hancock et al. (56) 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. (57)].

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 (21,23,24).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  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. (62,63), 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 [Demirtasȩ t al. (64); Kong et al. (65)].Importantly, these regional differences in dynamics covary with expression profiles for markers of glutamatergic and GABAergic neurotransmission and E-I balance [Burt et al. (66); Demirtaşet al. (64)].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.

Limitations
Patients in the sample used in this analysis were on typical and atypical antipsychotic medication with a mean dosage of 396.26 (CPZ-equivalent dosage).Antipsychotic medication, however, is known to affect functional connectivity.For example, risperidone treatment has been found to lead to abnormal functional and structural connectivity in striatal areas, prefrontal cortex, and limbic system components Hu et al. (67,68).Furthermore, Wang et al. (69) found increased FC in the default mode network and decreased FC in the salience network after antipsychotic treatment.Therefore, we cannot rule out that the FC alterations identified in our analysis are not a result of ScZ pathophysiology but rather an effect of chronic treatment with antipsychotic medication.
Another limitation of the participant sample analysed here is its moderate size.For the group comparisons of GBC a post-hoc analysis of achieved power [performed with GPower 3.1 Faul et al. (70,71)] resulted in sufficient achieved power, however, group comparisons of temporal dynamics suffered from lower power (see Supplementary Material for more details).Therefore, replication of our findings in larger independent samples is warranted.
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. (47)] is relatively coarse-grained with a number of 80 cortical regions.As further validation, a replication of the analyses provided here using different, more fine-grained anatomic parcellations is warranted.
The ALN model that was used to simulate regional activity has been demonstrated to approximate cortical resting-state activity [Cakan and Obermayer (39); Cakan et al. (43)].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. (72), but see Griffiths et al. (73)].
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.

Conclusion
The current study provides further evidence of large-scale changes in connectivity and temporal dynamics in ScZ through the analysis of resting-state fMRI.Furthermore, through computational modelling, it provides novel evidence that these changes might be explained solely by global reductions in coupling or increases noise levels, although we cannot rule out that local effects also contribute significantly.These findings emphasize the importance of global alterations in ScZ and might have possible implications for the development of treatments.

1
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. (53): SomMot, Somato-motor subnetwork; Cont, Control subnetwork; Def, Default mode subnetwork; Sal/VAttn, Salience/Ventral attention subnetwork; DAttn, Dorsal attention subnetwork; Lim, Limbic subnetwork; Vis, Visual subnetwork.

2
FIGURE 2Computational 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.

TABLE 1
Demographics and clinical characteristics.
-Illness duration (y) -17.19(12.61)- ).All other model parameters were set as given in Table 1 in Cakan et al. (43) and they are summarised in Table 2.

TABLE 2
Network parameters.Overview of the different parameter values for the whole-brain model employed here.excitation of inhibitory cells, 3) global interregional coupling, and 4) global noise levels.First, we systematically reduced GABAergic inhibition in the model.Postmortem gene expression studies have robustly demonstrated reduced levels of parvalbumine (PV) and somatostatin [SST) expression in PV (Hashimoto et al. (17)] and SST [Morris et al. (18)] interneurons together with a reduction of GAD 65 and GAD 67

TABLE 3
Local and global group differences.Overview of the global and local differences in functional connectivity and temporal dynamics between the healthy control and the ScZ patient group.

TABLE 5
ScZ-associated changes of noise parameters.

TABLE 4
ScZ-associated changes of global coupling.