The Winfree model with non-infinitesimal phase-response curve: Ott-Antonsen theory

A novel generalization of the Winfree model of globally coupled phase oscillators, representing phase reduction under finite coupling, is studied analytically. We consider interactions through a non-infinitesimal (or finite) phase-response curve (PRC), in contrast to the infinitesimal PRC of the original model. For a family of non-infinitesimal PRCs, the global 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 between collective synchronization and full desynchronization is observed in all cases.


I. INTRODUCTION
][13] The first successful attempt to model collective synchronization is due to Winfree. 1 Relying on his intuition, he devised a model where the only degrees of freedom were the oscillators' phases and the coupling was uniform and global (i.e., mean-field type).In the numerical simulations, a macroscopic cluster of synchronized oscillators emerged spontaneously when either the natural frequencies of the oscillators were narrowly distributed or the coupling was large enough.In mathematical language, the phases in the Winfree model are governed by a set of N ordinary differential equations (i = 1, . . ., N), θi = ω i + Q(θ i ) A, (1a) Here, ω i is the natural frequency of the ith oscillator and ε > 0 is a parameter controlling the coupling strength.The 2π -periodic function P specifies the pulse shape.Function Q is also 2π -periodic and is either called infinitesimal (or linear) phase-response curve (iPRC) or sensitivity function. 11,13,14s already mentioned, the Winfree model relies on two assumptions: weak coupling and all-to-all geometry.First of all, weak coupling permits ignoring the oscillators' amplitudes: the limit cycles are strongly attractive compared to perturbations, causing amplitudes to be strongly damped degrees of freedom.][15][16] In this work, we generalized the Winfree model considering nonlinear interactions.Mathematical tractability imposes certain restrictions on the distribution of the natural frequencies and on the class of "non-infinitesimal" (also called "finite" or "non-linear") PRCs, but we believe it is remarkable that such analytic solutions exist.This limited progress should be welcome given the relevance of the PRC theory in theoretical neuroscience, 17,18 and recent experiments evidence the insufficiency of the linear approximation. 19,20ur analysis is based on the so-called "Ott-Antonsen (OA) theory," which assumes a certain ansatz (the Poisson kernel) for the density of the phases in the thermodynamic limit (N → ∞).The OA ansatz was initially applied to the Kuramoto model and its variants 8,9 but eventually found application in several systems of pulse-coupled oscillators: the original Winfree model 5,6 (and a variant with heterogeneous iPRCs 21 ), ensembles of theta neurons, [22][23][24] quadratic integrate-and-fire neurons, [25][26][27] and excitable active rotators. 28,29

II. WINFREE MODEL WITH NON-INFINITESIMAL PRC
We consider a modification of the Winfree model ( 1), in which Eq. ( 1a) is replaced by where A is the mean field defined by (1b).At the lowest order in A, model (2) converges to the Winfree model (1): Assuming Q(θ , A) to be linear in A is equivalent to approximate the isochrons of a limit cycle by straight lines (or hyperplanes if the dimensionality is larger than two) in the phase reduction procedure. 13,14 Non-infinitesimal PRC Prior to specifying the PRC Q, we devote a few lines to iPRCs.Traditionally, iPRCs are classified as type I or type II. 30For type II, either an advance or a delay in the phase are possible depending upon the timing of the perturbation, while in the case of type I, the timing of the perturbation does not change the sign of the phase shift.The canonical examples of each type 13,15 are Q(θ ) ∝ 1 − cos θ for type I (e.g., the theta neuron) and Q(θ ) ∝ sin θ for type II (e.g., the Stuart-Landau oscillator).For non-infinitesimal PRCs, the previous classification falls short as the character of Q may change with the strength of the stimulus. 13,15he types of PRC we consider are conditioned by the applicability of the OA ansatz, as it enables a drastic dimensionality reduction.The OA ansatz imposes that no harmonics in θ beyond the first one are present in Q(θ , A).Still, the family of PRCs with only first harmonic in θ is wide enough to make the problem nontrivial.As we shall adopt pulses P(θ ) with peak value at θ = 0 (and multiples of 2π ), we impose the additional constraint Q(0, A) = 0 motivated by the fact that the PRC vanishes at spiking/flashing times for most neurons 31,32 and certain fireflies. 33,34Therefore, we restrict to a family of PRCs of this form, where f 1 and f 2 are arbitrary functions of A, provided f 1,2 (0) = 0 for obvious physical reasons.In similarity with the classification of iPRCs, we refer to the two terms in (3), proportional to (1 − cos θ ) and sin θ , as the type I and the type II components of the PRC, respectively.

