Low-dimensional dynamics of populations of pulse-coupled oscillators

Large communities of biological oscillators show a prevalent tendency to self-organize in time. This cooperative phenomenon inspired Winfree to formulate a mathematical model that originated the theory of macroscopic synchronization. Despite its fundamental importance, a complete mathematical analysis of the model proposed by Winfree ---consisting of a large population of all-to-all pulse-coupled oscillators--- is still missing. Here we show that the dynamics of the Winfree model evolves into the so-called Ott-Antonsen manifold. This important property allows for an exact description of this high-dimensional system in terms of a few macroscopic variables, and the full investigation of its dynamics. We find that brief pulses are capable of synchronizing heterogeneous ensembles which fail to synchronize with broad pulses, specially for certain phase response curves. Finally, to further illustrate the potential of our results, we investigate the possibility of `chimera' states in populations of identical pulse-coupled oscillators. Chimeras are self-organized states in which the symmetry of a population is broken into a synchronous and an asynchronous part. Here we derive three ordinary differential equations describing two coupled populations, and uncover a variety of chimera states, including a new class with chaotic dynamics.

Synchronization in large systems of oscillators is an ubiquitous phenomenon in nature [1][2][3]. In many cases, synchrony in these natural systems is achieved via brief pulsatile signals emitted by the individual oscillators. Well-known examples of pulsatile interactions producing synchrony are the action potentials emitted by neurons and other cells [1,4], the flashes of light emitted by fireflies [5], or the sound of hands in clapping audiences [6].
Successful theoretical accounts of collective synchronization have typically assumed weak interactions among the oscillators [7,8]. This allows to characterize the state of each oscillator with a single phase variable θ, which is the key idea of the prototypical phase model first proposed by Winfree [1,4,7]. This model considers an ensemble of N oscillators interacting all-to-all with a coupling strength ε: where i = 1, . . . , N , and the number of oscillators is assumed to be very large (N ≫ 1). The presence of heterogeneity in the population is usually modeled via the natural frequencies ω i , which are drawn from a certain probability distribution g(ω) [7,[9][10][11] (but also see [12]). The phase-resetting curve, or phase response curve (PRC) Q, measures the degree of advance or delay of the phases when the oscillators are perturbed. Here the PRC has a sinusoidal shape: A possible choice relating the offset σ and β is σ = sin β. Then the PRC vanishes at θ = 0, as it is naturally assumed in neuronal modeling. If β < π/2 neuronal oscillators are referred to as Type-II, whereas β = π/2 corresponds to a Type-I neuronal oscillator [12][13][14][15][16][17].
We complete the definition of system (1) with the pulse-like coupling: where the integer parameter n ≥ 1 allows to control the width of the pulses. The normalizing constant a n is chosen so that the integral of P (θ) equals 2π. Thus a 1 = 1, and for other values of n: a n = 2 n (n!) 2 /(2n)!. Note also that the n → ∞ limit of (3) is P (θ) = 2πδ(θ).
In this Letter we show that, as Kuramoto-like models, systems of the form (1) evolve into the so-called Ott-Antonsen (OA) manifold [20]. This permits us to describe the highdimensional model (1) by a system of two ODEs -for a specific distribution g(ω). Moreover, this important property obeyed by model (1) permits us to uncover the existence of a large variety of the so-called chimera states [30] in simple networks of pulse-coupled oscillators.
Model (1), with the PRC in Eq. (2) belongs to a family of models that can be rewritten as: In our case, B(t) = εσh(t) and H(t) = εe −iβ h(t), with the mean field h(t) = N −1 j P (θ j (t)). In the thermodynamic limit N → ∞, systems of type (4) have solutions in a reduced invariant manifold discovered by Ott and Antonsen in 2008 [20], which corresponds to a uniform distribution of certain constants of motion at each value of ω [25]. For B = 0, it has been proven [21,22] that, provided the ω's are drawn from a probability distribution function g(ω) which is differentiable and well-behaved in a certain way (see [22] for details), like Lorentzian or Gaussian functions, the dynamics of (4) converges to the OA manifold. One may verify that the proof in [21,22] also holds for B = 0. Hence, we use the OA ansatz with the certainty that it captures the asymptotic dynamics of Eq. (1).
Let F (θ, ω, t) dθ be the ratio of oscillators with phases between θ and θ + dθ and natural frequency ω at time t. The OA ansatz reads: where c.c. stands for complex conjugate. The dynamics of F is governed by the continuity equation since the number of oscillators is conserved. Subtituting the ansatz (5) in Eq. (6), we find that α(ω, t) necessarily obeys This is still an infinite set of equations, but it allows the study of basic problems with pulse-coupled oscillators that remain unsolved so far. As a simple testbed case, we assume a Lorentzian distribution To put our results in the proper context we note that applying the classical method of averaging [8] -valid at very small values of both ε and frequency dispersion-we can reduce the model defined by Eqs. (1), (2) and (3) to the Kuramoto-Sakaguchi model [19,29]: The effective coupling is ε eff = nε/(n + 1), and therefore, the larger n the stronger is the effective coupling. For the Lorentzian distribution (8), an analytic solution for the coupling at the emergence of the synchronous state exists [19]: Recall that this linear dependence of ε c on ∆ is an approximation, and as shown next, within the OA theory a complete description of the population and its transitions is possible. For the analysis that follows, it is convenient to write the generalized order parameters [23]: with m ∈ N. Recalling the ansatz (5), and noting that α admits an analytical continuation into the lower-half complex ω-plane [20], we can evaluate (11) applying the residue theorem. Since the Lorentzian function (8) has one simple pole ω p = ω 0 − i∆ inside the contour, we obtain that all order parameters depend on the value of α at the pole The dynamics of the Kuramoto order parameter Z 1 ≡ Re iΨ is governed by two ODEs obtained equating ω = ω p in (7): In order to make Eq. (12) meaningful for a particular interaction function P (θ), we express the mean field h in terms of R and Ψ. For n = 1 the result is trivial: h 1 = 1 + R cos Ψ. For n > 1, with the important observation that the generalized order parameters are powers of the Kuramoto order parameter Z m = Z m 1 [25], we obtain after some algebra: It is remarkable that, independently of the value of n, the dynamics is always governed by the two ODEs in Eq. (12). This equation cannot be solved analytically but the loci of the bifurcations, where the qualitative behavior changes, can be easily found by standard numerical continuation techniques.
We rescale ∆, ε, and time by ω 0 , so that ω 0 = 1 hereafter. Additionally, we select σ = sin β -see the insets in Fig. 2. We have observed that the results are qualitatively the same independently of n and β, hence the phase diagram for β = 0 and n = 10 in Fig. 1 accounts for all the phenomenology of the model. The case β = 0 was already studied in [9] for a uniform distribution g(ω) obtaining a similar result, albeit some essential differences show up due to the different support of the distributions. In Fig. 1 a Hopf bifurcation line emanates from the origin -with the slope predicted by Eq. (10)-limiting the region of synchronization together with the other solid lines. In the synchronous state a macroscopic cluster of oscillators rotates with the same coupling- modified frequency Ω, and as a result, the order parameter and the mean field oscillate; see the lateral panels of Fig. 1. Notice that, in addition, a cluster of oscillators with Ω = 0 and quivering near θ = 0 is present for all ε > 0. The region of synchronization is bounded at large values of ε by a homoclinic (hom) and a saddle-node on the invariant cycle (SNIC) bifurcations. The latter bifurcation line intercepts the ε-axis at (n + 1) n+1 /[a n (2n + 1) n+1/2 ], i.e. ε = 0.6735 . . . for n = 10. In the main panel of Fig. 1 we see that the Hopf line ends at a Takens-Bogdanov (TB) point [31], which with two other (codimension-two) points organize the region where Hopf and SNIC bifurcations meet; and this conveys bistability between the Synchronous and the Asynchronous states inside a small region bounded by the dashed (saddle-node bifurcation), Hopf, and homoclinic lines. After the preliminary introduction to the model dynamics, we focus on the influence of the width of the pulses P on synchronization. The boundaries in Fig. 2(a) for n = 1 and 10 evidence that the region of synchronization enlarges as n grows, as suggested by Eq. (10). It is interesting to note that the coordinate ε of the TB point diverges with n, while the ε values of the SNIC line decrease. As a result the region of bistability widens as n grows since the SNIC bifurcation at the ε-axis approaches the finite value e 2π = 0.6577 . . . as n → ∞, while the ε coordinate of the TB point progressively grows. The study of large n values is difficult due to the highly nonlinear character of function h, see Eq. (13). It is therefore useful from a mathematical perspective to consider the idealization P (θ) = 2πδ(θ). Using the trigonometric representation of the Dirac's delta function we obtain the mean field: Note that in this derivation the n → ∞ limit is taken after the N → ∞ limit, and therefore any subsequent result using h ∞ is expected to be a truly asymptotic one as n grows provided N is kept sufficiently large. On the contrary, implementing instantaneous interactions (n = ∞) with a finite population (N < ∞) cannot fit in the theory since the mentioned limits do not commute. Inserting h ∞ in Eq. (12) we obtain the boundaries [32] shown with dashed lines in Figs. 2(a) and 2(b). We see that for β = 0 there is already a noticeable similarity between the regions of synchronization for n = 10 and n = ∞. The main discrepancy is observed at high ε values, which is not particularly interesting since, in any case, almost the whole population does not rotate in that region (see [33] for a description of this effect). As said above, as n grows the TB point moves upwards, so that the synchronization regions eventually match at the n → ∞ limit (note nevertheless that the limit is somewhat singular because the bistability region disappears).
For β = 1, see Fig. 2(b), the difference between the results for n = 10 and n = ∞ becomes apparent, and more tangible than what could be naively expected from Eq. (10). In fact, the closer β approaches to π/2, the more favorable is a sharp P (θ) to achieve synchronization. We claim this is a general statement, since we have also observed it numerically with Gaussian g(ω). It may be conjectured that the effectiveness of sharp spikes to achieve synchronization is one reason for their common occurrence in nature.
Finally, as an illustration of the potential of the OA ansatz in analyzing ensembles of pulse-coupled oscillators, we focus next on the emergence of chimera states. In a chimera state the system splits into a fully synchronized subpopulation and a desynchronized one. The surprising feature of this state is that full synchronization of all oscillators is a stable solution and, therefore, the chimera state does not appear via a usual symmetry breaking mechanism. The simplest set-up capable of sustaining chimeras, in both experimental [34,35] and numerical [26,27] realizations, consists of two coupled subpopulations b = (1, 2) of identical oscillators (ω i = 1). Here, in contrast to previous theoretical works [24][25][26][27][28] with the Kuramoto-Sakaguchi coupling, we consider pulse-coupled oscillators: with b ′ = (2, 1). We can recognize in the equation for each subpopulation the structure of Eq. (4): In consequence, there is a solution in which each subsystem evolves in its own OA manifold (5). The absence of diversity in the populations makes the OA manifold to be neutrally stable and to be surrounded by a foliation of infinite neutrally stable solutions [24,25], each corresponding to different values of the Watanabe-Strogatz invariants [36] (introduced by the initial conditions). Nevertheless, the OA manifold becomes attracting as soon as a tiny amount of diversity is present [28]. Thus, in some sense, the OA manifold is the "skeleton" of the phase space, and it is legitimate to analyze the system with the OA ansatz.
It is straightforward to write the ODEs governing the dynamics of the order parameter of the b-th subpopulation As we are interested in chimera states where one subpopulation is fully synchronized, say the first one (R 1 = 1), the equation for R 1 disappears. We get then a system of only three ODEs (for Ψ 1 , R 2 , and Ψ 2 ) that makes possible to carry out an exhaustive exploration of the chimera states.
It is convenient for the analysis to define two parameters: A = (µ − ν)/(µ + ν) quantifying the unbalance between intra-and inter-population interactions, and S = µ + ν quantifying the coupling strength. Interestingly, in the limit of |µ|, |ν| ≪ 1, irrespective of the values of σ and n, the system reduces (via averaging) to two ODEs for R 2 and ψ = Ψ 1 −Ψ 2 identical to those in Eq. (12) of Abrams et al. [27] for oscillators of Kuramoto-Sakaguchi type. Hence, for S → 0 we can borrow the results in [27], in particular, the existence of chimeras only for β values not far from π/2. When S is not small the system behaves as genuinely threedimensional. This increase of dimensionality has a deep impact on the types of chimera states observed; see the rich phase diagram in Fig. 3 (obtained resorting to numerical continuation of the bifurcations of limit cycles [31] and computing the Lyapunov exponents [38]). The phase diagram is qualitatively generic since we have obtained similar results in a wide range of the parameters σ, S, and n. In fact, the structure the phase diagram is reminiscent of the one in Fig. 4 of Ref. [27], but now a much richer scenario emerges due to the additional degree of freedom. Above the Neimark-Sacker (or secondary Hopf) bifurcation, signaled by a thick (blue) line in Fig. 3, we find quasiperiodic chimeras, and the expected resonance tongues corresponding to limit cycles on the surface of the invariant torus. As we move away from the Neimark-Sacker bifurcation the torus breaks down [37] and the resonances merge giving rise to an intricate set of bifurcations (not shown, see [39]). Perhaps the most remarkable consequence of the torus break-down is the existence of chaotic chimera states in the dark (red) shaded region of Fig. 3. The upper left panel shows, for specific values of A and β, the projection of the chaotic trajectory onto the R 2 e iΨ2 plane.
In sum, this work shows that populations of pulse-coupled oscillators evolve in the Ott-Antonsen manifold. The main limitations for using the OA ansatz are the uniform all-to-all coupling in (1), and the sinusoidal shape of the PRC, Eq. (2). This leaves room for a number of possible extensions of our results that may be of relevance in neuroscience, among other fields.
We thank Juan M. López for a critical reading of the manuscript, and Arkady Pikovsky for interesting discussions. DP acknowledges support by Cantabria International Campus, and by MINECO (Spain) through a Ramón y Cajal fellowship and project No. FIS2009-12964-C05-05.