Physiologically motivated multiplex Kuramoto model describes phase diagram of cortical activity

We derive a two-layer multiplex Kuramoto model from Wilson-Cowan type physiological equations that describe neural activity on a network of interconnected cortical regions. This is mathematically possible due to the existence of a unique, stable limit cycle, weak coupling, and inhibitory synaptic time delays. We study the phase diagram of this model numerically as a function of the inter-regional connection strength that is related to cerebral blood flow, and a phase shift parameter that is associated with synaptic GABA concentrations. We find three macroscopic phases of cortical activity: background activity (unsynchronized oscillations), epileptiform activity (highly synchronized oscillations) and resting-state activity (synchronized clusters/chaotic behaviour). Previous network models could hitherto not explain the existence of all three phases. We further observe a shift of the average oscillation frequency towards lower values together with the appearance of coherent slow oscillations at the transition from resting-state to epileptiform activity. This observation is fully in line with experimental data and could explain the influence of GABAergic drugs both on gamma oscillations and epileptic states. Compared to previous models for gamma oscillations and resting-state activity, the multiplex Kuramoto model not only provides a unifying framework, but also has a direct connection to measurable physiological parameters.

Fast electrochemical processes taking place on a complicated cytoarchitectural network structure render the human brain a highly complex dynamical system. Brain activity, as measured directly via EEG or MEG, or indirectly by means of MRI recordings, reveals characteristic macroscopic patterns such as oscillations in various frequency bands 1 , synchronization 2-4 , or chaotic dynamics 5 . Generally it is believed that macroscopic activity (involving 10 10 is closely related to high-level functions such as cognition, attention, memory or task execution. To understand the mechanisms of this correspondence, both the overall network structure of the brain and the local properties of neural populations have to be taken into account 4,6 . Regarding the latter, neural inhibition seems to be essential for cortical processing 7 . Two physiological phenomena received much attention lately in terms of a mathematical understanding: resting-state activity, and gamma oscillations. Resting-state activity is spontaneous, highly structured activity of the brain during rest, and can be described in terms of networks of simultaneously active brain regions 8,9 . Models of resting-state networks often rely on anatomical networks derived from histological or imaging data, and on local interactions between populations of excitatory and inhibitory neurons [10][11][12][13] . Oscillatory neural activity in the gamma range (30 100 − Hz) is potentially related to consciousness and the binding problem although its precise function remains unclear 14 . To understand the origin of gamma oscillations, two mechanisms have been proposed 15 . One describes interactions between inhibitory neurons together with an external driving force 16,17 . The other mechanism is based on excitatory-inhibitory coupling with synaptic time delays [18][19][20] . The relation of gamma oscillations and inhibition is experimentally well established. In mice 21 , rats 22 and humans 3,23 , a decrease of GABA-concentrations (gamma-aminobutyric acid is the main inhibitory neurotransmitter in mammals) is accompanied by a strong attenuation of the gamma frequency band and sometimes by epileptiform activity.
Many existing network models for resting-state activity and gamma oscillations are based on single-neuron local dynamics 10,11,[16][17][18][19] . Since experimentally observed resting-state networks comprise individual regions containing about 10 8 to 10 9 individual neurons, we believe that a local description in terms of Wilson-Cowan equations is an attractive alternative. The subject of multiplex networks received recent attention with applications reaching from social and technological systems to economy and evolutionary games 24,25 .
In this work we derive a simple two-layer multiplex model from classical physiological equations that is able to capture the main features of cortical activity such as oscillations, synchronization and chaotic dynamics. This model unifies the roles of neural network topology, synaptic time delays, and excitation/ inhibition. It provides a closed framework for simultaneously understanding the origin of resting-state activity and gamma oscillations.

Results
Derivation of the multiplex Kuramoto model. We consider N cortical regions indexed by Fig. 1a. Each region is populated by ensembles of excitatory and inhibitory neurons (e.g. pyramidal cells and interneurons). We define the activity level of a region i as the fraction of firing excitatory (inhibitory) neurons of the total number of excitatory (inhibitory) neurons in that region at a unit time interval, and denote it by x i (y i ). Neglecting for the moment interactions between different regions, we assume that individual cortical regions obey the Wilson-Cowan-type dynamics 20 Here a ij , b ij , c ij , d ij are positive synaptic coefficients linking regions i and j. τ accounts for transmission delays at inhibitory synapses (not to be confused with axonal conduction delays). In physiology, τ can be altered by changing the synaptic concentration of GABA 18,19,21,22 . In the present model, we assume that τ is proportional to the average synaptic GABA concentration in the brain. To derive a multiplex Kuramoto model (MKM) from Eqs. (1) and (2), we make the following three assumptions: Figure 1. Schematic illustration of the model setup. a The cortical surface is divided into N macroscopic regions. Every region i (blue) comprises excitatory (green) and inhibitory (red) neural populations with activity levels x t i ( ) and y t i ( ), respectively. Activity levels quantify the ratio of firing neurons in the region at time t. b Region i (left) receives excitatory (green arrows) and inhibitory (red arrows) inputs plus selffeedback (blue arrows). Inputs from adjacent regions j (right) are weak (dashed arrows) or very weak (dotted arrow). c Variable transformation from activity variables x t i ( ) and y t i ( ) to phase deviation variables t i φ ε ( ). On a limit cycle g, ε-perturbations of the x y ( , )-dynamics at time t induce the phase deviations t φ ε ( ).
Scientific RepoRts | 5:10015 | DOi: 10.1038/srep10015 (i) Homogeneity. Cortical regions exhibit nearly identical dynamical behavior. We therefore assume the following parameters to be constant across regions, for all i, up to small perturbations, denoted by ε.
(ii) Stable local oscillations. We choose the parameters a b c d x y ρ ρ ( , , , , ( ) such that each uncoupled system Eq. (1), under the assumption given in Eq. (3), has a unique exponentially stable limit cycle 2  g ∈ . As a consequence, after a transient time solutions of Eq. (1) can be written as is an arbitrary solution of Eq. (1) on g 26,27 . i ϕ accounts for specific initial values. Let T denote the period of g . We assume that the frequency T 1/ lies in the physiological gamma range. (iii) Weak coupling. Interactions between adjacent regions are weak, and inhibitory-inhibitory interactions are very weak in the sense that, for all i j ≠ . These assumptions are justified because the number of synaptic connections within a cortical region is much larger than between regions, and excitatory neurons outnumber inhibitory neurons by approximately one order of magnitude 28,29 . Figure 1b summarizes the connectivity structure between regions i and j. Region i receives excitatory (green arrows) and inhibitory (red arrows) inputs plus feedback (blue arrows), magnitudes are indicated by the arrow labels. For the sake of clarity, arrows representing inputs of magnitude  1 ( ),  ε ( ) and  2 ε ( ) are drawn in continuous, dashed and dotted style, respectively. Under these assumptions, the system Eq.
describes the deviations from the uncoupled phases t T t 2 i i θ π ϕ ( ) ≡ / ( + ) that are associated with solutions of the uncoupled system Eq. (4). Accordingly, d dt i φ / describes the deviations from the uncoupled oscillation frequency T 2π/ . Time t has been rescaled, see SI. A ij is the adjacency matrix of the excitatory-excitatory interaction network as defined in Eq. (5), and A ij δ ( ) is a linear combination of the adjacency matrices B ij and C ij . A ij δ ( ) accounts for the interaction between excitatory and inhibitory populations, see SI.
, are the corresponding average degrees. δ is a phase shift parameter related to the time delay τ via T 2 δ π τ ≡ ( / ) . K is a global coupling constant that we assume to be proportional to the cerebral blood flow. This is reasonable because the latter is strongly correlated with the connection strengths of functional networks reconstructed in magnetic resonance imaging 30 . i ω , the so-called natural frequencies of the MKM, are the constant contribution to the frequency deviations d dt i φ / . We take i ω from a symmetric, unimodal random distribution g ω ( ), with mean 0 ω . Since the 1-parameter family of rotating-frame transformations − , leave Eq. (6) invariant for any Ω, without loss of generality we assume, 0 Physiological processes changing δ, K, and i ω occur on a much slower timescale than neural activity.
It is known that weakly coupled, nearly identical limit-cycle oscillators can be described in terms of phase variables 27,[31][32][33] . However, in terms of the new variables N 1 φ φ ( … ), interactions between cortical regions take place on two independent layers representing excitatory-excitatory and excitatory-inhibitory coupling, respectively, and the complicated connectivity structure of Fig. 1b reduces to a simple two-layer multiplex structure. Figure 1c shows the variable transformation from activity variables x and y, to the phase variable φ for any cortical region. In the unperturbed case, 0 ε = , the limit cycle g is parametrized by g ( ) = ( ( ), ( )) t x t y t . Since g is exponentially stable, ε-perturbations of activity dynamics ( ( ), ( ))  t x t y t lead to phase deviations g π φ ε ( + /( ) ( )) t T t 2 . For 0 δ= we recover the Kuramoto model on a single network, see SI and 34-39 . Order parameters. We characterize solutions of Eq. (6) by the following order parameters: Synchronization. We define the order parameter [32][33][34]40 It takes values between 0 (no synchronization) and 1 (full synchronization) 27 . Let r denote its time average ≡ ( ) r r t t .
Chaotic dynamics. The instantaneous largest Lyapunov exponent is given by Numerical simulation of the model. Synchronization. We find that synchronization r depends on the coupling strength K, and phase shift δ, Fig. 2a. For 0 δ= , we expect (see SI) a transition from an unsynchronized to a synchronized state at a critical value K 0 74 c ≈ . , which is confirmed by our simulations, Fig. 2a. With 0 δ > , stronger coupling K is required for this transition to occur. Above a value of approximately 2 δ π = / , a global synchronized state ceases to exist.
Chaotic dynamics. Above the synchronization threshold, K K c > , synchronization and chaotic dynamics are mutually exclusive, see  = . , at 0 δ= , a synchronization peak appears close to frequency zero. With increasing δ, this peak moves towards increasingly negative values, until 3 8 δ π ≈ / . Between 3 8 δ π ≈ / and 5 8 π/ , the distribution is rapidly becoming broader and shifts towards positive values. After reaching a maximum at 3 4 δ π ≈ / , it is finally centered around zero again, Fig. 3b. K 5 = , is similar, however larger positive and negative values for d dt i φ / occur, Fig. 3c. Figure 2d shows the average frequency deviation Ω as a function of K and δ. As expected (see SI), we find frequency suppression associated with synchronization in the region of large K and small δ, but also for large K and intermediate δ. For fixed K, maximal frequency suppression occurs at 3 8 δ π ≈ / . For large K and large δ (chaotic region) we find slightly positive Ω.

Robustness issues
Homogeneity. The derivation of the MKM is based on three key assumptions, see Eqs.

Stability of local oscillations and weak coupling.
In contrast, both the exponential stability of the limit cycles and the weak coupling assumption, Eq. (5), are strictly necessary for the derivation of the MKM, since they allow for a dimensional reduction from activity-to phase deviation variables (see SI). If the dimensional reduction can not be carried through, the full system Eq. (1) with Eq. (2) has to be studied, whose properties are much harder to access.
Numerical simulation. We tested the model for robustness with respect to the particular choice of parameters. As suggested by various brain atlases and cortical parcellation schemes, a number of N 100 150 ≈ − cortical regions seems reasonable 41,42 . We tested up to N 500 = and found no deviations from the presented qualitative picture. For the link density p, we find that as long as it exceeds the percolation threshold, p N N log > ( )/ , differences in simulations are marginal. Finally, we observe that like in the original Kuramoto model 32,33 , for different natural frequency distributions g ω ( ) the qualitative behaviour remains practically unchanged as long as g ω ( ) is unimodal and symmetric.

Discussion
In several variants of single-layer Kuramoto models with a phase shift or a time delay, frequency suppression appears 37,39,43,44 . In addition 39 , mentions chaotic behavior. However, the existence of the phase diagram with the three distinct macroscopic phases can not be inferred from any of those models to the best of our knowledge. In Ref. 45, several modes of synchronization are reported for a Kuramoto model on two interconnected networks with an inter-network time delay. While this computational model does not exhibit chaotic or unsynchronized phases, it suggests that a more complicated network topology can lead to a deeper structure within the synchronized phase in Kuramoto-type models.
Since the present work emphasizes the derivation of the MKM and the study of its stationary properties, we did not investigate the details of the synchronization transition. In this context we mention that the emergence of synchronization follows different paths in different types of networks 46 . Further, if a correlation between natural frequencies i ω and network properties is assumed, explosive synchronization and hysteretic effects may appear 47 .
Summarizing, we can show mathematically that a set of weakly coupled Wilson-Cowan oscillators on a cortical network with a synaptic time delay between excitatory and inhibitory neural populations is identical to a simple Kuramoto-type phase model on a two-layer multiplex network. Numerical investigations of this model reveal the presence of three distinct macroscopic phases in the space of control parameters K (associated with cerebral blood flow) and δ (associated with synaptic GABA concentration). For couplings K K c < , activities of individual cortical regions show independent oscillatory behavior (unsynchronized). Frequencies are distributed symmetrically around an average frequency that we assume to be located in the physiological gamma range. This dynamical state corresponds to Scientific RepoRts | 5:10015 | DOi: 10.1038/srep10015 "background activity" of the brain. For K K c > , two phases are possible: for small δ, the system becomes synchronized, which corresponds to "epileptic seizure activity" in physiology. For large δ, synchronized activity only appears in clusters; the system is chaotic in general. We identify this phase with "resting-state activity" in the brain. An important property of the present model is that the average oscillation frequency is shifted towards lower values when crossing the boundary to the synchronized phase. This could explain the experimental fact 2,21-23 that a decrease of the GABA concentration in the resting-state both triggers the appearance of epileptiform slow waves and diminishes gamma activity in the brain.

Methods
Equation (6) is integrated with a standard 4th-order Runge-Kutta algorithm with 500 time steps of size dt 0 1 = .
. The system size is N 100 = , both layers are chosen to be Erdös-Rényi networks with p 0 06 = . . Natural frequencies i ω are taken from a standard normal distribution, initial phase deviations 0 i φ ( ) from the interval [0 2 ] π , . The first 250 time steps are discarded to exclude transient effects. For the remaining time steps, r, max λ , and Ω are evaluated. All results are averaged over 100 identical, independent runs with different realizations of the initial conditions.