B. Pulse shape
In the study of the classical Winfree model, several pulse shapes can be considered, see Ref. 6.In this work, we adopt a "rectified-Poisson kernel", 6 This is a particularly convenient shape for the theoretical analysis below.P(θ ) is a symmetric unimodal function in the interval [−π , π ] [with the normalization π −π P(θ )dθ = 2π ] that peaks at θ = 0 and vanishes at θ = ±π .Parameter r is a real number allowing a continuous interpolation between a flat pulse for r = −1 and a Dirac-delta pulse, P(θ ) = 2πδ(θ ), for r = 1.In Fig. 1, the pulse function P(θ ) is depicted for three different values of r.

C. Natural frequencies
For the sake of achieving the maximal dimensionality reduction, we assume the natural frequencies to be distributed according to a Lorentzian distribution of half-width centered at ω 0 , (5)

III. OTT-ANTONSEN THEORY
Once the building blocks of the model have been introduced, we apply the OA theory. 8In this way, we derive a complex-valued ODE reproducing the long-time evolution of the model at the macroscopic level.As the procedure is standard, 6,8 the readers interested in the final result are pointed to Eqs. (10) and (11).
First of all, one must realize that our model (2) belongs to a general class of oscillator systems of the form

ARTICLE scitation.org/journal/cha
0]35 Functions B and H may depend explicitly on time or indirectly through a mean field.
For PRC (3), we have In the thermodynamic limit, we can define a phase density F(θ |ω, t) such that F(θ |ω, t)dθ is the fraction of oscillators of frequency ω at time t, with phases in the interval [θ , θ + dθ ].It is convenient to introduce the Fourier expansion of the density, We notice as well that, by conservation of the number of oscillators, F satisfies the continuity equation: where θ is the speed of an oscillator of natural frequency ω.Inserting the Fourier series of F into the continuity equation, we get A particular solution of this equation, the OA ansatz, is obtained equating the coefficient of mth mode to the mth power of the first mode: Hence, for the solution in this so-called OA manifold, we only need to consider the evolution of This is still an infinite set of coupled integrodifferential equations.A sharp reduction in the dimensionality of the problem is achieved for rational g(ω) and specially for the Lorentzian distribution. 8As the Kuramoto order parameter 2 Z = e iθ is related to α via Z * (t) = ∞ −∞ α(ω, t)g(ω)dω, we can evaluate this integral resorting to the residue theorem obtaining Z * (t) = α(ω 0 − i , t). (This is the result of performing an analytic continuation of α from real to complex ω and evaluating α at the pole of g(ω) in the lower half ω-plane.)Thus, setting ω = ω 0 − i in (9), we get a complex-valued ODE for the Kuramoto order parameter, where B and H have been written in terms of f 1 and f 2 according to Eq. (7).To close Eq.( 10), we need to express the mean field A as a function of Z.For the pulse shape in Eq. ( 4) and a Lorentzian frequency distribution, it can be proven (see Ref. 6 or the supplementary material of Ref. 4) that Note that 0 ≤ A ≤ A max , where the maximal value all oscillators exactly at θ j = 0).In addition to this, the central natural frequency ω 0 is hereafter set to 1, as this can always be achieved by rescaling time and f 1,2 .

