Birth and destruction of collective oscillations in a network of two populations of coupled type 1 neurons
BBirth and destruction of collective oscillationsin a network of two populations of coupled type 1 neurons
Benjamin J¨uttner a , Christian Henriksen a , and Erik A. Martens a a Department of Applied Mathematics and Computer Science, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark (Dated: September 25, 2020)We study the macroscopic dynamics of large networks of excitable type 1 neurons composed of two popula-tions interacting with disparate but symmetric intra- and inter-population coupling strengths. This nonuniformcoupling scheme facilitates symmetric equilibria, where both populations either display identical firing activity,characterized by either quiescent or spiking behavior; or asymmetric equilibria, where the firing activity of onepopulation exhibits quiescent but the other exhibits spiking behavior. Oscillations in the firing rate are possibleif neurons emit pulses with non-zero width, but are otherwise quenched. Here, we explore how collective os-cillations emerge for two identical neuron populations in the limit of an infinite number of neurons. A detailedanalysis reveals how collective oscillations are born and destroyed in various bifurcation scenarios and how theyare organized around higher codimension bifurcation points. Since both symmetric and asymmetric equilibriadisplay bistable behavior, a large configuration space with steady and oscillatory behavior is available. Switch-ing between configurations of neural activity is relevant in functional processes such as working memory, andthe onset of collective oscillations in motor control.
PACS numbers: 05.45.-a, 05.45.Gg, 05.45.Xt, 02.30.Yy
The Theta neuron model [1] is the normal form for thesaddle-node-on-invariant cycle bifurcation, i.e., it representsthe dynamic behavior near the excitation threshold of type 1neurons, and it is equivalent to the quadratic integrate-and-fire neuron [2–4]. These two neuron models have attractedmuch interest based on recently developed dimensional re-duction techniques [5, 6], allowing for an exact descriptionof neuron ensembles in terms of macroscopic collective vari-ables [6, 7], for reviews see also [8, 9]. Such neuron popu-lations mimic densely connected neural masses in the brain.Collective oscillations arising in the brain are important forgenerating rhythms in the brain, e.g., for motor control [10]and breathing [11]. The combination of excitatory and in-hibitory neurons is a known prerequisite for the generationof collective rhythms such as gamma rhythms [12]. In thisstudy we pursue the mathematical question of how collectiverhythms may arise in an even simpler model composed of twopopulations of identical excitatory neurons with nonuniformcoupling, and what their bifurcations are.
I. INTRODUCTION
The brain is a complex network of networks with hierar-chical structure [13, 14], thus organizing neurons into neuralmasses, communities with high connectivity, structures whichmay interact with one another [14, 15] to solve cognitive func-tions [16] by displaying different individual collective dy-namic behaviors. A prominent collective behavior observedin the brain occurs when a group of neurons synchronizes andoscillates in unison [17, 18]. Synchrony has been associatedwith solving functional tasks including memory [19], com-putational functions [20], cognition [21], attention [20, 22],processing and routing of information [23–26], control of gaitand motion [10], or breathing [11].Neural masses with densely connected neurons are inter- connected and form networks of modular structure. An im-portant functional aspect in such networks are situations underwhich each population may assume different collective dy-namic behaviors, such as low or high synchrony, or low andhigh firing activity. Thus, a network of oscillator populationsmay exhibit a large configuration space with different syn-chronization patterns, as is also exemplified by chimera statesin Kuramoto oscillator networks, where one or several pop-ulations are synchronized and the other desynchronized [26–36]. The dynamics of such networks with multi-populationstructure and their configurations has been explored in thecontext of neuroscience [7, 37, 38], including memory re-call [39], information processing via self-induced stochasticresonance [40], and deep brain stimulation [41].Many studies concern the modeling of neuronal processesat the microscopic scale of individual neurons. However,the number of neurons in the brain is enormous, and, con-sequently, mathematical models of the brain are very high di-mensional, so that analyzing the collective dynamic behaviorof large neuronal assemblies poses a prohibitive challenge; acoarse-grained description of the dynamics at the macroscopiclevel is desirable. Recently developed mathematical methods,based on the Ott-Antonsen [5] and Watanabe-Strogatz reduc-tions [42, 43] allow for an exact dimensional reduction, whichapplies to phase oscillator networks with sinusoidal coupling,including variants of the Kuramoto model, the theta neuronmodel and the equivalent quadratic integrate-and-fire neuronmodel. Unlike heuristic models [44, 45], the resulting modelequations exactly describe the collective dynamics for eachpopulation, and – connecting the microscopic to the macro-scopic description – accurately capture microscopic propertiesof the underlying system [8, 9, 46].Collective oscillations in neural activity occurs over a broadrange of frequencies and across many brain regions [47].Prominent are gamma frequency oscillations, relevant in con-nection with cognitive tasks [48], neuronal diseases [18], mo- a r X i v : . [ n li n . AO ] S e p tor control [10] and breathing [11]. Such collective oscilla-tions are known to occur in neuron networks with excitatoryand inhibitory coupling [46, 49]. Network models with (sta-tistically) identical neurons emitting infinitely ’sharp’ signalpulses as represented by Dirac distributions do not permit col-lective oscillations [50]. Conversely, collective oscillatory be-havior is possible when the pulse width is non-zero [7, 51].We study a network composed of two populations of in-hibitory type 1 neurons with non-uniform (but symmetric)coupling, interacting through pulses with non-zero width. Weconsider the dynamics in the continuum limit of infinitelymany neurons, allowing us to use aforementioned dimen-sional reduction methods [5, 8]. Rather than aiming at a highlevel of biophysical realism, we wish to elucidate how col-lective oscillations may get born and destroyed in a simplesetup and to explore their related bifurcation scenarios. Eventhough the coupling is symmetric and neurons are statisticallyidentical, the resulting dynamic behavior is surprisingly com-plicated. The neuronal activity in each population may as-sume distinct levels, thus resulting in multistable configura-tions, in similarity to synchronization patterns as those ob-served in chimera states [27, 36], or (non-oscillatory) neuralstates reported for models of working memory [39]. In partic-ular, one observes a rich structure of bifurcations producingcollective limit cycle oscillations for which we provide a de-tailed bifurcation analysis.The article is structured as follows. In Sec. II, we intro-duce our model of two populations of Theta neurons and itsequivalent form of quadratic integrate-and-fire neurons. Weoutline how an exact description of the macroscopic dynam-ics for populations of infinitely many neurons is obtained viathe Ott-Antonsen method, and how firing rate equations forthe equivalent QIF neurons are derived via a conformal map-ping [6]. In Sec. III, we summarize the known dynamicalbehavior for a single population, which represents a limitingcase for two populations with vanishing inter-population cou-pling or uniform coupling. In Sec. IV, we perform a detailedanalysis by using numerical continuation methods via Mat-Cont [52], and explain the various bifurcation scenarios thatare possible. Finally, we sum our findings up and concludewith a discussion in Sec. V. II. MODELA. Network of Theta neurons
We consider a model of M = 2 populations of N inter-acting Theta neurons, where the phase θ σ,k ∈ T := R / π Z of the k th neuron belonging to population σ = 1 , evolvesaccording to ˙ θ σ,k := d θ σ,k d t = 1 − cos θ σ,k + (1 + cos θ σ,k )( η σ,k + I σ ) (1)with excitability η σ,k of oscillator k in population σ sampledfrom a from a Lorentzian distribution g σ ( η ) with mean ˆ η σ and width ∆ σ . The Theta neuron (1) is the normal form of the saddle-node on an invariant circle (SNIC) (or saddle-node-infinite period) bifurcation [53] and is a canonical type 1 neu-ron [1]. The dynamics are as follows. For η k + I < , astable and unstable fixed point occur on the phase circle T ;for η k + I = 0 , these fixed points coalesce in a saddle-nodebifurcation; for η k + I > , the flow on the circle results ina cyclic/periodic motion. If η k + I < , the Theta neuron issaid to be excitable : in the absence of perturbations, the phaserelaxes to the stable fixed point on the phase circle T ; how-ever, a perturbation may lead to a single spike (at θ σ,k = π )before returning to the stable fixed point. This could happenin at least two ways: a perturbation of the phase across theunstable fixed point (constituting a threshold) is possible, ifone considers that the Theta model derives from a higher di-mensional model [1] so that the circle is embedded in a higherdimensional space; alternatively, a very short-lived (time scaleof a single cycle) increase in I σ momentarily pushes the sys-tem across the bifurcation threshold η k + I = 0 . If η k + I > ,the neuron is firing (or excited), i.e., it spikes periodically.The input current may result from a variety of interactions,for an overview see [8, 50]. Here, we assume that the inputcurrent is given by I σ = M (cid:88) τ =1 κ στ N N (cid:88) l =1 P s ( θ τ,l ) (2)where adjacent neurons interact via pulses, which we choose P s ( θ ) = a s (1 − cos θ ) s , (3)originally adopted by Ariaratnam and Strogatz [54], withshape parameter s ∈ N , see also Fig. 1, and coupling strengths κ στ between populations σ and τ . The normalization constant a s = 2 s ( s !) / (2 s )! is defined so that (cid:82) π P s ( θ ) dθ = 2 π . s = = = = = π / π π / π θ Figure 1. Pulse shape for varying pulse shape parameter s . Thepulse converges to Dirac delta as s → ∞ . The case of M = 2 populations results in eight parame-ters (excluding the pulse shape parameter s ). To reduce theproblem to a manageable number of parameters, we make thefollowing assumptions: (i) the oscillator properties in popu-lations σ = 1 , are identical so that ˆ η = ˆ η =: ˆ η and ∆ = ∆ =: ∆ ; and (ii) the coupling is symmetric withrespect to identical intra- and inter-coupling strengths, i.e., κ = κ =: κ . and κ = κ =: aκ . Unless stated oth-erwise, we keep (ˆ η, ∆ , s ) fixed and consider ( κ, a ) the mainbifurcation parameters. B. Network of Quadratic Integrate and Fire neurons
An equivalent description of the Theta neuron is thequadratic integrate-and-fire (QIF) neuron via the transfor-mation into the membrane potential V ( θ ) = tan ( θ/ ∈ ( −∞ , ∞ ) . The model equations then become [8]dd t V σ,k := ˙ V σ,k = V σ,k + η σ,k + I σ , (4)where V σ,k := V ( θ σ,k ) . In this formulation, the neuron fires(emits a spike) when the voltage reaches V k ( t − ) = ∞ (infinite time). It is then reset to V k ( t + ) = −∞ . QIF neuronshave been widely used in neuroscientific modeling; see [4, 55]for a general introduction and [3, 56] for a few examples ofapplications of QIF neurons. C. Exact macroscopic description for infinite neurons
We consider (1) in the limit N → ∞ , which allows us toexpress the ensemble dynamics in terms of a continuous neu-ron density ρ σ ( θ, η, t ) governed by the continuity equation ∂∂t ρ σ + ∂∂θ ( vρ σ ) = 0 (5)where v σ = 1 − cos θ + (1 + cos θ ) (cid:18) η + (cid:90) ∞−∞ (cid:90) π I σ (cid:19) . (6)The Ott-Antonsen method [5, 7] facilitates an exact reductionof the microscopic dynamics in (1) to a low-dimensional de-scription of the macroscopic dynamics in terms of the com-plex order parameter of each population, Z σ ( t ) = R σ ( t ) e − i Φ σ ( t ) = (cid:90) π (cid:90) ∞−∞ e iθ ρ σ ( θ, η, t ) d η d θ. (7)The absolute value of the order parameter informs us of thelevel of phase synchronization of the neuron population: when | Z σ | ≈ , phases are spread over the circle T , whereas | Z σ | ≈ implies phase synchronization, i.e., phases areclosely spread around the phase of the order parameter givenby Φ σ = − arg ( Z σ ) . The collective dynamics of population σ = 1 , is then given by [7, 8] ˙ Z σ = − (cid:2) (∆ σ − i ˆ η σ − iI σ )(1 + Z σ ) + i (1 − Z σ ) (cid:3) . (8)These equations are closed by the input current [7, 8] I ( s ) σ = M (cid:88) τ =1 κ στ P ( s ) σ , (9) with the average output from all other neurons in the network, P ( s ) σ = a s (cid:32) C + s (cid:88) q =1 C q ( Z q + ¯ Z q ) (cid:33) , (10) C q = s (cid:88) k =0 k (cid:88) m =0 s !( − k δ k − m,q k m !( n − k )!( k − m )! . (11)For details on this reduction method and theory in general in-cluding applications in neuroscience, see [8].Two cases are of particular interest to us: pulse shape pa-rameter s = 1 and s = ∞ (impulsive coupling) for which wehave P (1) σ = 1 −
12 ( Z σ + ¯ Z σ ) , (12)and P ( ∞ ) σ = 1 − | Z σ | (1 + Z σ )(1 + ¯ Z σ ) , (13)respectively. D. Firing rate equations
The model (8) has an equivalent formulation in terms of av-erage firing rate r σ and average membrane potential v σ calledthe Firing Rate Equations (FRE) [6]. Indeed, changing vari-ables via the (anti)conformal mapping Z = (1 − ¯ W ) / (1 + ¯ W ) or W = (1 − ¯ Z ) / (1 + ¯ Z ) , (14)gives ˙ W σ = ∆ σ + i ˆ η σ − iW σ + iI ( s ) σ (15)and P (1) σ = 1 − − | W σ | (1 + W σ )(1 + ¯ W σ ) , (16) P ( ∞ ) σ = 12 ( W σ + ¯ W σ ) . (17)Writing W σ = πr σ + iv σ , (16) and (17) take the form P (1) σ = 2 π r σ + πr σ + v σ ( πr σ + 1) + v σ , P ( ∞ ) σ = πr σ , (18)for s = 1 and s = ∞ , respectively. Taking real and imaginarypart of (15) yields the firing rate equations ˙ r σ = ∆ σ π + 2 r σ v σ , (19) ˙ v σ = v σ − π r σ + ˆ η σ + I ( s ) σ , (20)where I ( s ) σ = (cid:80) Mτ =1 κ στ πr τ + πr τ + v τ ( πr τ + 1) + v τ , s = 1 π (cid:80) Mτ =1 κ στ r τ , s = ∞ . (21)The microscopic and macroscopic description are related asfollows. A single Theta neuron fires when its phase crosses θ = π ; accordingly, the average firing rate r σ ( t ) of the net-work at time t is defined as the flux through θ = π (or equiv-alently, the flux at v σ = ∞ ), see for instance [8]. III. DYNAMIC BEHAVIOR OF ONE POPULATION
The dynamic behavior for the case of M = 1 populationhas already been studied previously [7, 51]. We briefly reviewthe dynamics observed for this case as it is instructive for un-derstanding the dynamic and oscillatory behavior exhibited by M = 2 populations. For two parameter choices, the modelequations (1) for M = 2 effectively reduce to the dynamics ofa single population, M = 1 . Recall that the intra- and inter-coupling strengths among the two populations are given via κ = κ = κ and κ = κ = κa . Thus, when a = 1 , allneurons experience identical coupling strength so that the twopopulations act like a single population consisting of twice thenumber of neurons; on the other hand, when a = 0 , the twopopulations are decoupled so that each of the two populationsin separation effectively corresponds to a M = 1 system. Forbrevity, we drop σ in (1) and all related equations.The bifurcation diagrams in Fig. 2 report minima and max-ima for the firing rate r while varying coupling strength κ withparameter values s = 1 , ∆ = 0 . fixed, and ˆ η = − . or ˆ η = 0 . in panels a) and b), respectively. Solution branchessometimes appear very close to each other for the firing rate r , therefore it is instructive to also report the magnitude ofthe order parameter, | Z | , which is related to the firing rate r via the (anti)conformal mapping (14). Equilibria and local bi-furcations (saddle-node, Hopf) can be computed analyticallyfrom (19) and (20); limit cycles and other bifurcations werecomputed and continued numerically using Matlab and Mat-Cont software [52].We first consider the case of excitable neurons ( ˆ η < )in Fig. 2a). For small κ , we observe a set of stable equi-libria (stable nodes) with | Z | ≈ ; the related microscopicstates are non-oscillatory, i.e., most of the neurons are quies-cent (Q), and so their spiking activity is negligible, r ≈ .This branch of equilibria may undergo two saddle-node bifur-cations ( SN and SN ) which are connected by a branch ofsaddles. Equilibria to the right of SN (larger κ ) are stablespirals and correspond to spiking neurons (S) with larger fir-ing rate r > . As the coupling strength κ increases, higherlevels of synchrony, eventually getting close to | Z | = 1 , maybe achieved.For the case of spiking (firing) neurons ( ˆ η > ), the bifur-cation diagram in Fig. 2b) reveals a similar bifurcation struc-ture with two saddle-node bifurcations. However, for certainvalues of ˆ η , an even more complicated bifurcation scenario κ m i n t ( r ) , m a x t ( r ) m i n t ( | Z | ) , m a x t ( | Z | ) SN SN SN SN Quiescence (Q) S p i k i ng ( S ) a) 6-6 κ m i n t ( r ) , m a x t ( r ) m i n t ( | Z | ) , m a x t ( | Z | ) SN SN HBHCSN SN HBHC
Quiescence (Q) S p i k i ng ( S ) b)Figure 2. Bifurcation diagrams in κ for M = 1 population of Thetaneurons. We display solution branches by reporting maxima andminima in the firing rate r (black) and synchrony level | Z | (gray),respectively. Stable and unstable branches of equilibria have co-inciding minima/maxima and are shown solid and dashed, respec-tively; minima/maxima corresponding to limit cycle behavior arehighlighted in blue. Bifurcations that may occur are: saddle node( SN , SN ), Hopf ( HB ) and homoclinic ( HC ). Fixed parametersare ∆ = 0 . and pulse shape parameter s = 1 , while ˆ η = − . and ˆ η = 0 . in panels a) and b), respectively. is possible along the branch to the right of SN : a supercrit-ical Hopf bifurcation ( HB ) gives birth to limit cycles whichultimately are destroyed in a homoclinic bifurcation ( HC ).In between the values of ˆ η = − . and ˆ η = 0 . shown inFig. 2, two distinct bifurcations of codimension 2 occur: (i) SN and SN merge in a cusp point, and (ii) the bifurcationcurves SN , HB and HC meet in a Bogdanov-Takens point.The scenario in which limit cycles occur is characteristic forspiking neurons ( ˆ η > ) with inhibitory coupling ( κ < ),as can be shown by further bifurcation analysis. For furtherdetails on these bifurcation structures, see [7, 51].Importantly, we note that collective oscillations emergingin the Hopf bifurcation HB cease to exist in the limit of pulsesdefined by (3) with zero width obtained in the limit of s → ∞ .While this was already noted in recent studies [50, 57] webriefly outline a derivation of this fact in Appendix A. Furtherinvestigations of ours show that Hopf bifurcations continue toexist for a large range of values of the pulse shape parameter, s . Our observations suggest that the Hopf bifurcations giv-ing birth to oscillations only vanish in the limit of s → ∞ ,prompting a degeneracy for this limit. The case of infinitelynarrow pulses, s = ∞ , provides the advantage that the fixedpoint conditions resulting from the corresponding FRE can besolved in closed form, enabling a simple mathematical analy-sis. However, since this case produces a degenerate bifurca-tion behavior where limit cycles are absent, we chose to fix s = 1 .Between the pair of fold bifurcations ( SN and SN ) a pa-rameter region of bistability arises, thus facilitating hystereticbehavior. This happens for excitable neurons, ˆ η < , with excitatory coupling , κ > , as well as for parameters cor-responding to firing neurons, ˆ η > , with inhibitory cou-pling, κ < . This bistable character of solutions observedfor M = 1 population translates to the case of M = 2 popu-lations, where each population may attain distinct stable con-figurations.In the following, we fix parameter values to ˆ η = − (ex-citable neurons) and ∆ = 0 . , while varying the intra-coupling strength, κ , and the inter-coupling strength, a > . IV. ANALYSIS FOR TWO POPULATIONS
1. Symmetric and asymmetric equilibria
It is instructive to begin the analysis by surveying the pos-sible asymptotic dynamic behavior for the firing rates r and r (or equivalently, Z and Z ) in the FRE (19) and (20) for M = 2 populations. We may distinguish two types of asymp-totic states as t → ∞ , namely (i) symmetric states character-ized by r ( t ) = r ( t ) and v ( t ) = v ( t ) ; and (ii) asymmet-ric states characterized by r ( t ) (cid:54) = r ( t ) and v ( t ) (cid:54) = v ( t ) .Furthermore, each neuron population may be in a state of qui-escence (Q) or spiking (S), depending on whether r σ reflectslow or high firing activity, respectively. In an asymmetric limitcycle, both populations oscillate around a distinct value cor-responding to quiescence or spiking, respectively. Fig. 3 il-lustrates the possible asymptotic states that may be observed,depending on parameter values and initial conditions chosen.Solution branches reported in Fig. 2 for M = 1 popula-tion translate to symmetric states in the model with M = 2 populations. To see this, let us first consider two special pa-rameter choices: a = 0 (decoupled populations) and a = 1 (two populations effectively act like one large population). Inthese cases, the system with M = 2 populations displays thesame bifurcation behavior as M = 1 population, as shown inFig. 4a) for a = 0 . The branch with low firing rate (QQ) corre-sponds to quiescent neurons with coherent stationary phases;whereas the branch with high firing rate (SS) corresponds tospiking populations whose synchronization level and firingrate grow with increasing coupling strength κ . Just as for M = 1 population, the system exhibits bistable regions inwhich both configurations (quiescence and spiking) are possi-ble. However, note that in the case of a = 1 , both populationsmay only attain identical (symmetric) configurations of quies-cence or spiking, namely SS or QQ; in contrast, the decoupledcase with a = 0 trivially additionally allows for the two popu-
100 110 120 13000.10.20.30.4 r , t QQ r = r a)
100 110 120 13000.10.20.30.4 r , t SS r = r b)
100 110 120 13000.10.20.30.4 r , t QS / SQ r , r , c)
100 110 120 13010 -2 -1 r , t QS o / SQ o r , r , d)Figure 3. Asymptotic dynamic behaviors for the two populationmodel after transients. Low and high levels of the firing rate r ( t ) indicate quiescent (Q) or spiking (S) behavior. a) Symmetric (sta-ble) equilibrium where both populations are quiescent (QQ). b) Sym-metric (stable) equilibrium where both populations are spiking (SS)(Transients leading up to the equilibrium are oscillatory). c) Asym-metric (stable) equilibrium where one population is quiescent and theother spiking (QS or SQ) (Transients leading up to the equilibriumare oscillatory only for the spiking population). d) Asymmetric (sta-ble) limit cycle (QS o or SQ o ) where both populations oscillate withthe same period (but with different amplitudes). The coupling pa-rameter is κ = 1 . for panels a) through c) and κ = 2 . for panel d);parameters a = 0 . , ˆ η = − , ∆ = 0 . , s = 1 are fixed through-out. lations to attain distinct (asymmetric) configurations, namelySQ or QS. Importantly, symmetric states persist even when a (cid:54) = 0 or a (cid:54) = 1 since parameters are symmetric across thetwo populations. Specifically, if r is an equilibrium of the M = 1 population system, then so is ( r, r ) an equilibriumof the M = 2 population system but now with κ replaced by κ/ (1 + a ) . Asymmetric states, corresponding to QS or SQconfigurations with r (cid:54) = r , are trivially possible when thetwo populations are decoupled ( a = 0 ). However, their rangeof existence and stability off the degenerate cases a (cid:54) = 0 and a (cid:54) = 1 deserves further exploration, and we therefore considersmall perturbations a (cid:54) = 0 . Considering the case of a = 0 . in Fig. 4b) we observe that unstable asymmetric states (lightred) branch off the unstable symmetric state (gray) in twopitchfork bifurcations, denoted PF and PF (see also inset).The set of asymmetric equilibria forms a loop in ( κ, r ) -spacewith two folds, i.e., the equilibria undergo saddle-node bifur-cations in SN and SN . As a result, bistable asymmetricconfiguration states QS and SQ (red) become available. Thesestable branches co-exist with the bistable symmetric solutionbranches QQ and SS (black). However, for a ≥ , asym-metric states exist only in the range of intermediate couplingstrength κ . Note that, considering the case of decoupled popu-lations with a = 0 , symmetric and asymmetric branches mayappear like they coincide when inspecting Fig. 4, however,the two types of solution branches are not identical: Whilethe projections Z ∈ C and Z ∈ C indeed share identicalvalues for symmetric and asymmetric equilibria, this cannothold true in the full phase space for ( Z , Z ) ∈ C , where the -2 -1 SN SN SSQQ r σ κ a) -2 -1 SN SN SN SN SN PF PF SN SS QQQS(SQ)QS(SQ) r σ κ b)Figure 4. Bifurcation diagrams for M = 2 populations in κ forthe firing rate r σ reveal symmetric equilibria, r = r (black, gray)and asymmetric equilibria, r (cid:54) = r , (red, light red), emerging inbifurcations as follows. a) a = 0 : Both populations are decou-pled. Symmetric equilibria are folded in two saddle-node bifurca-tions ( SN , SN ) where the lower and upper branches (black) corre-sponding to (Q)uiescence and (S)piking, respectively, are stable; themiddle branch (gray) is unstable. Thus, four states are possible, SSQQ QS and SQ, facilitating multistability and hysteretic behavior. b) a = 0 . : Symmetric equilibria with r = r seen for a = 0 are stillpresent (black, gray). However, the unstable branch (gray) undergoesa pitchfork bifurcation in PF and PF , giving rise to unstable asym-metric equilibria QS / SQ with r (cid:54) = r (light red). These equilibriaare connected in a loop, folded twice in two saddle-node bifurcations( SN and SN ), giving rise to the co-existence of two stable asym-metric states (red) where one population is quiescent while the otheris spiking, (SQ or QS). Parameters are ∆ = 0 . , ˆ η = − , s = 1 ev-erywhere. Values of a chosen in this Figure are indicated as dashedlines in Fig. 6. definitions for symmetric ( r = r , v = v ) and asymmetricstates ( r (cid:54) = r , v (cid:54) = v ) are obeyed.
2. Birth and destruction of limit cycle oscillations
For larger values of the inter-coupling strength, a , asym-metric equilibria QS (SQ) may undergo Hopf bifurcationsgiving rise to limit cycle oscillations (QS o , SQ o ), indicatedby their minima/maxima (blue) in Fig. 5a) through d). Sincethese limit cycles branch off asymmetric equilibria (red), theycorrespond to asymmetric configurations characterized by fir-ing rates r ( t ) (cid:54) = r ( t ) . These limit cycles are created anddestroyed in various bifurcations, as outlined in the following. a. Birth of stable limit cycles ( HB − ). Stable limit cy-cles (blue minima/maxima) are born in the supercritical Hopfbifurcation denoted by HB − , as shown in Fig. 5a) for a =0 . . As κ increases, the amplitude waxes and wanes, as thebifurcation HB − is intersected twice in the direction of vary-ing κ , see also Fig. 6). b. Birth of stable/unstable limit cycles and annihilation insaddle-node-of-limit-cycle bifurcation ( HB − , HB + , SNLC ). Stable limit cycles (blue) are still born in a supercritical Hopfbifurcation at HB − , but now an unstable limit cycle (lightblue) of smaller amplitude emerges for greater κ in the sub-critical Hopf at HB + . The continuum of cycles folds over ina saddle-node of limit cycles bifurcation at SNLC , where thestable and unstable limit cycles coalesce and disappear, seeFig. 5b) for a = 0 . . c. Stabilization of unstable limit cycle in secondarysaddle-node-of-limit-cycle bifurcation ( SNLC ). Stable andunstable limit cycles are created in HB − and HB + . While thestable limit cycle is destroyed in the homoclinic bifurcation HC − , the unstable limit cycle is subject to a more compli-cated series of bifurcations: It undergoes not only one, buttwo saddle-node of cycles bifurcations, SNLC and SNLC .The unstable limit cycle emerging from SNLC collides withthe saddle equilibrium of the asymmetric branch in the homo-clinic bifurcation HC + and is destroyed, as shown in Fig. 5c)for a = 0 . . d. Simple birth and destruction of stable/unstable limitcycles ( HB − , HC − , HB + , HC + ). The stable and unstablelimit cycles are born in the Hopf bifurcations HB − and HB + and are destroyed in the homoclinic bifurcations HC − and HC + , respectively. The complicated scenario including twosaddle-node of limit cycles bifurcations from c) is entirely ab-sent. This simple scenario is shown in Fig. 5d) for a = 0 . .
3. Stability diagram
We now explain how the various bifurcation scenarios arerelated, i.e., how stability boundaries are connected in the ( κ, a ) -parameter plane and how bifurcation curves are struc-tured around bifurcation points of higher co-dimension.Let us first consider the overall bifurcation structure for alarger parameter range ( κ, a ) as displayed in Fig. 6a), mainlyfocusing on symmetric (QQ, SS) and asymmetric equilibria -2 -1 PF PF HB - HB - HB - HB - SN SN SN SN HB - SN m i n t ( r σ ) , m a x t ( r σ ) κ a) -2 -1 PF PF HB - HB - HB + HB + SN SN SN SN SNLC SNLC SNLC SNLC m i n t ( r σ ) , m a x t ( r σ ) κ b) SN -2 -1 PF PF HB - HB - HB + HB + SN SN SN HC - HC - HC - HC - HC + HC + HC + HC + SNLC SNLC SNLC SNLC SNLC SNLC SNLC SNLC HC + SNLC SNLC m i n t ( r σ ) , m a x t ( r σ ) κ c) -2 -1 PF PF HB - HB - HB + HB + SN SN SN SN HC - HC - HC - HC - HC + HC + HC + HC + m i n t ( r σ ) , m a x t ( r σ ) κ d)Figure 5. Birth and destruction of asymmetric limit cycle oscillations (blue) in the firing rate for M = 2 populations while varying κ . Possiblebifurcation scenarios are illustrated for various values of a (also indicated as dashed lines in Fig. 6). a) a = 0 . : stable limit cycles areborn in the supercritical Hopf bifurcation HB − (the bifurcation curve is intersected twice, see Fig. 6). b) a = 0 . : stable and unstable limitcycles are born in HB − and HB + , respectively, and annihilate in a saddle-node of limit cycles bifurcation ( SNLC ). c) a = 0 . : stableand unstable limit cycles are destroyed in a homoclinic bifurcation HC − and HC + , respectively. Moreover, the unstable cycle is subject totwo saddle-node-of-limit-cycles bifurcations at SNLC and SNLC , in between which it gains stability (see inset); the unstable limit cyclebranch emerging from SNLC gets destroyed in HC + . d) a = 0 . : the unstable limit cycle born in HB + is destroyed in the homoclinic HC + ; saddle-node-of-limit-cycles bifurcations are now absent. Symmetric equilibria connecting to PF and PF are omitted for simplicity.Stable and unstable solution branches are shown as dark and light colored shades, respectively. Parameters are ∆ = 0 . , ˆ η = − , s = 1 everywhere. QS (or SQ). On the branches of symmetric equilibria (QQ,SS)two saddle-node bifurcations occur, SN and SN (black).These bifurcations delineate the gray shaded region of bista-bility between QQ and SS. On the left of SN , stable symmet-ric equilibria correspond to low firing activity, i.e. quiescentneurons (QQ); to the right of SN , stable symmetric equi-libria correspond to high firing activity, i.e. spiking neurons (SS). Unstable saddle branches of the symmetric equilibriabetween SN and SN (gray shaded region) undergo furtherpitchfork bifurcations PF and PF (not shown for simplic-ity, see Fig. 4b) and Fig. 5 a) through d)), which give rise tounstable asymmetric unstable branches (light red). These un-stable branches gain stability on the saddle-node bifurcationcurves SN and SN (red curves in Fig. 6a)), which meet inthe codimension 2 cusp point CP . The resulting asymmetricstable configurations (QS or SQ) reside inside the red shadedregion bounded by the saddle-node curves SN and SN .Stable asymmetric equilibria QS and SQ lose stability inthe supercritical Hopf bifurcations HB − or HB’ , thus result-ing in stable asymmetric limit cycles QS o (or SQ o ) within theblue shaded regions; these limit cycles may get destroyed inthe homoclinic bifurcations denoted by HC − and HC’ (violet).Hopf ( HB − and HB’ ) and homoclinic bifurcation curves asso-ciated with the emergence and destruction of these limit cy-cles ( HC − and HC’ ) meet with the (asymmetric) saddle-nodebifurcation curve SN in two other bifurcation points of codi-mension 2, namely the Bogdanov-Takens points BT and BT’ ,respectively, characterized by double zero eigenvalues [58].The bifurcations pertaining to the asymmetric limit cyclesare structured around further, more complicated bifurcationcurves and bifurcation points of higher co-dimension, seeFig. 6b) and c). Following the Hopf bifurcation curve HB − in panel b), we arrive at a Generalized Hopf bifurcation point( GH ) of codimension 2 [58, 59]. Such a point not only has apair of purely imaginary eigenvalues, but also the first Lya-punov coefficient for the Hopf bifurcation changes sign atthis point so that subcritical ( HB + ) and supercritical ( HB − )Hopf bifurcations are separated in GH ; in addition, a branchof saddle-node of limit cycle bifurcations, SNLC , emergesfrom GH where the stable and unstable limit cycles born in HB − and HB + are annihilated.Following the bifurcation curve HB + the associated sub-critical Hopf bifurcation tangentially intersects the saddle-node bifurcation SN in the Zero-Hopf bifurcation ZH (orsaddle-node Hopf bifurcation) [58, 60], characterized by azero eigenvalue and a pair of purely imaginary eigenvalues.At ZH , the first Lyapunov coefficient vanishes once more andchanges sign. Hopf bifurcations HB + − above the ZH pointare supercritical (i.e., having a negative first Lyapunov coef-ficient), but continue to produce unstable limit cycles as thecenter manifold (of the Hopf bifurcation) is repelling.Following the saddle-node bifurcation of limit cycles curve SNLC , we observe that it terminates in another bifurca-tion point of codimension 2, Cusp of Cycles ( CPC ), whereit collides with a second saddle-node bifurcation of limitcycles curve
SNLC . This latter bifurcation curve mergeswith the homoclinic bifurcation curve HC + in a codimension ≥ point SLH , see Fig. 6c). The point
SLH separates twobranches of the homoclinic curves, HC − and HC + , and tan-gentially intersects with SNLC . Homoclinic bifurcations on HC − ( HC + ) destroy stable (unstable) limit cycles as κ ap-proaches the homoclinic bifurcation point from above (be-low).We come to the following conclusion. In similarity withthe case of M = 1 population, the cases of very small andvery strong coupling κ result in regimes with quiescent andspiking activity, respectively; both are characterized by highlevels of synchrony. In the intermediate regime, the dynamicbehavior is more complicated. We thus find the following fivestability regions: (i) for small coupling strength κ , both pop-ulations are quiescent, corresponding to the symmetric con-figuration QQ (white region); (ii) for large coupling strength BT' BT ZHGHCPCSLHCPSN SN SN SN HC'HC + HB_ HB'SNLC S p i k i ng SS QQ o r SS Q S o r S Q Q S o o r S Q o Q u i e sc en c e QQ a κ a) BT GHCPC1.6 1.8 2 2.2 2.4 2.6 2.80.20.250.30.35
SLH ZHHB_HC_ SN SN HB +– SNLCHC + SNLC HB + QQ or SS Q S o o r S Q o QS or SQ a κ b) SNLC CPC HC + QQ or SS Q S o o r S Q o SLH a κ c)Figure 6. Stability diagram for M = 2 populations in the ( K, a ) -parameter plane with ∆ = 0 . , ˆ η = − , s = 1 fixed. Gray shadingindicates a bistable region of the two stable symmetric equilibria,QQand SS. Red shading indicates the bistable region of stable asym-metric equilibria, that is, QS, or, equivalently, SQ. Blue shading indi-cates the bistable region of stable asymmetric limit cycles (QS o andSQ o ). Pitchfork bifurcations PF and PF giving rise to asymmetricunstable branches are omitted for simplicity. Dashed curves delin-eate the choices of parameter a for the bifurcation diagrams in Fig. 4and Fig. 5. κ , both populations are spiking, corresponding to the sym-metric configuration SS (white region); (iii) for intermediatecoupling strengths, we find a region of bistability between theconfigurations QQ and SS; this region of bistability co-existswith asymmetric configurations of either (iv) stationary firingrate, SQ or QS (red region), or (v) oscillatory firing rates, SQ o or QS o (blue region). Note that the permutation symmetrybetween populations 1 ↔ V. DISCUSSION
Collective oscillations in neural ensembles are responsiblefor the rhythm generation required for solving functionallyrelevant tasks in the brain [10, 18, 47]. Collective oscillationsmay be facilitated by a variety of network setups, includingheterogeneous networks with excitatory and inhibitory cou-pling leading to gamma rhythms [49]. Here, we investigatedthe emergence of collective oscillations in a simple modelconsisting of a homogeneous network composed of two (sta-tistically) identical populations of type 1 neurons with non-uniform but symmetric coupling, i.e., neurons are coupledwith strength κ and aκ (with a (cid:54) = 1 ) within and between thetwo populations, respectively.In this model, each population may assume states corre-sponding to quiescent (Q) or spiking (S) firing activity. Thus,we may distinguish symmetric configurations, where bothpopulations are either quiescent or spiking (QQ, SS), andasymmetric configurations, where one population is quies-cent but the other is spiking (SQ, QS). We found that stablesymmetric configurations may coexist for certain parameterchoices (see Fig. 4a)). We did not find that symmetric config-urations are oscillatory except for uniform coupling ( a = 1 ) orfor absent inter-coupling ( a = 0 ). As we deviate from uniformcoupling, a (cid:54) = 1 , unstable asymmetric equilibria emerge fromsymmetric configurations in symmetry-breaking pitchfork bi-furcations. Along these solution branches, asymmetric equi-libria may further undergo saddle-node bifurcations and thusgain stability (see Fig. 4b)). Asymmetric oscillatory configu-rations (QS o , SQ o ) emerge in Hopf bifurcations (Fig. 5) thatare organized around higher codimension bifurcation points.Depending on parameters, symmetric and asymmetric con-figurations may be stable and co-exist, resulting in multista-bility between either stationary configurations only (QQ, SS,QS, SQ); or between stationary and oscillatory configurations(QQ, SS, QS o , SQ o ). For these regions of stability we havedetermined valid parameter regions and stability boundaries(Fig. 6).Oscillator networks with such modular network structureare known to exhibit a high degree of multistability, i.e., de-pending on initial conditions, a variety of dynamic config-urations for the collective states may be assumed in eachpopulation. A prominent example are synchronization pat-terns known as chimera states in Kuramoto oscillator net-works [27, 29, 36], which may be employed to store memoryor perform computations [61] or direct the flow of informa-tion between populations [8, 26]. However, compared to Ku- ramoto networks with rigidly rotating oscillators, the excitablenature of neurons intrinsically leads to more complicated dy-namics and synchronization behavior [62]. While compli-cated dynamics may arise in networks composed of identicalKuramoto oscillators arranged with at least two populations(as well as for broken parameter symmetries) [28, 63], ex-citable type 1 neurons produce rich bifurcation behavior andbistability already for a single population, as illustrated inFig. 4 and discussed in [7, 8]. Such multistability is of greatinterest in applications, e.g. in neuroscience. A recent studymodeled networks of type 1 neurons and demonstrated howthe bistability between low and high firing activity — result-ing in a large configuration space that scales with the numberof populations — may be employed to solve cognitive taskssuch as memory storage and recall [39].Several studies considered networks of type 1 neurons interms of their macroscopic behavior. The collective dy-namics of a single population was studied in terms of non-identical Theta neurons with non-zero pulse width [7, 51], ofthe response to an external (rigid) forcing [38], of quadraticintegrate-and-fire neurons [6], and of different coupling func-tions, oscillations and aging transitions [57], and of the roleof distributed delay in the coupling function [64]. Luke etal [38] studied a two population model similar to ours; how-ever, they considered unidirectional coupling. This driver-response system exhibits some of the bifurcation structuresand collective macroscopic behaviors that we reported here,i.e., the response population exhibits equilibrium states, limitcycle states and multistable behavior. However, unlike theirstudy, we did not allow for asymmetric coupling and did notobserve quasiperiodic and chaotic dynamics. Future researchmight address the question if breaking parameter symmetriesbetween the populations (such as the coupling strength) mayinduce bifurcations leading to chaos. For such cases, onemay for instance envision torus bifurcations emerging fromthe Zero-Hopf bifurcations [58], offering a route to chaosvia bifurcations of Shil’nikov homoclinic orbits to saddlefoci. Ratas and Pyragas [65] studied a network of quadraticintegrate-and-fire neurons with two populations. While theirsystem is similar to ours, it differs in some important aspects.Firstly, neurons are considered to be strongly heterogeneouswith an excitability spread around ˆ η = 0 , thus resulting ina network including both excitable and spiking neurons; weconsider excitable neurons only. Secondly, for the couplingfunction, they use a threshold modulation coupling functioncorresponding to a Heaviside function. This system exhibitssteady and oscillatory states with symmetric and asymmetriccharacter, but unlike our system, also chaotic behavior andstates characterized by anti-phase configurations.To study collective oscillations of firing activity in ourmodel, it is necessary to deviate from the case of instantaneouspulse coupling ( s → ∞ ) where collective oscillations are ab-sent (see Appendix A and [50, 57]). The pulse width given byEq. (3) with s = 1 was large; other pulse shape models [4, 6]may be more realistic and consider that incident pulses arriveinstantaneously in order to decay exponentially fast upon ar-rival over a characteristic time scale τ . It is then frequentlyassumed that τ → , resulting in time-symmetric and instan-0taneous pulses. This strategy certainly simplifies analysis; yet,it appears that this limit biophysically is no more realistic, es-pecially since it results in the same macroscopic equations asgiven by (19) and (20) for the limit of s → ∞ ; this again rulesout the potential to produce any macroscopic oscillations. Fora future study it might be interesting to examine how the spe-cific choice of pulse shape in terms of width and asymmetryaffects the unfolding of bifurcations.Many questions remain. For instance, breaking parame-ter symmetry may result in richer dynamics [30] includingchaos [63]; is chaotic motion feasible if excitability param-eters ( ˆ η σ and ∆ σ ) are non-identical, or if small delay is intro-duced in the coupling? In terms of switching between config-urations and devising a control method to do this, it may beuseful to determine basins of attraction for the various config-urations or responses to directed perturbations [28, 39]. Fur-thermore, networks with larger population number M > provide a larger set of dynamic configurations [39, 66, 67];but how large is the set of configurations as a function of thepopulation number, and which of the configurations are sta-ble, and which oscillatory? Future studies may address suchand further questions. ACKNOWLEDGEMENTS
The authors would like to thank C. Bick and B. Pietras forhelpful discussions, and A. Torcini and P. So for helpful cor-respondence.
Appendix A: Collective oscillations for non-zero pulse width
We briefly discuss the existence of Hopf bifurcations andresulting limit cycle oscillations in the firing rate r ( t ) for vary-ing pulse shape parameter s , for the simple case of M = 1 population. To determine the presence of Hopf instabilities,we examine eigenvalues of the Jacobian of (19), J = (cid:18) v r − π r + ∂∂r I v + ∂∂v I (cid:19) . (A1)Steady state of (19) implies v ∗ = − ∆2 π r ∗ so thattr ( J ) = − π r ∗ + ∂∂v I | ( r,v )=( r ∗ ,v ∗ ) . (A2)A necessary condition for a Hopf bifurcation is that tr ( J ) =0 . For the case of infinitely narrow pulses, s = ∞ , Hopfbifurcations are impossible: we have ∂∂v I = κ ∂∂v P ( ∞ ) = 0 thus and tr ( J ) = − /π/r ∗ < for all r ∗ > . Hence, Hopfbifurcations and resulting limit cycles regardless of the choiceof parameters can be ruled out for this case.Conversely, we know that Hopf bifurcations are possiblefor s = 1 (see Fig. 2) and s = 2 (see Luke et al. [7]). Indeed,the trace for < s < ∞ involves more complicated terms,and Hopf bifurcation cannot easily be ruled out. While ananalytical proof remains elusive, using a numerical analysisbased on solving the zero trace condition, one finds that limitcycle oscillations are feasible for a large range of pulse shapeparameters, including at least ≤ s ≤ . [1] G. B. Ermentrout and N. Kopell. Parabolic Bursting in an ex-citable system coupled with a slow oscillation. SIAM Journalon Applied Mathematics , 46(2):233–253, 1986.[2] P. E. Latham, B. J. Richmond, P. G. Nelson, and S. Nirenberg.Intrinsic dynamics in neuronal networks. I. Theory.
Journal ofNeurophysiology , 83(2):808–827, 2000.[3] D. Hansel and G. Mato. Existence and stability of persistentstates in large neuronal networks.
Physical Review Letters ,86(18):4175–4178, 2001.[4] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski.
Neu-ronal dynamics: From single neurons to networks and modelsof cognition . Cambridge University Press, 2014.[5] E. Ott and T. M. Antonsen. Low dimensional behavior oflarge systems of globally coupled oscillators.
Chaos (Wood-bury, N.Y.) , 18(3):037113, sep 2008.[6] E. Montbri´o, D. Paz´o, and A. Roxin. Macroscopic descriptionfor networks of spiking neurons.
Physical Review X , 5(2):1–15,2015.[7] T. B. Luke, E. Barreto, and P. So. Complete Classification of theMacroscopic Behavior of a Heterogeneous Network of ThetaNeurons.
Neural Computation , 25:3207–3234, 2013.[8] C. Bick, C. Laing, M. Goodfellow, and E.A. Martens. Under-standing the dynamics of biological and neural oscillator net-works through exact mean-field reductions: a review.
Journalof Mathematical Neuroscience , 9(10), 2020.[9] A. Byrne, D. Avitabile, and S. Coombes. Next-generation neu-ral field model: The evolution of synchrony within patterns and waves.
Physical Review E , 2019.[10] E. Marder and D. Bucher. Central pattern generators and thecontrol of rythmic movements.
Current Biology , 11:R986–R996, 2001.[11] J. C. Smith, H. H. Ellenberger, K. Ballanyi, D. W. Richter, andJ. L. Feldman. Pre-B¨otzinger Complex: A Brainstem RegionThat May Generate Respiratory Rhythm in Mammals.
Science ,254(5032):726–729, 1991.[12] G. Buzs´aki and X.-J. Wang. Mechanisms of Gamma Oscilla-tions.
Annual Review of Neuroscience , 35(1):203–225, 2012.[13] E. T. Bullmore and O. Sporns. Complex brain networks: graphtheoretical analysis of structural and functional systems.
Naturereviews. Neuroscience , 10(3):186–98, 2009.[14] D. Meunier, R. Lambiotte, and E. T. Bullmore. Modular and hi-erarchically modular organization of brain networks.
Frontiersin Neuroscience , 4(DEC):1–11, 2010.[15] K. D. Harris. Neural signatures of cell assembly organization: Article : Nature Reviews Neuroscience.
Nature Reviews ,6(May):399–407, 2005.[16] C. W. Lynn and D. S. Bassett. The physics of brain networkstructure, function and control.
Nature Reviews Physics , 2019.[17] L. Glass. Synchronization and rhythmic processes in physiol-ogy.
Nature , 410(March):277–284, 2001.[18] P. J. Uhlhaas and W. Singer. Neural synchrony in brain disor-ders: relevance for cognitive dysfunctions and pathophysiology.
Neuron , 52(1):155–68, 2006.[19] J. Fell and N. Axmacher. The role of phase synchronization in memory processes. Nature Reviews Neuroscience , 12(2):105–118, 2011.[20] P. Fries. Neuronal gamma-band synchronization as a funda-mental process in cortical computation.
Annual Review of Neu-roscience , 32:209–24, 2009.[21] X. J. Wang. Neurophysiological and Computational Princi-ples of Cortical Rhythms in Cognition.
Physiological Reviews ,90(3):1195–1268, 2010.[22] W. Singer and C. M. Gray. Visual feature integration and thetemporal correlation hypothesis.
Ann Rev Neurosci , 18:555–586, 1995.[23] P. Fries. A mechanism for cognitive dynamics: neuronal com-munication through neuronal coherence.
Trends in CognitiveSciences , 9(10):474–80, 2005.[24] M. I. Rabinovich, V. S. Afraimovich, C. Bick, and P. Varona. In-formation flow dynamics in the brain.
Physics of Life Reviews ,9(1):51–73, 2012.[25] C. Kirst, M. Timme, and D. Battaglia. Dynamic informa-tion routing in complex networks.
Nature Communications ,7:11061, 2016.[26] N. Deschle, A. Daffertshofer, D. Battaglia, and E.A. Martens.Directed Flow of Information in Chimera States.
Frontiers inApplied Mathematics and Statistics , 5, 2019.[27] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley.Solvable model for chimera states of coupled oscillators.
Phys.Rev. Lett , 101:084103, 2008.[28] E. A. Martens, M. J. Panaggio, and D. M. Abrams. Basins ofAttraction for Chimera States.
New Journal of Physics, FastTrack Communication , 18:022002, 2016.[29] E. A. Martens. Bistable Chimera Attractors on a Triangu-lar Network of Oscillator Populations.
Physical Review E ,82(1):016216, jul 2010.[30] E. A. Martens, C. Bick, and M. J. Panaggio. Chimera statesin two populations with heterogeneous phase-lag.
Chaos ,26(9):094819, 2016.[31] E. A. Martens, S. Thutupalli, A. Fourri`ere, and O. Hallatschek.Chimera States in Mechanical Oscillator Networks.
Proc. Natl.Acad. Sci. , 110(26):10563–10567, 2013.[32] S. Olmi, E. A. Martens, S. Thutupalli, and A. Torcini. Intermit-tent chaotic chimeras for coupled rotators.
Phys. Rev. E RapidCommunications , 92:030901, 2015.[33] C. R. Laing. The dynamics of chimera states in heteroge-neous Kuramoto networks.
Physica D: Nonlinear Phenomena ,238(16):1569–1588, aug 2009.[34] C. R. Laing, K. Rajendran, and I. G. Kevrekidis. Chimeras inrandom non-complete networks of phase oscillators.
Chaos: AnInterdisciplinary Journal of Nonlinear Science , 22(1):013132,2012.[35] E. Sch¨oll. Synchronization patterns and chimera states in com-plex networks: Interplay of topology and dynamics.
The Euro-pean Physical Journal Special Topics , 225(6-7):891–919, 2016.[36] M. J. Panaggio and D. M. Abrams. Chimera states: coexistenceof coherence and incoherence in networks of coupled oscilla-tors.
Nonlinearity , 28(3):R67, 2015.[37] C. R. Laing. Phase oscillator network models of brain dy-namics. In Ahmed A. Moustafa, editor,
Computational Mod-els of Brain and Behavior , chapter 37, pages 505–518. Wiley-Blackwell, 2017.[38] T. B. Luke, E. Barreto, and P. So. Macroscopic complexity froman autonomous network of networks of theta neurons.
Frontiersin Computational Neuroscience , 8(NOV):1–11, 2014.[39] H. Schmidt, D. Avitabile, E. Montbri´o, and A. Roxin. Net-work mechanisms underlying the role of oscillations in cogni-tive tasks.
PLoS Computational Biology , 14(9):1–24, 2018. [40] M. Yamakou, P.G. Hjorth, and E. A. Martens. Optimal self-induced stochastic resonance in multiplex neural networks:electrical versus chemical synapses.
Frontiers in Neuroscience,arXiv:2002.12443 , 2020 (provisionally accepted.[41] G. Weerasinghe, B. Duchet, H. Cagnan, P.r Brown, C. Bick,and R. Bogacz. Predicting the effects of deep brain stimula-tion using a reduced coupled oscillator model. bioRxiv:448290 ,2018.[42] S. Watanabe and S. H. Strogatz. Constants of motion for su-perconducting Josephson arrays.
Physica D , 74(3-4):197–253,1994.[43] C. R. Laing. The Dynamics of Networks of Identical ThetaNeurons.
Journal of Mathematical Neuroscience , 8(1), 2018.[44] H. R. Wilson and J. D. Cowan. Excitatory and inhibitory inter-actions in localized populations of model neurons.
Biophysicaljournal , 12(1):1–24, 1972.[45] S. Amari. Dynamics of pattern formation in lateral-inhibitiontype neural fields.
Biological Cybernetics , 27(2):77–87, 1977.[46] L. Lin, E. Barreto, and P. So. Synaptic Diversity SuppressesComplex Collective Behavior in Networks of Theta Neu-rons.
Frontiers in Computational Neuroscience , 14(May):1–17,2020.[47] G. Buzs´aki and B. O. Watson. Brain rhythms and neural syn-tax: Implications for efficient coding of cognitive content andneuropsychiatric disease.
Dialogues in Clinical Neuroscience ,2012.[48] P. Fries, D. Nikoli´c, and W. Singer. The gamma cycle.
Trendsin neurosciences , 30(7):309–316, 2007.[49] M. Segneri, H. Bi, S. Olmi, and A. Torcini. Theta-nestedgamma oscillations in next generation neural mass models.
Frontiers in Computational Neuroscience , pages 1–28, 2020.[50] F. Devalle, A. Roxin, and E. Montbri´o. Firing rate equationsrequire a spike synchrony mechanism to correctly describe fastoscillations in inhibitory networks.
PLoS Computational Biol-ogy , 13(12):1–21, 2017.[51] P. So, T. B. Luke, and E. Barreto. Networks of theta neuronswith time-varying excitability: Macroscopic chaos, multistabil-ity, and final-state uncertainty.
Physica D: Nonlinear Phenom-ena , 267:16–26, may 2013.[52] A. Dhooge, W. Govaerts, Yu A. Kuznetsov, H. G.E. Meijer, andB. Sautois. New features of the software MatCont for bifurca-tion analysis of dynamical systems.
Mathematical and Com-puter Modelling of Dynamical Systems , 2008.[53] B. Ermentrout. Ermentrout-Kopell canonical model.
Scholar-pedia , 3(3):1398, 2008.[54] J. T. Ariaratnam and S. H. Strogatz. Phase diagram for the Win-free model of coupled nonlinear oscillators.
Physical ReviewLetters , 86(19):4278, 2001.[55] G. B. Ermentrout and D. H. Terman.
Mathematical Founda-tions of Neuroscience , volume 35 of
Interdisciplinary AppliedMathematics . Springer, New York, NY, 2010.[56] N. Brunel and P. E. Latham. Firing rate of the noisy quadraticintegrate-and-fire neuron.
Neural Computation , 15(10):2281–2306, 2003.[57] I. Ratas and K. Pyragas. Macroscopic self-oscillations and ag-ing transition in a network of synaptically coupled quadraticintegrate-and-fire neurons.
Physical Review E , 94(3):1–11,2016.[58] Y. A. Kuznetsov.
Elements of applied bifurcation theory .Springer-Verlag, New York, 1998.[59] J. Guckenheimer and Y. Kuznetsov. Bautin bifurcation.
Schol-arpedia , 2007.[60] J. Guckenheimer and Y. Kuznetsov. Fold-Hopf bifurcation.
Scholarpedia , 2007. [61] C. Bick and E. A. Martens. Controlling Chimeras. New Journalof Physics , 17:033030, 2015.[62] D. C˘alug˘aru, J. F. Totz, E. A. Martens, and H. Engel. First-ordersynchronization transition in a large population of relaxationoscillators.
Science Advances, arxiv:1812.04727 , 6, 2020.[63] C. Bick, M. J. Panaggio, and E. A. Martens. Chaos in KuramotoOscillator Networks.
Chaos , 28:071102, 2018.[64] I. Ratas and K. Pyragas. Macroscopic oscillations of a quadraticintegrate-and-fire neuron network with global distributed-delaycoupling.
Physical Review E , 98(5):1–11, 2018. [65] I. Ratas and K. Pyragas. Symmetry breaking in two interactingpopulations of quadratic integrate-and-fire neurons.
PhysicalReview E , 96(4):1–9, 2017.[66] M. Shanahan. Metastable chimera states in community-structured oscillator networks.
Chaos (Woodbury, N.Y.) ,20(1):013108, mar 2010.[67] M. Wildie and M. Shanahan. Metastability and chimerastates in modular delay and pulse-coupled oscillator networks.