The Winfree model with non-infinitesimal phase-response curve: Ott-Antonsen theory
TThe Winfree model with non-infinitesimal PRC
The Winfree model with non-infinitesimal phase-response curve:Ott-Antonsen theory
Diego Pazó a) and Rafael Gallego Instituto de Física de Cantabria (IFCA), CSIC-Universidad de Cantabria, 39005 Santander,Spain Departamento de Matemáticas, Universidad de Oviedo, Campus de Viesques, 33203 Gijón,Spain (Dated: 28 July 2020)
A novel generalization of the Winfree model of globally coupled phase oscillators, representing phase reduction underfinite coupling, is studied analytically. We consider interactions through a non-infinitesimal (or finite) phase-responsecurve (PRC), in contrast to the infinitesimal PRC of the original model. For a family of non-infinitesimal PRCs, theglobal dynamics is captured by one complex-valued ordinary differential equation resorting to the Ott-Antonsen ansatz.The phase diagrams are thereupon obtained for four illustrative cases of non-infinitesimal PRC. Bistability betweencollective synchronization and full desynchronization is observed in all cases.
In 1967 Winfree proposed a model for the sponta-neous synchronization of large ensembles of biologicaloscillators . The Winfree model played a seminal role inthe field of collective synchrony, inspiring the Kuramotomodel as well as promoting recent advances in theoreti-cal neuroscience . In spite of the simplifying assumptionsof the Winfree model, uniform all-to-all weak coupling,analytical solutions have been found only recently usingthe Ott-Antonsen ansatz , see also . Weak coupling isimplicit in the use of phase oscillators as the units of themodel. Moreover, their interactions are modeled by the so-called infinitesimal phase-response curve (iPRC), whichis only valid in the limit of vanishing coupling. In thispaper we extend the Winfree model considering a non-infinitesimal (also called finite) PRC, such that the phaseshift of one oscillator is not proportional to the magni-tude of the input. For a family of non-infinitesimal PRCs,and a Lorentzian distribution of natural frequencies, theglobal dynamics is captured by one complex-valued ordi-nary differential equation by means of the Ott-Antonsenansatz . We obtain the phase diagrams for four instruc-tive cases. I. INTRODUCTION
Collective synchronization in large ensembles of self-sustained oscillators is a pervasive phenomenon in nature andtechnology . The first successful attempt to model collec-tive synchronization is due to Winfree . Relying on his intu-ition he devised a model where the only degrees of freedomwere the oscillators’ phases, and the coupling was uniformand global (i.e. mean-field type). In the numerical simulationsa macroscopic cluster of synchronized oscillators emergedspontaneously when either the natural frequencies of the os-cillators were narrowly distributed or the coupling was large a) Author to whom correspondence should be addressed: [email protected] enough. In mathematical language, the phases in the Winfreemodel are governed by a set of N ordinary differential equa-tions ( i = , . . . , N ): ˙ θ i = ω i + ˜ Q ( θ i ) A , (1a) A = ε N N ∑ j = P ( θ j ) . (1b)Here ω i is the natural frequency of the i -th oscillator, and ε > π -periodic function P specifies the pulse shape. The function ˜ Q is also 2 π -periodic and is either called infinitesimal (or linear)phase-response curve (iPRC), or sensitivity function .As already mentioned, the Winfree model relies on two as-sumptions: weak coupling and all-to-all geometry. Weak cou-pling permits, first of all, ignoring the oscillators’ amplitudes:the limit cycles are strongly attracting compared to perturba-tions causing amplitudes to be strongly damped degrees offreedom. In addition, the effect of the mean field A on thephase is exactly proportional to A —higher powers of A areabsent in Eq. (1a)—, which only holds in the limit of asymp-totically small interactions .In this work we generalized the Winfree model consider-ing nonlinear interactions. Mathematical tractability imposescertain restrictions on the distribution of the natural frequen-cies and on the class of “non-infinitesimal” (also called “fi-nite” or “non-linear”) PRCs, but we believe it is remark-able that such analytic solutions exist. This limited progressshould be welcome given the relevance of PRC theory in the-oretical neuroscience , and recent experiments evidencingthe insufficiency of the linear approximation . Our anal-ysis is based on the so-called “Ott-Antonsen (OA) theory”,which assumes a certain ansatz (the Poisson kernel) for thedensity of the phases in the thermodynamic limit ( N → ∞ ).The OA ansatz was initially applied to the Kuramoto modeland its variants , but eventually found application in sev-eral systems of pulse-coupled oscillators: the original Win-free model (and a variant with heterogeneous iPRCs ),ensembles of theta neurons , quadratic integrate-and-fireneurons , and excitable active rotators . a r X i v : . [ n li n . AO ] J u l he Winfree model with non-infinitesimal PRC 2 II. WINFREE MODEL WITH NON-INFINITESIMAL PRC
We consider a modification of the Winfree model (1), inwhich Eq. (1a) is replaced by˙ θ i = ω i + Q ( θ i , A ) , i = , . . . , N , (2)where A is the mean field defined by (1b). At the lowest or-der in A , the model (2) converges to the Winfree model (1): dQ ( θ i , A ) / dA | A = = ˜ Q ( θ i ) . Assuming Q ( θ , A ) linear in A isequivalent to approximate the isochrons of a limit cycle bystraight lines (or hyperplanes if the dimensionality is largerthan two) in the phase reduction procedure . A. Non-infinitesimal PRC
Prior to specifying the PRC Q , we devote a few lines to theiPRCs. Traditionally, iPRCs are classified as type I or typeII . For type II, either an advance or a delay in the phase arepossible depending upon the timing of the perturbation, whilein the case of type I the timing of the perturbation does notchange the sign of the phase shift. The canonical examplesof each type are ˜ Q ( θ ) ∝ − cos θ for type I (e.g. thetheta neuron), and ˜ Q ( θ ) ∝ sin θ for type II (e.g. the Stuart-Landau oscillator). For non-infinitesimal PRCs, the previousclassification falls short as the character of Q may change withthe strength of the stimulus .The types of PRC we consider are conditioned by the appli-cability of the OA ansatz, as it enables a drastic dimensionalityreduction. The OA ansatz imposes that no harmonics in θ be-yond the first one are present in Q ( θ , A ) . Still, the family ofPRCs with only first harmonic in θ is wide enough to makethe problem nontrivial. As we shall adopt pulses P ( θ ) withpeak value at θ = π ), we impose the ad-ditional constraint Q ( , A ) = and certain fireflies . Therefore, we restrict to a family ofPRCs of this form: Q ( θ , A ) = f ( A )( − cos θ ) − f ( A ) sin θ , (3)where f and f are arbitrary functions of A , provided that f , ( ) = ( − cos θ ) and sin θ , as the type-I and thetype-II components of the PRC, respectively. B. Pulse shape
In the study of the classical Winfree model several pulseshapes can be considered, see . In this work, we adopt a “rec-tified Poisson kernel” : P ( θ ) = ( − r )( + cos θ ) − r cos θ + r . (4)This is a particularly convenient shape for the theoretical anal-ysis below. P ( θ ) is a symmetric unimodal function in the in-terval [ − π , π ] (with the normalization (cid:82) π − π P ( θ ) d θ = π ) that − π π θ P ( θ ) r = 0 . r = 0 . r = 0 . FIG. 1. Rectified-Poisson pulse (4) for three values of parameter r . peaks at θ = θ = ± π . Parameter r is areal number allowing a continuous interpolation between a flatpulse for r = − P ( θ ) = πδ ( θ ) , for r =
1. In Fig. 1 the pulse function P ( θ ) is depicted for threedifferent values of r . C. Natural frequencies
For the sake of achieving the maximal dimensionality re-duction, we assume the natural frequencies to be distributedaccording to a Lorentzian distribution of half-width ∆ cen-tered at ω : g ( ω ) = ∆ / π ( ω − ω ) + ∆ . (5) III. OTT-ANTONSEN THEORY
Once the building blocks of the model have been intro-duced, we apply the OA theory . In this way we derivea complex-valued ODE reproducing the long-time evolutionof the model at the macroscopic level. As the procedure isstandard , the readers interested in the final result are pointedto Eqs. (10) and (11).First of all, one must realize that our model (2) belongs to ageneral class of oscillator systems of the form˙ θ i ( t ) = ω i + B ( t ) + Im (cid:104) H ( t ) e − i θ i ( t ) (cid:105) , (6)which can be analyzed with the OA ansatz . Functions B and H may depend explicitly on time or indirectly througha mean field. For the PRC (3) we have B ( t ) = f ( A ) , H ( t ) = f ( A ) − i f ( A ) . (7)In the thermodynamic limit we can define a phase density F ( θ | ω , t ) , such that F ( θ | ω , t ) d θ is the fraction of oscillatorsof frequency ω at time t , with phases in the interval [ θ , θ + he Winfree model with non-infinitesimal PRC 3 d θ ] . It is convenient to introduce the Fourier expansion of thedensity F ( θ | ω , t ) = ∞ ∑ m = − ∞ α m ( ω , t ) e im θ , with α − m = α ∗ m . We notice as well that, by conservation ofthe number of oscillators, F satisfies the continuity equation: ∂ t F + ∂ θ ( F ˙ θ ) =
0, where ˙ θ is the speed of an oscillator ofnatural frequency ω . Inserting the Fourier series of F into thecontinuity equation we get: ∂ t α m ( ω , t ) = − im ( ω + B ) α m + m ( H ∗ α m − − H α m + ) . (8)A particular solution of this equation, the OA ansatz, is ob-tained equating the coefficient of m -th mode to the m -th powerof the first mode: α m = α m . Hence, for the solution in this so-called OA manifold, we only need to consider the evolutionof α ≡ α : ∂ t α ( ω , t ) = − i ( ω + B ) α + (cid:0) H ∗ − H α (cid:1) . (9)This is still an infinite set of coupled integro-differential equa-tions. A sharp reduction in the dimensionality of the problemis achieved for rational g ( ω ) , and specially for the Lorentziandistribution . As the Kuramoto order parameter Z = e i θ isrelated to α via Z ∗ ( t ) = (cid:82) ∞ − ∞ α ( ω , t ) g ( ω ) d ω , we can eval-uate this integral resorting to the residue theorem obtaining Z ∗ ( t ) = α ( ω − i ∆ , t ) . (This is the result of performing ananalytic continuation of α from real to complex ω , and eval-uating α at the pole of g ( ω ) in the lower half ω -plane.) Thus,setting ω = ω − i ∆ in (9), we get a complex-valued ODE forthe Kuramoto order parameter:˙ Z = ( − ∆ + i ω ) Z − i f ( A )( − Z ) + f ( A )( − Z ) , (10)where B and H have been written in terms of f and f ac-cording to Eq. (7). To close Eq. (10) we need to express themean field A as a function of Z . For the pulse shape in Eq. (4)and a Lorentzian frequency distribution it can be proven (see or the supplemental material of ) that: A = ε Re (cid:18) + Z − r Z (cid:19) . (11)Note that 0 ≤ A ≤ A max , where the maximal value A max = ε / ( − r ) is achieved if Z = θ j = ω ishereafter set to 1, as this can always be achieved by rescalingtime and f , . IV. FOUR ILLUSTRATIVE PRCs
Among the infinite set of functions f ( A ) and f ( A ) we se-lected a few illustrative case studies. In each of these cases thecharacter of the PRC undergoes a crossover as A grows: from TABLE I. The four cases of non-infinitesimal PRCs analyzed in thispaper. Functions f and f determine the PRC in (3) ( σ ( A ) is acrossover function, see (12)). The code in the last column indicatesthe iPRC and the asymptotic PRC at large A (see text).Case f ( A ) f ( A ) Codea σ ( A ) A σ ( A ) I − II s b A σ ( A ) σ ( A ) II s − Ic 0 ( − A ) σ ( A ) II s − II r d 0 − ( − A ) σ ( A ) II r − II s − − (a) I − II s Q ( A , θ ) / A A = 0 (iPRC) A = 1 / A = 2 A = 8 (b) II s − I0 π π − − (c) II s − II r θ Q ( A , θ ) / A π π (d) II r − II s θ FIG. 2. The non-infinitesimal PRCs analyzed in this work as a func-tion of θ for four representative values of A , including A = + (iPRC)and A = one iPRC type to a different PRC type for large A . We de-note the limiting PRC at A → ∞ as ‘asymptotic PRC’ (aPRC).From now on, we apply the classical distinction between typesI and II to both iPRCs and aPRCs. Recall that if the sign isthe same for all θ we call the iPRC (or the aPRC) type I (im-plying f = f = f . In turn, we dis-tinguish between two subclasses of the type II: II s ( f >
0) andII r ( f <
0) corresponding to the synchronizing and repulsiveinteractions, respectively.As we are interested in introducing one crossover in thePRC between the iPRC and the aPRC, and have three funda-mentally different types (I, II s , and II r ) this gives 6 possiblecombinations. However, we shall consider only four of theseiPRC-aPRC pairs, since only type II s favors synchrony and isto be included either in the iPRC or in the aPRC. Otherwiseno synchronization phenomena are expected: type I is neutraland type II r is repulsive. Hence, we focus on the four caseslisted in Table I, in which different PRC types characterizesmall and large A regimes. As a guide, in the fourth columnhe Winfree model with non-infinitesimal PRC 4 Time (a)
Time (c) (cid:0) (cid:0) Re .Z/ (cid:0) (cid:0) I m . Z / (b) (cid:0) (cid:0) Re .Z/ (cid:0) (cid:0) I m . Z / (d) FIG. 3. (a) Raster plot —a dot indicates at which time one oscillatorphase crosses a multiple of 2 π — for a population of N =
500 and thenon-infinitesimal PRC of case d, see Table I and Fig. 2(d). The initialcondition is uniform θ j ( t = ) = .
01 and parameters are ∆ = . r = . ε = .
4. The frequencies are deterministically drawn froma Lorentzian distribution: ω i = ω + ∆ tan [ π ( i − N − ) / ( N )] . (c)The same as (a) but for a random initial distribution of phases. (b)and (d) depict the Kuramoto order parameter Z ( t ) = N − ∑ j e i θ j ( t ) for 50 t.u., once the simulations in (a) and (c), respectively, reachedthe stationary state. The red dashed line and the red cross in panels(b) and (d) are the periodic and fixed point attractors of Eq. (10),coexisting at the same parameter values. of the Table we write a code X-Y, where X refers to the iPRCand Y to the aPRC. The saturation function σ ( A ) in the Tablehas positive slope at A =
0, and saturates at large A . In partic-ular, we chose this specific saturation function in our study: σ ( A ) = A + A . (12)(Our results have been occasionally tested against anotherchoice σ ( A ) = tanh ( A ) , finding no qualitative difference.)Graphical representations of the four PRCs (cases a to d), forfour representative A values, are shown in Figs. 2(a)-(d). Ineach panel the PRC appears divided by A as usual , and thelack of overlapping between different lines evidences its non-linearity.In the next section we obtain the phase diagrams corre-sponding to each of the four cases introduced here, based onthe analysis of the complex-valued ordinary differential equa-tion (10). But before doing so, it is worth making direct sim-ulations of the full system (2), and test (and understand) thecorrespondence with the solutions of Eq. (10). We simulatedthe full model in case d with ω =
1, heterogeneity param-eter ∆ = .
01, pulse-shape parameter r = . ε = .
4. As may be seen in Fig. 3, the populationexhibits bistability between a desynchronized state and a syn-chronized state with some oscillators oscillating with the samefrequency. This bistability is not surprising as the system is“more synchronizing” when already synchronized since theaPRC is of type II s , while it is hardly synchronizable when al-ready desynchronized by virtue of the type II r iPRC. In termsof Z , the synchronous solution is (approximately) a periodic orbit, while the desynchronized state exhibits only small fluc-tuations around a point due to finite size effects ( N = V. PHASE DIAGRAMS
In the remainder of this paper we obtain the phase diagramsfor the four reference cases by means of Eq. (10). As (10)is a (generic) planar system, the only possible attractors arefixed points and limit cycles. Their bifurcation loci, depictedin the phase diagrams below, have been obtained using the
MATCONT toolbox of MATLAB . Moreover, we recall that Z is only physically meaningful inside the unit disk | Z | ≤
1, andtherefore attractors and bifurcations occurring outside it areignored. As seen in Fig. 3, limit cycles correspond to synchro-nized solutions, in which a macroscopic part of the populationrotates at the same average frequency.
A. Case a: I − II s In Fig. 4(a) we show the phase diagram spanned by param-eters ∆ and ε . Bifurcation lines for three values of parameter r , controlling the pulse width, are depicted. The results al-most replicate those in Ref. for the standard Winfree modelwith type II s iPRC. Synchronization is found in two adjacentregions, in one of them (dark shaded) coexisting with a desyn-chronized state. (There exists a region (not shown) besidesthe bistability region where two desynchronized states coex-ist, see ). In contrast to the averaging approximation (theKuramoto-Sakaguchi model), valid at small ε and ∆ , synchro-nization becomes impossible if the population is too hetero-geneous (large ∆ ).For small coupling (and heterogeneity), synchronizationemerges from a supercritical Hopf bifurcation undergone bythe desynchronized state, akin to the classical Kuramototransition . This Hopf bifurcation line terminates at a doublezero eigenvalue (Bogdanov-Takens, BT) point. A homoclinic(Hom) line emanates from the BT point limiting the coexis-tence region. As observed for the regular Winfree model synchronization is more efficient for narrow pulses. The pulsewidth does not qualitatively change the phase diagram.The phase diagram only differs appreciably from those in at the origin. We see that, due to the type-I iPRC, the Hopfline approaches the origin with an infinite slope. In particular,the asymptotic dependence of the critical ε H on ∆ follows anunusual square-root law with the frequency dispersion ∆ : ε H = (cid:114) ∆ + r . (13)We can deduce this result deriving the associated Kuramoto-Sakaguchi model of model (2) via averaging. Or, alternatively,preserving in (10) only linear, rotationally invariant terms in Z , and equating the linear coefficient to i Ω .he Winfree model with non-infinitesimal PRC 5 Synchronousstate SNICSNSNSL (a) (b)(d)(c)
Hopf supHopf subHopf supHopf subHopf sup
SNLCTC BTSNSLTC
Hopf supHopf sub
SNICSNSNLCHomHom SNSL
FIG. 4. Synchronization regions of the Winfree model in the ( ∆ , ε ) -plane for the four cases of PRCs described in Table I, and three differentvalues of the parameter r ∈ { . , . , . } . Panels (a), (b), (c) and (d) correspond to cases a, b, c, and d, respectively. For the value r = .
9, lightshaded regions indicate where there is a stable limit cycle, corresponding to a macroscopic synchronized state. In the dark shaded regions, thelimit cycle (synchronous state) coexists with a stable fixed point (asynchronous state). The codimension-two points are depicted by specificsymbols: Generalized Hopf (GH- (cid:78) ), saddle-node separatrix-loop (SNSL- (cid:4) ), Bogdanov-Takens (BT- • ), and transcritical (TC- (cid:70) ) bifurcations.The bifurcation corresponding to each line type is indicated in the legend of the respective panel. Insets in panels (c) and (d) are magnificationsof the regions inside the respective rectangles. B. Case b: II s − I In case b, iPRC and aPRC are interchanged with respect tocase a. This means that synchronization is favored at smallcoupling, but becomes increasingly difficult as the couplinggrows. Accordingly, the phase diagram in Fig. 4(b) showsthe expected supercritical Hopf bifurcation line emanating asa straight line from the origin : ε H ∝ ∆ + O ( ∆ ) .At large ε there is a bistability region such that the synchro-nized state disappears in a saddle-node bifurcation of limit cy-cles (SNLC). The locus of the SNLC is a line that emanatesfrom a generalized Hopf (or Bautin) point (GH), and termi-nates at the ε axis at a point marked with a star on the ε -axis ofthe phase diagram. The stars pinpoint the (equivariant) tran-scritical (TC) bifurcation , in which the fully synchronizedstate ( θ i ( t ) = θ j ( t ) ) of identical oscillators ( ∆ =
0) becomes unstable. For r = .
9, the instability of full synchronizationtakes place at ε c = . . . . , far above the range of ε displayedin the phase diagram. The location of ε c was not calculated us-ing (10), but by directly looking for the stability threshold tothe fully synchronized state, see Appendix.Finally, note that the synchronization region shrinks as thepulse becomes wider, but there is not a qualitative change inthe phase diagram whatsoever. C. Case c: II s − II r In this case the aPRC is repulsive, in contrast to case bwhere the aPRC is type I (i.e. neutral in terms of synchro-nization). In turn the phase diagram in Fig. 4(c) shows aquite small synchronization region (notice the scale of thehe Winfree model with non-infinitesimal PRC 6axes). Synchronization is bounded exclusively by a supercrit-ical Hopf bifurcation, save for broad pulses. In the latter casea GH point is found, and the Hopf bifurcation is subcritical atthe left of it. Accordingly, we find a bistability region boundedby a line of saddle-node bifurcation of limit cycles (SNLC)and a subcritical Hopf bifurcation, as in case b. The precisevalue of r below which the bistability region exists (i.e. theGH point is present) is r ∗ (cid:39) . ∆ =
0, above which full synchrony destabilizes . The tran-scritical bifurcation is not structurally stable, see e.g. Fig. 11in , and increasing ∆ from 0 may either leave no trace of bi-furcation or “decay” into two saddle-node bifurcation of limitcycles. The latter scenario occurs for r < . . . . , see thebifurcation lines for r = . | Z | > D. Case d: II r − II s Case d exhibits the most complex phase diagram amongall those obtained here. The aPRC is of type II s , as in casea, and (accordingly) the large ε region is organized by twocodimension-two points: The Bogdanov-Takens (BT), and thesaddle-node separatrix-loop (SNSL) codimension-two points.The associated region of bistability between synchrony andasynchrony is bounded by homoclinic, saddle-node and Hopfbifurcations.Remarkably, there is also a bistability region at small ε val-ues for r > r ∗ (cid:39) . r → VI. CONCLUSIONS
In this work we have studied a non-trivial extension of theWinfree model in which the PRC is nonlinear in the meanfield. If the PRC contains only the first harmonic of the an-gle, the OA ansatz permits a sharp dimensionality reduction.Among all possible dependencies of the PRC on the meanfield, we have considered only those with a crossover be-tween two different canonical components. In particular, wehave analyzed four cases in which an attractive type-II com-ponent competes either against a repulsive type-II componentor against a type-I component. Synchronization regions arepeculiar for each case. Bistability between macroscopic syn-chronization and complete desynchronization are found in allcases (in case c, only for broad pulses), but in different relativelocations in the ∆ − ε plane.Our results indicate that the nonlinearity of the PRC withthe forcing, by itself, is not enough to generate complexcollective phenomena. This is certain for a Lorentzian dis-tribution of frequencies since the reduced system is only two dimensional, irrespective of the exact form of f ( A ) and f ( A ) . As happens in Kuramoto-like models, phe-nomena such as clustering or glassy dynamics may requiremultiple Fourier components (in the PRC) or strongerheterogeneity , respectively. Concerning collective chaos,other ingredients such as a time-varying coupling , two inter-acting populations or multimodal frequency distributions appear to be imperative.Needless to say, our study is only a drop in the ocean ofpossible PRCs and model generalizations. For instance, relax-ation oscillators and bursting (neuronal) oscillators havePRCs very different from the first-harmonic shape function inEq. (3). Nevertheless, in spite of its limitations, we regard themodel defined by Eqs. (2) and (3) as a noteworthy example ofsystem in which the OA theory can be fully applied. ACKNOWLEDGMENTS
We acknowledge support by the Agencia Estatal de Inves-tigación and Fondo Europeo de Desarrollo Regional underProject No. FIS2016-74957-P (AEI/FEDER, EU).
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
APPENDIX: IDENTICAL OSCILLATORS
If the oscillators are identical there is a fully synchronizedsolution θ j ( t ) = Ψ ( t ) . The dynamics of Ψ obeys:˙ Ψ = ω + f [ ε P ( Ψ )]( − cos Ψ ) − f [ ε P ( Ψ )] sin Ψ . Next, we calculate the stability threshold of full synchrony,fixing ω = N → ∞ ), we may perturb one oscillator, say the first one,without changing the mean field. Hence, one infinitesimalperturbation δ θ = θ − Ψ obeys:˙ δ θ = λ ( Ψ ) δ θ , where the multiplicative factor λ ( Ψ ) = f [ ε P ( Ψ )] sin Ψ − f [ ε P ( Ψ )] cos Ψ depends on time through Ψ ( t ) . In order toknow the average exponential growth (or contraction) rate of δ θ we need to integrate over variable Ψ , taking into accountits density ρ ( Ψ ) . These means that the sign of constant λ ,given by λ = (cid:90) π − π λ ( Ψ ) ρ ( Ψ ) d Ψ , determines the stability of the fully synchronized solution. If λ is positive, the oscillator “evaporates” from the main cluster,i.e. full synchrony is unstable.he Winfree model with non-infinitesimal PRC 7The density ρ ( Ψ ) is proportional to the inverse of thespeed: ρ ( Ψ ) ∝ ˙ Ψ − . Imposing λ =
0, we obtain the conditionfor the stability threshold of full synchrony: (cid:90) π − π f [ ε c P ( Ψ )] sin Ψ − f [ ε c P ( Ψ )] cos Ψ + f [ ε c P ( Ψ )]( − cos Ψ ) − f [ ε c P ( Ψ )] sin Ψ d Ψ = . This integral cannot be solved analytically, but the thresholdcoupling ε c is easily found numerically. REFERENCES A. T. Winfree, “Biological rhythms and the behavior of populations of cou-pled oscillators.” J. Theor. Biol. , 15–42 (1967). Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscil-lators,” in
International Symposium on Mathematical Problems in Theoreti-cal Physics , Lecture Notes in Physics, Vol. 39, edited by H. Araki (Springer,Berlin, 1975) pp. 420–422. S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of syn-chronization in populations of coupled oscillators,” Physica D , 1–20(2000). E. Montbrió and D. Pazó, “Kuramoto model for excitation-inhibition-basedoscillations,” Phys. Rev. Lett. , 244101 (2018). D. Pazó and E. Montbrió, “Low-dimensional dynamics of populations ofpulse-coupled oscillators,” Phys. Rev. X , 011009 (2014). R. Gallego, E. Montbrió, and D. Pazó, “Synchronization scenarios in theWinfree model of coupled oscillators,” Phys. Rev. E , 042208 (2017). J. T. Ariaratnam and S. H. Strogatz, “Phase diagram for the Winfree modelof coupled nonlinear oscillators,” Phys. Rev. Lett. , 4278–4281 (2001). E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems ofglobally coupled oscillators,” Chaos , 037113 (2008). E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator sys-tems,” Chaos , 023117 (2009). E. Ott, B. R. Hunt, and T. M. Antonsen, “Comment on “long time evolutionof phase oscillators systems”,” Chaos , 025112 (2011). A. T. Winfree,
The Geometry of Biological Time (Springer, New York,1980). A. S. Pikovsky, M. G. Rosenblum, and J. Kurths,
Synchronization, a Uni-versal Concept in Nonlinear Sciences (Cambridge University Press, Cam-bridge, 2001). E. M. Izhikevich,
Dynamical Systems in Neuroscience (The MIT Press,Cambridge, Massachusetts, 2007) Chap. 10. Y. Kuramoto,
Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984). P. Sacré and R. Sepulchre, “Sensitivity analysis of oscillator models in thespace of phase-response curves: Oscillators as open systems,” IEEE Con-trol Systems Magazine , 50–74 (2014). B. Pietras and A. Daffertshofer, “Network dynamics of coupled oscillatorsand phase reduction techniques,” Phys. Rep. , 1 – 105 (2019), networkdynamics of coupled oscillators and phase reduction techniques. G. B. Ermentrout and D. H. Terman,
Mathematical Foundations of Neuro-science , Vol. 64 (Springer, 2010). C. Börgers,
An Introduction to Modeling Neuronal Dynamics , Texts in Ap-plied Mathematics (Springer International Publishing, 2017). J. Rode, J. F. Totz, E. Fengler, and H. Engel, “Chimera states on a ring ofstrongly coupled relaxation oscillators,” Frontiers in Applied Mathematicsand Statistics , 31 (2019). D. C˘alug˘aru, J. F. Totz, E. A. Martens, and H. Engel, “First-order syn-chronization transition in a large population of relaxation oscillators,”arXiv:1812.04727. D. Pazó, E. Montbrió, and R. Gallego, “The Winfree model with hetero-geneous phase-response curves: analytical results,” J. Phys. A: Math. andTheor. , 154001 (2019). T. B. Luke, E. Barreto, and P. So, “Complete classification of the macro-scopic behavior of a heterogeneous network of theta neurons,” Neural Com-put. , 3207–3234 (2013). C. R. Laing, “Derivation of a neural field model from a network of thetaneurons,” Phys. Rev. E , 010901 (2014). P. So, T. B. Luke, and E. Barreto, “Networks of theta neurons with time-varying excitability: Macroscopic chaos, multistability, and final-state un-certainty,” Physica D , 16–26 (2014). E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for net-works of spiking neurons,” Phys. Rev. X , 021028 (2015). D. Pazó and E. Montbrió, “From quasiperiodic partial synchronization tocollective chaos in populations of inhibitory neurons with delay,” Phys. Rev.Lett. , 238101 (2016). I. Ratas and K. Pyragas, “Macroscopic self-oscillations and aging transitionin a network of synaptically coupled quadratic integrate-and-fire neurons,”Phys. Rev. E , 032215 (2016). K. P. O’Keeffe and S. H. Strogatz, “Dynamics of a population of oscillatoryand excitable elements,” Phys. Rev. E , 062203 (2016). J. Roulet and G. B. Mindlin, “Average activity of excitatory and inhibitoryneural populations,” Chaos , 093104 (2016). D. Hansel, G. Mato, and C. Meunier, “Synchrony in excitatory neural net-works,” Neural Comput. , 307–337 (1995). A. D. Reyes and E. E. Fetz, “Two modes of interspike interval shortening bybrief transient depolarizations in cat neocortical neurons,” J. Neurophysiol. , 1661–1672 (1993). T. I. Netoff, M. I. Banks, A. D. Dorval, C. D. Acker, J. S. Haas, N. Kopell,and J. A. White, “Synchronization in hybrid neuronal networks of the hip-pocampal formation,” J. Neurophysiol. , 1197–1208 (2005). J. Buck, “Synchronous rhythmic flashing of fireflies II,” Q. Rev. Biol. ,265–289 (1988). F. E. Hanson, “Comparative studies of firefly pacemakers,” Fed. Proc. ,2158–2164 (1978). B. Pietras and A. Daffertshofer, “Ott-Antonsen attractiveness for parameter-dependent oscillatory systems,” Chaos , 103101 (2016). A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “MATCONT: A matlabpackage for numerical bifurcation analysis of ODEs,” SIGSAM Bull. ,21–22 (2004). P. Ashwin and J. W. Swift, “The dynamics of n weakly coupled identicaloscillators,” J. Nonlin. Sci. , 69–108 (1992). As the “OA manifold” is not attracting for identical oscillators , the re-sulting dynamics depends on the initial conditions. (For initial conditionsin the OA manifold, such as purely random phases ( Z = , a state which Z oscillates periodically while the individual oscillators behave quasiperi-odically.). J. D. Crawford, “Introduction to bifurcation theory,” Rev. Mod. Phys. ,991–1037 (1991). K. Okuda, “Variety and generality of clustering in globally coupled oscilla-tors,” Physica D: Nonlinear Phenomena , 424 – 436 (1993). D. Iatsenko, P. V. E. McClintock, and A. Stefanovska, “Oscillator glassin the generalized Kuramoto model: synchronous disorder and two-steprelaxation,” Nat. Commun. , 4188 (2014). P. So and E. Barreto, “Generating macroscopic chaos in a network of glob-ally coupled phase oscillators,” Chaos , 033127 (2011). C. Bick, M. J. Panaggio, and E. A. Martens, “Chaos in Kuramoto oscillatornetworks,” Chaos , 071102 (2018). H. Cheng, S. Guo, Q. Dai, H. Li, and J. Yang, “Collective chaos and period-doubling bifurcation in globally coupled phase oscillators,” Nonlinear Dyn. , 2273–2281 (2017). P. Sacré and A. Franci, “Singularly perturbed phase response curves for re-laxation oscillators,” in (2016) pp. 4680–4685. W. E. Sherwood and J. Guckenheimer, “Dissecting the phase response of amodel bursting neuron,” SIAM J. Appl. Dyn. Syst. , 659–703 (2010). C. van Vreeswijk, “Partial synchronization in populations of pulse-coupledoscillators,” Phys. Rev. E , 5522–5537 (1996). A. Politi and M. Rosenblum, “Equivalence of phase-oscillator andintegrate-and-fire models,” Phys. Rev. E91