IV. FOUR ILLUSTRATIVE PRCs
Among the infinite set of functions f 1 (A) and f 2 (A), we selected a few illustrative case studies.In each of these cases, the character of the PRC undergoes a crossover as A grows: from one iPRC  (12)].The code in the last column indicates the iPRC and the asymptotic PRC at large A (see text).

Case
type to a different PRC type for large A. We denote the limiting PRC at A → ∞ as "asymptotic PRC" (aPRC).From now on, we apply the classical distinction between types I and II to both iPRCs and aPRCs.Recall that if the sign is the same for all θ , we call the iPRC (or the aPRC) as type I (implying f 2 = 0), while in the complementary case with f 1 = 0, we refer to the iPRC (or to the aPRC) as canonical type II or simply type II.Notably, type II may either promote or impede synchronization depending on the sign of f 2 .In turn, we distinguish between two subclasses of type II: II s ( f 2 > 0) and II r ( f 2 < 0) corresponding to the synchronizing and repulsive interactions, respectively.As we are interested in introducing one crossover in the PRC between the iPRC and the aPRC and have three fundamentally different types (I, II s , and II r ), this gives six possible combinations.However, we shall consider only four of these iPRC-aPRC pairs, since only type II s favors synchrony and is to be included either in the iPRC or in the aPRC.Otherwise, no synchronization phenomena are expected: type I is neutral and type II r is repulsive.Hence, we focus on the four cases listed in Table I, in which different PRC types characterize small and large A regimes.As a guide, in the fourth column of the table, we write a code X-Y, where X refers to the iPRC and Y to the aPRC.The saturation function σ (A) in the table has positive slope at A = 0 and saturates at large A. In particular, we chose this specific saturation function in our study, In each panel, the PRC appears divided by A as usual, 13 and the lack of overlapping between different lines evidences its nonlinearity.In Sec.V, we obtain the phase diagrams corresponding to each of the four cases introduced here, based on the analysis of the complex-valued ordinary differential Eq. (10).But before doing so, it is worth making direct simulations of full system (2) and test (and understand) the correspondence with the solutions of Eq. (10).We simulated the full model in case d with ω 0 = 1, heterogeneity parameter = 0.01, pulse-shape parameter r = 0.9, and coupling constant ε = 0.4.As may be seen in Fig. 3, the population exhibits bistability between a desynchronized state and a synchronized state with some oscillators oscillating with the same frequency.This bistability is not surprising as the system is "more synchronizing" when already synchronized since the aPRC is of type II s , while it is hardly  I).The code iPRC-aPRC is indicated in each panel.
synchronizable when already desynchronized by virtue of the type II r iPRC.In terms of Z, the synchronous solution is (approximately) a periodic orbit, while the desynchronized state exhibits only small fluctuations around a point due to finite size effects (N = 500).The agreement with the stable fixed point and the stable limit cycle of Eq. ( 10), also represented in Figs.3(b) and 3(d), is excellent.

V. PHASE DIAGRAMS
In the remainder of this paper, we obtain the phase diagrams for the four reference cases by means of Eq. (10).As (10) is a (generic) planar system, the only possible attractors are fixed points and limit cycles.Their bifurcation loci, depicted in the phase diagrams below, have been obtained using the MATCONT toolbox 36 of MATLAB.Moreover, we recall that Z is only physically meaningful inside the unit disk |Z| ≤ 1 and, therefore, attractors and bifurcations occurring outside it are ignored.As seen in Fig. 3, limit cycles correspond to synchronized solutions in which a macroscopic part of the population rotates at the same average frequency.

