Modulations of local synchrony over time lead to resting-state functional connectivity in a parsimonious large-scale brain model

Biophysical models of large-scale brain activity are a fundamental tool for understanding the mechanisms underlying the patterns observed with neuroimaging. These models combine a macroscopic description of the within- and between-ensemble dynamics of neurons within a single architecture. A challenge for these models is accounting for modulations of within-ensemble synchrony over time. Such modulations in local synchrony are fundamental for modeling behavioral tasks and resting-state activity. Another challenge comes from the difficulty in parametrizing large scale brain models which hinders researching principles related with between-ensembles differences. Here we derive a parsimonious large scale brain model that can describe fluctuations of local synchrony. Crucially, we do not reduce within-ensemble dynamics to macroscopic variables first, instead we consider within and between-ensemble interactions similarly while preserving their physiological differences. The dynamics of within-ensemble synchrony can be tuned with a parameter which manipulates local connectivity strength. We simulated resting-state static and time-resolved functional connectivity of alpha band envelopes in models with identical and dissimilar local connectivities. We show that functional connectivity emerges when there are high fluctuations of local and global synchrony simultaneously (i.e. metastable dynamics). We also show that for most ensembles, leaning towards local asynchrony or synchrony correlates with the functional connectivity with other ensembles, with the exception of some regions belonging to the default-mode network.

neural-mass models), this mean-field reduction provides an exact bridge between the microscopic properties of individual elements and their macroscopic collective properties. For this reason, neural mean-field reductions of this sort have been referred to as 'Next generation neural mass models' [2].
Here we outline briefly the mean-field reduction of Montbrio et al. [3]. This reduction uses a quadratic integrate-and-fire (QIF) neuron model as elementary unit. It does however not rely crucially on this particular neuron type since similar reductions have been done with a theta neuron model [1,2,4]. The macroscopic membrane potential and the firing-rate of an ensemble of neurons are obtained from a heterogeneous all-to-all-connected population of QIF neurons. This reduction is exact in the thermodynamic limit, i.e., when the number of QIF neurons is infinite. The method and the steps to obtain the macroscopic variables are not within the scope of this paper, but can be found in Montbrio et al. [3]. Next we summarize this procedure.
The QIF neuron is a class I neuron model derived to simplify the Hodgkin-Huxley neuronal model. Class I neuron models go from quiescence state to periodic firing when the input current is increased, and the neuron passes through saddle-node bifurcation. The state of a QIF neuron, j, is given by an ordinary differential equation of its membrane potential, V j .
The overdot represents the temporal derivative. I j represents a trans-membrane input current. When the membrane potential reaches a threshold V th , the neuron emits a spike and the membrane potential is reset to a resting value V r . The trans-membrane input current has the form , where η j is a constant external current, I(t) is a time varying input current common to all neurons, J is a synaptic weight, and s(t) is the mean synaptic activation.
[3] demonstrated that the macroscopic dynamics of a fully connected ensemble of neurons in the thermodynamic limit follow exactly two ordinary differential equations in which f represents the firing-rate and v the average membrane potential.
The heterogeneity of the external currents is represented by a Lorentzian distribution with halfwidth at half-maximum γ and maximum peak at η.

Derivation of the LSBM equations
For the sake of simplicity we obtain the KOP from a fully connected network of Kuramoto oscillators [7,8] instead of equations 3S, 4S, and 5S. The Kuramoto model has been widely studied theoretically and it has also been applied to many biological examples of collective synchronization [8]. In a Kuramoto model the synchronization is controlled by a coupling parameter. Natural frequencies are drawn form a probability distribution g(ω) and initial phases are drawn from a uniform distribution (r = 0). When the coupling is low, each oscillator follows its own natural frequency. When the coupling increases above a critical value the oscillators began to synchronize.
As the coupling continues increasing their synchrony increases until they reach full synchrony if the coupling continues increasing.
Next, we derive the equations of the LSBM. We consider a network with E ensembles labeled n = We are interested in describing each ensemble by the KOP introduced in equation S5. The KOP in an ensemble of oscillators can be defined as follows [7]: where r n and ψ n represent the average phase coherence and the mean phase in ensemble n respectively. Substituting equation S7 into equation S6, we can rewrite equation S6 in terms of the KOP as follows: Now we consider the case in which the number of oscillators in one ensemble is very large i.e.
N n →∞, and we define ρ n (θ n , ω n , t) as the probability density of oscillators in n with phase θ and frequency ω at time t. Since the number of oscillators on each ensemble remains unchanged, ρ n satisfies the local continuity equation, Next we consider the Fourier series of ρ n ρ n (θ,ω,t ) = g n (ω) where c.c. stands for complex conjugate. Following the Ott and Antonsen ansatzt [9], we consider a class of c (ω,t ) such that c n,k (ω,t ) = a n (ω,t ) k , where |a n (ω,t )| ≤ 1 . Then, inserting the series expansion a n (ω,t) into equation S9 we obtain ∂a n ∂t + iω n a n + 1 (S11) Letting the distribution of frequencies g n (ω) to be a Lorentzian distribution with center at Ωn and spread of Δ n , i.e.
we can calculate z n in the continuum limit as Thus, by evaluating equation S11 at ω = Ω n -Δ n , we can obtain the dynamics of z n aṡ Therefore the LSBM is defined by E complex delayed differential equations. Because of the first equality of equation S7, the LSBM with E equations can be converted to 2E real delayed differential equations with the formṡ r n = − Δ n r n + 1 − r n ψ n = Ω n + r n 2 + 1 whereṙ n and ψ n describe the temporal evolution of the KOP ( z n =r n e iψ n ) in terms of the level of synchrony and the mean phase in ensemble n respectively. The equations 1a and 1b in the main text are obtained by extracting the global coupling scaling factor G/E, the local coupling L n , and the anatomical network of neural fibers A np from the coupling matrix K np in equations S15a and S15b.