From quasiperiodic partial synchronization to collective chaos in populations of inhibitory neurons with delay
aa r X i v : . [ n li n . C D ] J un From quasiperiodic partial synchronization to collective chaos in populations of inhibitory neuronswith delay
Diego Paz´o and Ernest Montbri´o Instituto de F´ısica de Cantabria (IFCA), CSIC-Universidad de Cantabria, 39005 Santander, Spain Center for Brain and Cognition, Department of Information and CommunicationTechnologies, Universitat Pompeu Fabra, 08018 Barcelona, Spain (Dated: October 24, 2018)Collective chaos is shown to emerge, via a period-doubling cascade, from quasiperiodic partial synchroniza-tion in a population of identical inhibitory neurons with delayed global coupling. This system is thoroughlyinvestigated by means of an exact model of the macroscopic dynamics, valid in the thermodynamic limit. Thecollective chaotic state is reproduced numerically with a finite population, and persists in the presence of weakheterogeneities. Finally, the relationship of the model’s dynamics with fast neuronal oscillations is discussed.
PACS numbers: 05.45.Xt 87.19.lm 02.30.Ks
Electrical measurements of brain activity display a broadspectrum of oscillations, reflecting the complex coordinationof spike discharges across large neuronal populations [1]. Aparticularly fruitful theoretical framework for investigatingneuronal rhythms is to model networks of neurons as pop-ulations of heterogeneous oscillators [2–4]. These modelsexhibit a prevalent transition from incoherence to partial co-herence, when a fraction of the oscillators becomes entrainedto a common frequency. As a result a macroscopic oscillatorymode appears with the same frequency as that of the synchro-nized cluster [2, 5].Yet, even populations of globally coupled identical oscil-lators are capable of exhibiting a much wider diversity ofcomplex oscillatory states, see [6] for a recent survey. Ingeneral, this is due to the complexity of the coupling func-tions and of the individual oscillators. A relevant exampleis the so-called quasiperiodic partial synchronization (QPS),which has been extensively investigated in networks of exci-tatory leaky integrate-and-fire (LIF) neurons [7–11], as wellas in populations of limit-cycle oscillators and phase oscilla-tors [12–18]. In QPS, the network sets into a nontrivial dy-namical regime in which oscillators display quasiperiodic dy-namics while the collective observables oscillate periodically.Remarkably, the period of these oscillations differs from themean period of the individual oscillators. As pointed out re-cently [17], this interesting property of QPS is shared by the collective chaos observed in populations of globally coupledlimit-cycle oscillators [19–25]. Here, the collective chaoticmode is typically accompanied by microscopic chaotic dy-namics at the level of the individual oscillators. However, asnoticed in [20], populations of limit-cycle oscillators may alsodisplay pure collective chaos without trace of orbital instabil-ity at the microscopic level. In this state the coordinates of theoscillators fall on a smooth closed curve and no mixing oc-curs, what points to the existence of collective chaos in popu-lations of oscillators governed by a single phase-like variable.In this Letter we uncover the spontaneous emergence ofpure collective chaos from QPS, via a cascade of period-doubling bifurcations. Notably, this is found in a simple pop-ulation of identical integrate-and-fire oscillators with time- delayed pulse coupling, which is thoroughly analyzed withinthe framework of the so-called Ott-Antonsen theory [26–29].Moreover, we show that pure collective chaos persists whenweak heterogeneities are considered. This suggests that cer-tain forms of irregular collective motion observed in large net-works of heterogeneous LIF neurons with delay [30] may bealready found for identical neurons.We investigate a model consisting of a population of N ≫ neurons, with membrane potentials { V j } j =1 ,...,N . The evo-lution of V j is governed by the so-called quadratic integrate-and-fire (QIF) model, which obeys the nonlinear differentialequation [31–33] τ ˙ V j = V j + I j , (1)where τ is the neuron’s membrane time constant. When V j reaches the value V p , the QIF neuron emits a spike, and V j isreset to V r . Thereafter we consider V p = − V r = ∞ [34]. Inthis case the model (1) can be exactly transformed to a phasemodel called theta-neuron [31, 33]. The external inputs I j have the form I j = η j + J s D , (2)where parameters η j determine the dynamics of each uncou-pled neuron, J = 0 : Those neurons with η j < are excitable,whereas neurons with η j > behave as self-sustained oscil-lators with period, or interspike interval ISI j = πτ / √ η j . InEq. (2), the delayed mean activity s D ≡ s ( t − D ) is definedsumming the spikes of all neurons: s D = τN τ s N X j =1 X k Z t − Dt − D − τ s δ ( t ′ − t kj ) dt ′ . (3)In this equation, t kj is the time of the k th spike of j th neuron,and δ ( t ) is the Dirac delta function. We assume the thermo-dynamic limit N → ∞ , so that a second limit in the tem-poral window τ s → leads to the relationship s D = τ r D ,where r D ≡ r ( t − D ) is the time-delayed firing rate, i.e. thepopulation-averaged number of spikes per unit time. Thestrength of the interactions is controlled in Eq. (2) by the n . i nde x I S I ( k + ) I S I ( k + ) n . i nde x
20 40 60 80 t n . i nde x ISI(k) I S I ( k + ) (a)(c)(e) (f)(d)(b) FIG. 1. Quasiperiodic partial synchronization (a-d) and collectivechaos (e,f) in an inhibitory network of N = 1000 identical QIF neu-rons with η j = ¯ η = τ = 1 , τ s = 10 − , and: (a,b) D = 2 . , J = − . , (c,d) D = 2 . , J = − . , and (e,f) D = 3 , J = − . . (a,c,e) Raster plots of 200 randomly selected neurons.Dots corresponds to a firing events, and neurons are indexed accord-ing to their firing order. (b,d,f) Return plots for interspike inter-vals ISI j ( k ) = t k +1 j − t kj of an arbitrary neuron j . synaptic weight constant J , which can be either positive ornegative for excitatory or inhibitory synapses, respectively.We start performing numerical simulations of an inhibitory( J < ) population of identical neurons with η j = ¯ η > . Inpanels (a) and (c) of Fig. 1, showing raster plots for two valuesof J , the system exhibits QPS. In fact, the return plots in pan-els (b,d) show a closed line indicating quasiperiodic single-neuron dynamics, see [7]. Remarkably, for certain values ofthe time delay D , see Fig. 1(e,f), increasing inhibition leadsto a different macroscopic state, where neurons exhibit irreg-ular dynamics whereas the macroscopic dynamics is chaotic,as shown below.Where and how QPS and collective chaos emerge is inves-tigated next. To this aim, we follow [29] and using the Ott-Antonsen theory (by means of a Lorentzian ansatz) derive theso-called firing-rate equations (FREs) governing the dynam-ics of the firing rate r , and the population’s mean membranepotential v . Considering that currents η j are distributed ac-cording to a Lorentzian distribution of half-width ∆ , centeredat ¯ η , g ( η ) = (∆ /π )[( η − ¯ η ) + ∆ ] − , we obtain a system ofone ordinary and one delay differential equations [35] τ ˙ r = ∆ πτ + 2 rv, (4a) τ ˙ v = v + ¯ η + Jτ r D − τ π r , (4b)which exactly describe the macroscopic dynamics of the sys-tem in the infinite N limit [36]. Hereafter we set τ = ¯ η = 1 in Eq. (4), without lack of generality [37].Figures 2(a,c,e) display the time series of the population fir-ing rate, using the parameters of Fig. 1, for both the networkof spiking neurons (1) and the FREs (4). The attractor of theFREs in Fig. 2(b,d,f) is in perfect agreement with the globalbehavior of the population. Figures 1(b,d) and 2(a,c) displaythe fingerprint of QPS: oscillations of the mean field, with a v r v r t r r -20020 v (a)(c)(e) (b)(d)(f) FIG. 2. Macroscopic dynamics of quasiperiodic partial synchroniza-tion (a-d), and collective chaos (e,f); same parameters as in Fig. 1are used. Black curves are obtained integrating the FREs (4), with ∆ = 0 . Red curves in (a,c,e) are the firing rates obtained fromnumerical simulation of the population, Eqs. (1)-(3). Right pan-els (b,d,f): Phase portraits of the FREs. The unstable fixed point(brown circle) and the unstable orbit (brown, dashed line) corre-spond to incoherence and full synchronization, respectively. Thethree largest Lyapunov exponents of the chaotic attractor in panel(f) are { . , , − . } . different period (in the present case, longer) than the individ-ual neurons ISIs. Noteworthy, the two oscillations shown inFigs. 2(a,c) have exactly the same period: T = 2 D . Thisis the consequence of the symmetry of the limit cycle under v → − v , see Figs. 2(b,d). Indeed, using Eqs. (4) with ∆ = 0 ,one finds that this symmetric cycle is only possible if the pe-riod of the oscillation satisfies: T m = 2 Dm , with m = 1 , , . . . (5)As parameters are varied, the reflection symmetry of the limitcycle breaks down at a period-doubling bifurcation. More-over, the inset of Fig. 3 shows that this bifurcation is fol-lowed by a period-doubling cascade as parameter J is varied,giving rise to a state of collective chaos as that of Fig. 2(f).Remarkably, though the collective dynamics is chaotic, thesingle-neuron evolution is not. Indeed, as a consequence ofthe mean-field character of the model and its first order kinet-ics, the firing order of the neurons is preserved (i.e. neuron j always fires just before neuron j − ) and mixing is not pos-sible.In the following, we analyze the FREs (4) in detail, whatpermits to elucidate why collective chaos is found only in acertain range of delays, and only for inhibitory coupling. Foridentical neurons, ∆ = 0 , the only fixed point is ( r s , v s ) = (cid:0) ( J + √ J + 4 π ) / (2 π ) , (cid:1) , corresponding to an incoher-ent state. Its stability can be determined linearizing aroundthe fixed point r ( t ) = r s + δr e λt and v ( t ) = δv e λt , and im-posing the condition of marginal stability: λ = i Ω . We find afamily of Hopf instabilities at J ( n ) H = π (Ω n − × ( (6Ω n + 12) − / for odd n (2Ω n − − / for even n (6) D π π π π J -4-2024 -6 -2-.2 0 D = 3 J L y a p . e xp . FIG. 3. Stability regions of the incoherent and fully synchronizedstates for ∆ = 0 . Shaded region: Incoherence is stable. Hatchedregion: Full synchrony is unstable . Dark gray (blue) and the lightgray (red) lines are the loci of sub- and super-critical Hopf bifurca-tions of incoherence, respectively. The approximate periodicity ofthe phase diagram with D stems from the ISI = π of the uncou-pled neurons. Inset: Three largest (collective) Lyapunov exponentsin the range − < J < − . for D = 3 , computed numericallyfrom Eq. (4) using the method in [43]. Note the supercritical Hopfbifurcation at J (1) H = − . . . . and the accumulation of period-doubling bifurcations at J ≈ . . with associated frequencies Ω n = nπ/D . The line with sev-eral cusps depicted in Fig. 3 correspond to the boundaries ofincoherence given by Eq. (6). The blue and red colors indicatethe sub- and super-critical character of the bifurcation, respec-tively, and have been calculated perturbatively [38]. The sta-bility region of incoherence (shaded) closely resembles thatof the Kuramoto model of coupled oscillators, with alternat-ing domains at positive and negative J values as time delayis increased [39–42]. However, the presence of supercriticalHopf bifurcations in some ranges of the inhibitory part of thediagram is a distinct and important feature of model (1)-(2),as we show below.We also calculated the stability boundaries of the fully syn-chronized states, V j ( t ) = v ( t ) , which are given by the familyof functions J ( m ) c = 2 cot (cid:18) Dm (cid:19) , with m = 1 , , , . . . (7)and by evenly spaced vertical lines at D = nπ , with n =1 , . . . [44]. Accordingly, the regions of unstable full syn-chrony correspond to the hatched regions of the phase diagramFig. 3. Note that for weak coupling, i.e. close to the J = 0 axis, the phase diagram in Fig. 3 is fully consistent with that ofthe Kuramoto model with delay [39], as it can be proven ap-plying the averaging approximation to model (1) with ∆ = 0 ,see [5, 28]. Specifically, we observe three qualitatively differ-ent regions at small | J | : Incoherence (shaded-hatched), oneor more fully synchronized states (white-unhatched), and co-existence between incoherence and full synchrony (shaded-unhatched). (a) ∆ = 0 π D (1) c D Full sync. Q P S (b) ∆ = 0 D FIG. 4. (a) Sketch of the transition between full synchronizationand QPS as D is changed ( J < − . , ∆ = 0 ); the y -axis is anarbitrary coordinate like e.g. the minimal value attained by the firingrate. Full synchrony undergoes two transcritical bifurcations. Thefirst one, at D (1) c , gives rise to QPS which may, eventually, undergosecondary instabilities (not depicted) leading to chaos. At D = π full synchrony recovers its stability abruptly. (b) Same sketch for ∆ = 0 . Left transcritical bifurcation evaporates while the right onebecomes a saddle-node bifurcation. Away from the weak coupling regime, the system displayscollective phenomena unseen in the Kuramoto model. Insidethe unshaded-hatched region, located below the Hopf curve J (1) H , both incoherence and synchronization are simultane-ously unstable . Moreover, due to the supercritical characterof the Hopf boundary J (1) H in the range . < D < . [45], QPS emerges as a stable, small-amplitude oscillatory so-lution —as that of Fig. 2(a)— bifurcating from incoherencewith period T = 2 D .Additionally, QPS can also emerge via the destabilizationof full synchronization at J (1) c . The simulation of the FREsconfirms the prediction of Eq. (7), and allows to complete asomewhat peculiar picture: The fully synchronous state is adegenerate, infinitely long trajectory along the v -axis, and thelimit cycle corresponding to QPS emanates from it with anunbounded size —see Fig. 2(d), for a situation not far awayfrom the bifurcation point. In Fig. 4(a), a sketch of the bi-furcation diagram (valid for J < − . and D around π ) isdepicted. Stable QPS bifurcates from the fully synchronousstate at D (1) c = arctan(2 /J ) < π , through a transcritical bi-furcation of limit cycles. Then, in a second transcritical bifur-cation at D = π (involving unstable QPS), the synchronizedstate recovers its stability. This scenario implies the existenceof a region of bistability between QPS (or collective chaos)and full synchronization for D > π —in consistence, again,with the supercritical character of the Hopf bifurcation J (1) H for D < . .So far, we have concentrated on identical QIF neurons.Our final results concern the robustness of QPS and collectivechaos against heterogeneity. In the presence of heterogeneityfull synchronization and QPS cannot be observed, but statesreminiscent of them persist, as sketched in Fig. 4(b). Indeed,as the transcritical bifurcation is fragile, the bifurcation origi-nally located at D = π is replaced by a saddle-node bifurca-tion, whereas the other bifurcation at D = D (1) c vanishes.Regarding collective chaos, Fig. 5(a-c) shows numericalsimulations of the heterogeneous QIF neurons (1), with pa-rameter values close to those of Fig. 1(e). We observe inFig. 5(c) synchronized clusters at different average ISIs. Us-ing the FREs (4) we checked that (i) the macroscopic infinite- neu r on i nde x < I S I j > t r neuron index ( j ) . . . - λ j (a) (c)(d)(b) FIG. 5. (a) Raster plot and (b) Time series of the mean firing rate ofa population of N = 1000 heterogeneous ( ∆ = 0 . ) QIF neu-rons, with inhibitory coupling ( J = − . ) and delay ( D = 3 . ).(c) Time-averaged ISIs of individual neurons vs. neuron index (redpoints)—labels are associated to values of η j . As a double-check,the FREs (4) are simulated for the same parameters and used to driveindividual neurons. We obtain in this way the h ISI j i depicted byblack circles. Note the horizontal segments corresponding to clus-ters of neurons with identical h ISI i . (d) Lyapunov exponents λ j ofindividual neurons driven by FREs (4). The λ j ’s are all negative,while λ j = 0 for the ∆ = 0 case in Fig. 1(e). N dynamics of the model is indeed chaotic with leading Lya-punov exponents { . , , − . } ; (ii) the microscopic dy-namics is stable as revealed by the Lyapunov exponents ob-tained forcing each neuron by Eq. (4), see Fig. 5(d). Interest-ingly, a similar state was numerically uncovered in [30], andits chaotic nature was attributed to the presence of quenchedheterogeneity. However our conclusion is quite the opposite:the chaotic state in Fig. 5 can be regarded as a perturbedversion of the collective chaos in Figs. 1(e) and 2(e), andtherefore heterogeneity is not essential for observing collec-tive chaos.The fact that the model studied here exhibits nontrivial dy-namics precisely for inhibitory coupling —in contrast to theprevious studies using LIF neurons [7–11]—, deserves to beemphasized. A large body of data demonstrate that brain os-cillations in the gamma and fast frequency ranges (30-200 Hz)are inextricably linked to the behavior of populations of in-hibitory neurons [1, 46–48]. Moreover, theoretical and com-putational studies indicate that these oscillations emerge asa consequence of the interplay between inhibition and thesignificant time delays produced by synaptic processing, seee.g. [46, 49]. Our results add to these body of work, showingthat QPS and collective chaos also arise in simple inhibitorypopulations of phase oscillators with delayed pulse coupling.Plausible values for the synaptic delays are of the order of D ∼ ms, so that the QPS state studied here necessarily hasa frequency f ∼ (2 D ) − = 100 Hz, corresponding to fastbrain oscillations. This is in agreement with the frequency ofthe oscillations displayed by heuristic firing rate models withfixed time delays and inhibitory coupling [49–52]. Exactlythe same range of frequencies is also observed in networks of identical, noise-driven inhibitory neurons with synaptic de-lays, in the so-called sparsely synchronized state [49, 53–55].Remarkably, sparse synchronization also displays a macro-scopic/microscopic dichotomy, similar to that of the QPS andcollective-chaos states analyzed here.The analysis of the thermodynamic limit of the model (1)-(2) by means of the firing-rate equations (4), permits to dissectmacroscopic from microscopic dynamics in that limit. Thisstrategy seems to be particularly useful for investigating col-lective chaos [19–25] as well as irregular activity states in het-erogeneous neuronal ensembles [30, 56].We thank Hugues Chat´e for valuable comments. DP ac-knowledges support by MINECO (Spain) under the Ram´ony Cajal programme. We acknowledge support by MINECO(Spain) under project No. FIS2014-59462-P, and by the Eu-ropean Union’s Horizon 2020 research and innovation pro-gramme under the Marie Skłodowska-Curie grant agreementNo. 642563. [1] X.-J. Wang, Physiological reviews , 1195 (2010).[2] A. T. Winfree, The Geometry of Biological Time (Springer, NewYork, 1980).[3] F. C. Hoppensteadt and E. M. Izhikevich,
Weakly connectedneural networks. (Spinger Verlag, N.Y., 1997).[4] P. Ashwin, S. Coombes, and R. Nicks, The Journal of Mathe-matical Neuroscience , 2 (2016).[5] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984).[6] A. Pikovsky and M. Rosenblum, Chaos , 097616 (2015).[7] C. van Vreeswijk, Phys. Rev. E , 5522 (1996).[8] P. K. Mohanty and A. Politi, J. Phys. A: Math. Gen. , L415(2006).[9] S. Olmi, A. Politi, and A. Torcini, EPL (Europhys. Lett.) ,60007 (2010).[10] S. Luccioli, S. Olmi, A. Politi, and A. Torcini, Phys. Rev. Lett. , 138103 (2012).[11] R. Burioni, S. di Santo, M. di Volo, and A. Vezzani, Phys. Rev.E , 042918 (2014).[12] A. Vilfan and T. Duke, Phys. Rev. Lett. , 114101 (2003).[13] M. Rosenblum and A. Pikovsky, Phys. Rev. Lett. , 064101(2007).[14] A. Pikovsky and M. Rosenblum, Physica D: Nonlinear Phe-nomena , 27 (2009).[15] A. A. Temirbayev, Z. Z. Zhanabaev, S. B. Tarasov, V. I. Pono-marenko, and M. Rosenblum, Phys. Rev. E , 015204 (2012).[16] A. A. Temirbayev, Y. D. Nalibayev, Z. Z. Zhanabaev, V. I. Pono-marenko, and M. Rosenblum, Phys. Rev. E , 062917 (2013).[17] M. Rosenblum and A. Pikovsky, Phys. Rev. E , 012919(2015).[18] A. Politi and M. Rosenblum, Phys. Rev. E , 042916 (2015).[19] V. Hakim and W. J. Rappel, Phys. Rev. A , R7347 (1992).[20] N. Nakagawa and Y. Kuramoto, Prog. Theor. Phys. , 313(1993).[21] N. Nakagawa and Y. Kuramoto, Physica D , 74 (1994).[22] N. Nakagawa and Y. Kuramoto, Physica D , 307 (1995).[23] K. A. Takeuchi, F. Ginelli, and H. Chat´e, Phys. Rev. Lett. ,154103 (2009).[24] K. A. Takeuchi, H. Chat´e, F. Ginelli, A. Politi, and A. Torcini, Phys. Rev. Lett. , 124101 (2011).[25] W. L. Ku, M. Girvan, and E. Ott, Chaos , 123122 (2015).[26] E. Ott and T. M. Antonsen, Chaos , 037113 (2008).[27] T. B. Luke, E. Barreto, and P. So, Neural Comput. , 3207(2013).[28] D. Paz´o and E. Montbri´o, Phys. Rev. X , 011009 (2014).[29] E. Montbri´o, D. Paz´o, and A. Roxin, Phys. Rev. X , 021028(2015).[30] S. Luccioli and A. Politi, Phys. Rev. Lett. , 158104 (2010).[31] B. Ermentrout and N. Kopell, SIAM J. Appl. Math. , 233(1986).[32] P. Latham, B. Richmond, P. Nelson, and S. Nirenberg, Journalof Neurophysiology , 808 (2000).[33] E. M. Izhikevich, Dynamical Systems in Neuroscience (TheMIT Press, Cambridge, Massachusetts, 2007).[34] See Supplementary Material [url] for details of the numericalimplementation.[35] See Supplementary Material [url] for the exact derivation of theFREs corresponding to Eqs. (1)-(2).[36] For identical neurons the dynamics of the model is degenerate(and described by the Watanabe-Strogatz theory [57]), but thepresence of a tiny amount of noise attracts the dynamics to theLorentzian manifold, making Eq. (4) with ∆ = 0 asymptot-ically correct [58]. Accordingly, the initial conditions for thenumerical simulations in Figs. 1 and 2 are taken to representa Lorentzian density with arbitrary values of the center v withhalf-width πr : V j (0) = v + πr tan[( π/ j − N − / ( N +1)] .[37] This can always be achieved (for ¯ η > ) under rescaling oftime ˜ t = t √ ¯ η/τ , the variables ˜ r = rτ / √ ¯ η, ˜ v = v/ √ ¯ η and theparameters ˜ J = J/ √ ¯ η , ˜ D = D √ ¯ η/τ , ˜∆ = ∆ / ¯ η .[38] See Supplementary Material [url]. [39] M. K. S. Yeung and S. H. Strogatz, Phys. Rev. Lett. , 648(1999).[40] M. Y. Choi, H. J. Kim, D. Kim, and H. Hong, Phys. Rev. E ,371 (2000).[41] E. Montbri´o, D. Paz´o, and J. Schmidt, Phys. Rev. E , 056201(2006).[42] W. S. Lee, E. Ott, and T. M. Antonsen, Phys. Rev. Lett. ,044101 (2009).[43] J. D. Farmer, Physica D , 366 (1982).[44] See Supplementary Material [url].[45] Supercriticality is also present for the line J (3) H in a tiny intervalabove D = 2 π .[46] X.-J. Wang and G. Buzs´aki, J. Neurosci. , 6402 (1996).[47] M. Whittington, R. Traub, N. Kopell, B. Ermentrout, andE. Buhl, Int. Journal of Psychophysiol. , 315 (2000).[48] M. Bartos, I. Vida, and P. Jonas, Nature reviews neuroscience , 45 (2007).[49] N. Brunel and V. Hakim, Chaos , 015113 (2008).[50] A. Roxin, N. Brunel, and D. Hansel, Phys. Rev. Lett. ,238103 (2005).[51] D. Battaglia, N. Brunel, and D. Hansel, Phys. Rev. Lett. ,238106 (2007).[52] A. Roxin and E. Montbri´o, Physica D , 323 (2011).[53] N. Brunel and V. Hakim, Neural Comput. , 1621 (1999).[54] N. Brunel and X.-J. Wang, Journal of neurophysiology , 415(2003).[55] N. Brunel and D. Hansel, Neural Comput. , 1066 (2006).[56] E. Ullner and A. Politi, Phys. Rev. X , 011015 (2016).[57] S. Watanabe and S. H. Strogatz, Physica D74