A. Case a: I − II s
In Fig. 4(a For small coupling (and heterogeneity), synchronization emerges from a supercritical Hopf bifurcation undergone by the desynchronized state, akin to the classical Kuramoto transition. 8his Hopf bifurcation line terminates at a double zero eigenvalue (Bogdanov-Takens, BT) point.A homoclinic (Hom) line emanates from the BT point limiting the coexistence region.As observed for the regular Winfree model, 5,6 synchronization is more efficient for narrow pulses.The pulse width does not qualitatively change the phase diagram.
The phase diagram only differs appreciably from those in Ref. 6 at the origin.We see that, due to the type I iPRC, the Hopf line approaches the origin with an infinite slope.In particular, the asymptotic dependence of the critical ε H on follows an unusual square-root law with the frequency dispersion , 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 .

B. Case b: II s − I
In case b, iPRC and aPRC are interchanged with respect to case a.This means that synchronization is favored at small coupling but becomes increasingly difficult as the coupling grows.Accordingly, At large ε, there is a bistability region such that the synchronized state disappears in a saddle-node bifurcation of limit cycles (SNLCs).The locus of the SNLC is a line that emanates from a generalized Hopf (or Bautin) point (GH) and terminates at the ε axis at a point marked with a star on the ε-axis of the phase diagram.The stars pinpoint the (equivariant) transcritical (TC) bifurcation, 37 in which the fully synchronized state [θ i (t) = θ j (t)] of identical oscillators ( = 0) becomes unstable.For r = 0.9, the instability of full synchronization takes place at ε c = 9.555 . .., far above the range of ε displayed in the phase diagram.The location of ε c was not calculated using (10) but by directly looking for the stability threshold to the fully synchronized state, see the Appendix.
Finally, note that the synchronization region shrinks as the pulse becomes wider, but there is not a qualitative change in the phase diagram whatsoever.

C. Case c: II s − II r
In this case, the aPRC is repulsive, in contrast to case b, where the aPRC is type I (i.e., neutral in terms of synchronization).In turn, the phase diagram in Fig. 4(c) shows a quite small synchronization region (notice the scale of the axes).Synchronization is bounded exclusively by a supercritical Hopf bifurcation, save for broad pulses.
In the latter case, a GH point is found, and the Hopf bifurcation is subcritical at the left of it.Accordingly, we find a bistability region bounded by a line of saddle-node bifurcation of limit cycles (SNLCs) and a subcritical Hopf bifurcation, as in case b.The precise value of r below which the bistability region exists (i.e., the GH point is present) is r * 0.278 91.
Note also the presence of a TC point in the phase diagram at = 0, above which full synchrony destabilizes. 38The transcritical bifurcation is not structurally stable, see, e.g., Fig. 1 in Ref. 39  and increasing from 0 may either leave no trace of bifurcation or "decay" into two saddle-node bifurcation of limit cycles.The latter scenario occurs for r < 0.275 77 . .., see the bifurcation lines for r = 0.1 in Fig. 4(c), but in our case, one of the bifurcations is not shown as it entails |Z| > 1.

D. Case d: II r − II s
Case d exhibits the most complex phase diagram among all those obtained here.The aPRC is of type II s , as in case a, and (accordingly) the large region is organized by two codimensiontwo points: the Bogdanov-Takens (BT) and the saddle-node separatrix-loop (SNSL) codimension-two points.The associated region of bistability between synchrony and asynchrony is bounded by homoclinic, saddle-node, and Hopf bifurcations.
Remarkably, there is also a bistability region at small ε values for r > r * 0.278 91 (recall the simulations in Fig. 3), which is bounded by a subcritical Hopf and a saddle-node of limit cycles bifurcations.In contrast to previous cases, this synchronization region is detached from the origin due to the repulsive character of the iPRC.To be more precise, the bottom corner of the lower bistability region located at point TC approaches the origin as r → 1.

VI. CONCLUSIONS
In this work, we have studied a non-trivial extension of the Winfree model in which the PRC is nonlinear in the mean field.If the PRC contains only the first harmonic of the angle, the OA ansatz permits a sharp dimensionality reduction.Among all possible dependencies of the PRC on the mean field, we have considered only those with a crossover between two different canonical components.In particular, we have analyzed four cases in which an attractive type II component competes either against a repulsive type II component or against a type I component.Synchronization regions are peculiar for each case.Bistability between macroscopic synchronization and complete desynchronization is found in all cases (in case c, only for broad pulses), but in different relative locations in the -ε plane.
Our results indicate that the nonlinearity of the PRC with the forcing, by itself, is not enough to generate complex collective phenomena.This is certain for a Lorentzian distribution of frequencies since the reduced system is only two dimensional, irrespective of the exact form of f 1 (A) and f 2 (A).As happens in Kuramotolike models, phenomena such as clustering or glassy dynamics may require multiple Fourier components 40 (in the PRC) or stronger heterogeneity, 41 respectively.Concerning collective chaos, other ingredients such as a time-varying coupling, 42 two interacting populations, 43 or multimodal frequency distributions 44 appear to be imperative.
Needless to say, our study is only a drop in the ocean of possible PRCs and model generalizations.For instance, relaxation oscillators 45 and bursting (neuronal) oscillators 46 have PRCs very different from the first-harmonic shape function in Eq. (3).Nevertheless, in spite of its limitations, we regard the model defined by Eqs. ( 2) and (3) as a noteworthy example of a system in which the OA theory can be fully applied.

)(
Our results have been occasionally tested against another choice σ (A) = tanh(A), finding no qualitative difference.)Graphical representations of the four PRCs (cases a-d), for four representative A values, are shown in Figs.2(a)-2(d).

FIG. 2 .
FIG. 2. The non-infinitesimal PRCs analyzed in this work as a function of θ for four representative values of A, including A = 0 + (iPRC) and A = 8 (resembling the aPRC).Panels (a)-(d) correspond to cases a-d, respectively (see TableI).The code iPRC-aPRC is indicated in each panel.
), we show the phase diagram spanned by parameters and ε.Bifurcation lines for three values of parameter r, controlling the pulse width, are depicted.The results almost replicate those in Ref. 6 for the standard Winfree model with type II s iPRC.Synchronization is found in two adjacent regions, in one of them (dark shaded) coexisting with a desynchronized state.[There exists a region (not shown) besides the bistability region where two desynchronized states coexist, see Refs. 5 and 6].In contrast to the averaging approximation (the Kuramoto-Sakaguchi model), valid at small ε and , synchronization becomes impossible if the population is too heterogeneous (large ).

FIG. 3 .
FIG. 3. (a) Raster plot-a dot indicates at which time one oscillator phase crosses a multiple of 2π-for a population of N = 500 and the non-infinitesimal PRC of case d, see Table I and Fig. 2(d).The initial condition is uniform θ j (t = 0) = 0.01, and parameters are = 0.01, r = 0.9, ε = 0.4.The frequencies are deterministically drawn from a Lorentzian distribution: ω i = ω 0 + tan[π(2i − N − 1)/(2N)].(c) The same as (a) but for random initial distribution of phases.(b) and (d) depict the Kuramoto order parameter Z(t) = N −1 j e iθ j (t) for 50 t.u., once the simulations in (a) and (c), respectively, reached the 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.

FIG. 4 .
FIG. 4. Synchronization regions of the Winfree model in the ( , ε)-plane for the four cases of PRCs described in Table I, and three different values of the parameter r ∈ {0.1, 0.5, 0.9}.Panels (a)-(d) correspond to cases a-d, respectively.For the value r = 0.9, light shaded regions indicate where there is a stable limit cycle, corresponding to a macroscopic synchronized state.In the dark shaded regions, the limit cycle (synchronous state) coexists with a stable fixed point (asynchronous state).The codimension-two points are depicted by specific symbols: generalized Hopf (GH-), saddle-node separatrix-loop (SNSL-), Bogdanov-Takens (BT-•), and transcritical (TC-) bifurcations.The bifurcation corresponding to each line type is indicated in the legend of the respective panel.Insets in panels (c) and (d) are magnifications of the regions inside the respective rectangles.

TABLE I .
The four cases of non-infinitesimal PRCs analyzed in this paper.Functions f 1 and f 2 determine the PRC in (3) [σ (A) is a crossover function, see