Synchronization scenarios in the Winfree model of coupled oscillators

TL;DR

This paper derives an exact low-dimensional model for the Winfree oscillator system, revealing two distinct synchronization scenarios linked by a bifurcation, and compares these with Kuramoto-Sakaguchi predictions.

Contribution

It provides an exact low-dimensional representation of the Winfree model and identifies two different synchronization scenarios based on pulse shape and PRC offset.

Findings

Two distinct synchronization scenarios identified

A bifurcation links the scenarios via a Bogdanov-Takens point

Discrepancy with Kuramoto-Sakaguchi prediction behaves as an odd function of PRC offset

Abstract

The emergence of collective synchronization was reproduced long ago by Winfree in a classical model consisting of an ensemble of pulse-coupled phase oscillators. By means of the Ott-Antonsen ansatz, we derive an exact low-dimensional representation which is exhaustively investigated for a variety of pulse types and phase response curves (PRCs). Two structurally different synchronization scenarios are found, which are linked via the mutation of a Bogdanov-Takens point. From our results, we infer a general rule of thumb relating pulse shape and PRC offset with each scenario. Finally, we compare the exact synchronization threshold with the prediction of the averaging approximation given by the Kuramoto-Sakaguchi model. At the leading order, the discrepancy appears to behave as an odd function of the PRC offset.

Tables

Table 1. Table 1: Summary of the various pulse functions considered. The normalizing constant of the AS pulse is a n = 2 n ​ ( n ! ) 2 / ( 2 ​ n ) ! = n ! / ( 2 ​ n − 1 ) !! subscript 𝑎 𝑛 superscript 2 𝑛 superscript 𝑛 2 2 𝑛 𝑛 double-factorial 2 𝑛 1 a_{n}=2^{n}(n!)^{2}/(2n)!=n!/(2n-1)!! . The fourth column shows the mean field h ​ ( Z ) ℎ 𝑍 h(Z) , which is the function entering in Eq. ( 11 ) describing the system’s mean field dynamics exactly. In the last column, the shape factor Π Π \Pi quantifies the effective strength of each pulse under the averaging approximation, see Eq. ( 16 ). In Fig.  1 , lines of the same style share the same Π Π \Pi value.

Pulse name	$P\left(\theta \right)$	Parameter	Mean field: $h\left(Z\right)$	Shape factor: $\Pi$
Ariaratnam-Strogatz (AS)	$a_n{\left(1+\cos \theta \right)}^n$	$n\in {\mathbb{Z}}^{+}$	$1+{\left(n!\right)}^2\sum _{k=1}^n\frac{Z^k+{\left(Z^{\ast }\right)}^k}{\left(n+k\right)!\left(n-k\right)!}$	$\frac{n}{n+1}$
Rectified-Poisson (RP)	$\frac{\left(1-r\right)\left(1+\cos \theta \right)}{1-2r\cos \theta +r^2}$	$r\in (-1,1)$	$\text{Re}\left[\frac{1+Z}{1-rZ}\right]$	$\frac{1+r}{2}$
Square	$\{\begin{array}{cc}\pi /b\text{for}\left|\theta \right|\le b& \\ 0\text{otherwise}& \end{array}$	$b\in (0,\pi )$	$1-\frac{1}{b}\text{Im}[\ln (1-Z{e}^{ib})-\ln (1-Z{e}^{-ib})]]$	$\frac{\sin b}{b}$
Dirac delta	$2\pi \delta \left(\theta \right)$	$-$	$\text{Re}\left[\frac{1+Z}{1-Z}\right]$	$1$

Equations

$${\dot{\theta }}_i={\omega }_i+Q({\theta }_i)\frac{\epsilon }{N}\sum _{j=1}^{N}P({\theta }_j).$$

$$Q(\theta )=\frac{\sin \beta -\sin (\theta +\beta )}{\cos \beta }=q(1-\cos \theta )-\sin \theta ,$$

$$g(\omega )=\frac{\Delta /\pi }{(\omega -{\omega }_0{)}^2+{\Delta }^2}.$$

$${\partial }_{t}F+{\partial }_{\theta }\left(F\dot{\theta }\right)=0$$

$$F(\theta |\omega ,q,t)=\frac{1}{2\pi }\left\{1+\left[\sum _{m=1}^{\infty }{\alpha }_m(\omega ,t)e^{im\theta }+\text{mbox}c.c.\right]\right\},$$

$${\partial }_t\alpha =-i(\omega +\epsilon hq)\alpha +\frac{\epsilon h}{2}\left[(1+iq)-(1-iq){\alpha }^2\right].$$

$$h(t)={\int }_{-\infty }^{\infty }g(\omega ){\int }_0^{2\pi }F(\theta |\omega ,t)P(\theta )d\theta d\omega$$

$$Z(t)={\int }_{-\infty }^{\infty }g(\omega ){\int }_0^{2\pi }F(\theta ,\omega ,t)e^{i\theta }d\theta d\omega .$$

$$Z^{\ast }(t)={\int }_{-\infty }^{\infty }g(\omega )\alpha (\omega ,t)d\omega .$$

$$Z^{\ast }(t)=\alpha ({\omega }_p,t),$$

$$\dot{Z}=(-\Delta +i{\omega }_0)Z+\frac{\epsilon h}{2}\left[1-Z^2-iq(1-Z{)}^2\right].$$

$$P(\theta )=\sum _{m=-\infty }^{\infty }{c}_m{e}^{im\theta }$$

$$h=1+\sum _{k=1}^{\infty }{c}_k{\int }_{-\infty }^{\infty }g(\omega )\left\{{\left[\alpha (\omega ,t{)}^{\ast }\right]}^k+{\left[\alpha (\omega ,t)\right]}^k\right\}d\omega ,$$

$$h=1+\sum _{k=1}^{\infty }{c}_k\left[Z^k+{\left(Z^{\ast }\right)}^k\right].$$

$$\zeta =\overline{|Z-\bar{Z}|},$$

$$P(\pi \pm \delta \theta )\simeq \frac{1-r}{2(1+r{)}^2}\delta {\theta }^2,\frac{(n!{)}^2}{(2n)!}\delta {\theta }^{2n},0,$$

$${\Delta }_H(y)$$

$${\epsilon }_H(y)$$

$$g(y,q)$$

$$f(y,q)$$

$${\dot{\theta }}_i={\omega }_i+\epsilon q+\Pi \frac{\epsilon }{N}\sum _{j=1}^N\left[\sin ({\theta }_j-{\theta }_i)-q\cos ({\theta }_j-{\theta }_i)\right]$$

$${\epsilon }_c^{(\text{av})}=\frac{2\Delta }{\Pi },$$

$$\rho (\Delta )=\frac{{\epsilon }_H-{\epsilon }_c^{(\text{av})}}{{\epsilon }_c^{(\text{av})}}$$

$$\rho (\Delta )=\xi (q)\Delta +O({\Delta }^2)$$

Full text

I Introduction

Macroscopic synchronization is a well-known emergent phenomenon arising in ensembles of oscillators when, despite their unavoidable differences, some fraction of the oscillators spontaneously lock to one another and oscillate together with exactly the same frequency Winfree (2002); Strogatz (2000); Pikovsky et al. (2001). Examples of collective synchronization are abundant and surprisingly diverse, see e.g. Strogatz (2003). They include the synchronous flashing of fireflies Buck and Buck (1976), circadian Liu et al. (1997) and cardiac Glass and Mackey (1988) rhythms, the spontaneous transitions to synchronous stepping Strogatz et al. (2005) and to synchronous clapping Néda et al. (2000), or the collective synchronization of chemical oscillators Kiss et al. (2002), and arrays of optomechanical cells Heinrich et al. (2011), and Josephson junctions Wiesenfeld et al. (1996).

The first successful attempt to model macroscopic synchronization is due to Arthur Winfree. In 1967, Winfree proposed a mathematical model consisting of a large population of globally coupled oscillators. Assuming weak coupling, Winfree postulated the dynamics of the individual oscillators to be well described by a single phase variable. Interactions are modeled by means of pulses that are emitted by each oscillator and perturb the phase velocity of all the other oscillators. Mathematically, this is expressed through two independent functions: The (infinitesimal) Phase Response Curve (PRC), determining how the phase of an oscillator changes under perturbations; and a function specifying for the shape of the pulses. Numerical simulations in  Winfree (1967, 1980) showed that, under suitable conditions, the Winfree model displayed a transition from a totally disordered state to collective synchronization, analogously to phase transition in statistical mechanics. Though the Winfree model was later investigated in a few more papers Ariaratnam and Strogatz (2001); Quinn et al. (2007); Basnarkov and Urumov (2009), the interest soon turned to the simpler and renowned Kuramoto model  Strogatz (2000); Pikovsky et al. (2001); Strogatz (2003); Acebrón et al. (2005).

A new boost in the theoretical understanding of phase-oscillator populations models occurred in 2008, when Ott and Antonsen (OA) discovered an exact dimensionality reduction of the Kuramoto model, called OA ansatz Ott and Antonsen (2008, 2009); Ott et al. (2011). The discovery of the OA ansatz opened up the possibility of tackling unresolved problems and investigate novel variants and extensions of the Kuramoto model, see e.g. Antonsen Jr. et al. (2008); Pikovsky and Rosenblum (2008); Childs and Strogatz (2008); Abrams et al. (2008); Abdulrehem and Ott (2009); Pazó and Montbrió (2009); Martens et al. (2009); Lee et al. (2009); Lafuerza et al. (2010); Alonso and Mindlin (2011); Hong and Strogatz (2011); Montbrió and Pazó (2011); Barlev and Ott (2011); Omel’chenko and Wolfrum (2012); Iatsenko et al. (2013); Restrepo and Ott (2014); Pietras et al. (2016); Ott and AntonsenJr. (2017). Remarkably, the OA ansatz is also applicable to pulse-coupled oscillators  Luke et al. (2013); So et al. (2014); Laing (2014, 2015); Montbrió et al. (2015); Pazó and Montbrió (2016); O’Keeffe and Strogatz (2016); Chandra et al. (2017) and, in particular, to the original Winfree model Pazó and Montbrió (2014). This allows to investigate synchronization phenomena which are not accessible using Kuramoto-like models. Specifically, the advantage of the Winfree model is that permits to investigate separately how pulse shape and PRC type influence collective synchronization. Note that the PRC of cells, such as neurons Preyer and Butera (2005); Tateno and Robinson (2007) and cardiac cells Kralemann et al. (2013), can be measured experimentally. Therefore understanding better the Winfree model should contribute to narrow the gap between mathematical models and biological phenomena.

Here we build on our previous work Pazó and Montbrió (2014), and systematically analyze the impact of (i) pulse shapes and (ii) PRC offsets, onto collective synchronization in the Winfree model. We find that the phase diagram obtained in Pazó and Montbrió (2014) is not unique, and that a novel synchronization scenario emerges for certain pulse types via the mutation of a codimension-two Bogdanov-Takens (BT) point. We end investigating the limit in which the oscillators are nearly identical and very weakly coupled. In that situation the averaging approximation is valid, and a Kuramoto-like model captures the dynamics with a level of accuracy that is measured numerically.

The paper is organized as follows: In Sec. II we introduce the Winfree model, and in Sec. III an exact reduction to two ordinary differential equations (ODEs) is derived using the OA ansatz. In Sec. IV we present the results obtained from those two ODEs, for a variety of pulse shapes and PRCs. In Sec. V we compare the Winfree model with its averaging approximation. Finally, in Sec. VI we address the conclusions of this work.

II The Winfree Model

The Winfree model Winfree (1967, 1980) consists of an ensemble of $N\gg 1$ globally coupled phase oscillators with heterogeneous natural frequencies $\omega_{i}$, $i=1,\ldots,N$. The phases $\theta_{i}$ are governed by the set of $N$ ordinary differential equations (ODEs)

[TABLE]

Here, each oscillator receives the input of the mean field $h=N^{-1}\sum_{j=1}^{N}P(\theta_{j}),$ and its response to it depends on its own phase through the PRC function $Q(\theta)$. Both $P$ and $Q$ are $2\pi$-periodic functions on the real line, and hence can be defined either in the range $[0,2\pi)$ or in the range $[-\pi,\pi)$. Finally, the global coupling strength is controlled by the parameter $\varepsilon>0$.

II.1 Pulse shape, $P(\theta)$

The function $P$ in Eq. (1) specifies the form of the pulses. We only consider pulses with the following properties:

1. (i)

$P$ is unimodal and symmetric around $\theta=0$. 2. (ii)

$P$ vanishes at $\theta=\pi$. 3. (iii)

$P$ has a normalized area: $\int_{-\pi}^{\pi}P(\theta)d\theta=2\pi$.

We consider the three pulse types with finite width shown in Fig. 1, and defined in Table 1. The first pulse, labeled as AS, was originally adopted by Ariaratnam and Strogatz Ariaratnam and Strogatz (2001), and is commonly used in recent studies of pulse coupled-phase oscillators Goel and Ermentrout (2002); Pazó and Montbrió (2014); Luke et al. (2013); So et al. (2014); Laing (2014, 2015); O’Keeffe and Strogatz (2016); Chandra et al. (2017). Additionally, we consider a variant of the pulse used by O’Keeffe and Strogatz in O’Keeffe and Strogatz (2016) equal to the Poisson kernel, but with an offset so that it fulfills the condition (ii). We term this pulse as “Rectified Poisson kernel” (RP). Finally we consider a square pulse with a flat profile and vanishing in a finite interval of theta: $[-\pi,-b)\cup(b,\pi)$.

Concerning the macroscopic dynamics of the Winfree model, the precise value of $N$ becomes irrelevant provided it is large enough (i.e. only trivial finite-size fluctuations are observed). However, for Dirac delta pulses, this is not the case, as we discuss in Sec. IV.4. The Dirac delta is the limiting case of the pulse types considered, i.e. $n\to\infty$, $r\to 1$, and $b\to 0$ for the AS, RP, and square pulses, respectively.

II.2 Phase-response curve (PRC), $Q(\theta)$

The influence of a certain (small) perturbation on the phase of an oscillator is determined by the PRC, $Q(\theta)$. Here we assume that (i) the PRC vanishes at the phase where the pulses peak, i.e. at $Q(\theta=0)=0$; and (ii) the PRC has a sinusoidal shape. This latter condition is crucial, for the OA theory to be applicable. The constraints (i) and (ii) lead us to the following one-parameter family of PRCs

[TABLE]

where parameter $q=\tan\beta$ determines the degree of asymmetry of the PRC. As illustrated in Figure 2 , $Q$ is more positive (advancing) than negative for $q>0$, while it is more negative (retarding) for $q<0$. The case $q=0$ corresponds to a perfectly balanced PRC. Hence, we call $q$ ‘offset parameter’ hereafter. Note that in Pazó and Montbrió (2014) the PRC is defined in a slightly different manner: Here $\varepsilon$ is equivalent to $\varepsilon\cos\beta$ in our previous work.

II.3 Frequency distribution, $g(\omega)$

As indicated above, heterogeneity in the population enters through the set of natural frequencies $\omega_{i}$. As we show in the next section, to simplify the analysis of the Winfree model (1) it is convenient to adopt a Lorentzian distribution centered at $\omega_{0}$ with half-width $\Delta$:

[TABLE]

III Low-dimensional dynamics of the Winfree model

In the following we present a reduction of the dimensionality of our problem using the so-called Ott-Antonsen theory, what permits to determine the system dynamics exactly in the thermodynamic limit $N\to\infty$ for all parameter values. Hence, we introduce the density $F$ function, such that $F(\theta|\omega,t)\,d\theta$ represents the fraction of oscillators with phases between $\theta$ and $\theta+d\theta$, and natural frequency $\omega$ at a time $t$. The conservation of the number of oscillators imposes $F$ to obey the continuity equation:

[TABLE]

To solve this equation it is natural to use a Fourier expansion of $F$:

[TABLE]

where c.c. stands for complex conjugate. We adopt the OA ansatz assuming that the $m$-th mode is the $m$ power of the first mode: $\alpha_{m}(\omega,t)=[\alpha(\omega,t)]^{m}$. This drastic reduction of dimensionality was justified in Pazó and Montbrió (2014) following Ott and Antonsen (2009); Ott et al. (2011), see also Pietras and Daffertshofer (2016). Now, inserting (5) into the continuity equation (4) we get an equation for $\alpha(\omega,t)$:

[TABLE]

In this equation the mean field

[TABLE]

couples every $\alpha(\omega,t)$ with all others $\alpha(\omega^{\prime},t)$.

We use the Kuramoto order parameter $Z$ to monitor the macroscopic dynamics of the system. It quantifies the amplitude of first Fourier mode of the density $F$, and reads

[TABLE]

Under the assumption that the system evolves in the OA manifold:

[TABLE]

(The asterisk denotes complex conjugation.) For Lorentzian $g(\omega)$ this integral over the real line can be computed by performing an analytical continuation of $\alpha(\omega,t)$ from real $\omega$ into complex $\omega=\omega_{r}+i\omega_{i}$, see Ott and Antonsen (2008) for details. Closing the integral by a half-circle at infinity in the lower complex $\omega$ half-plane, permits to apply the residue’s theorem, obtaining

[TABLE]

where $\omega_{p}=\omega_{0}-i\Delta$ is the simple pole of $g(\omega)$ inside the integration contour.

The exact, low-dimensional form of the Winfree model with Lorentzian frequency distribution (3) and sinusoidal PRC (2), is obtained setting $\omega=\omega_{p}$ in Eq. (6). Then, we obtain the complex-valued ODE

[TABLE]

To use this equation, one needs to express the mean-field as $h(Z)$. With that aim, it is convenient to expand $P$ in Fourier series:

[TABLE]

with $c_{m}=c_{-m}\in\mathbb{R}$ and $c_{0}=1$, because of the properties (i) and (iii) stipulated in Sec. II.1. Inserting Eqs. (12) and (5) into Eq. (7) we get

[TABLE]

which again can be simplified applying the residue’s theorem, and recalling Eq. (10) allows to express the result in terms of $Z$ only

[TABLE]

This relation permits to achieve, after some algebra, the desired expressions of $h(Z)$ for the set of pulse types. They are listed in the fourth column of Table 1.

Note that, compared to the AS pulse used in previous studies Ariaratnam and Strogatz (2001); Goel and Ermentrout (2002); Pazó and Montbrió (2014); Luke et al. (2013); So et al. (2014); Laing (2014, 2015); O’Keeffe and Strogatz (2016); Chandra et al. (2017), the RP pulse has the advantage that $h(Z)$ remains a simple function of $Z$, no matter how much the pulse width is decreased, see Table I. Moreover, though the mean field function $h(Z)$ for square pulses is more cumbersome than that of the RP pulses, it still permits to investigate the Winfree model with pulses of arbitrary small width, without the drawback of dealing with the long sums of the AS pulse’s mean field.

IV PHASE DIAGRAMS AND PHASE PORTRAITS

Next we exploit the low-dimensional character of Eqs. (11), to fully investigate the bifurcations of the Winfree model for the pulse functions described in Table 1, and various PRC offsets $q$. First note that the dynamics of the Winfree model (11) depends on five parameters: The coupling strength ($\varepsilon$), the pulse width (through $n$, $r$ or $b$), the center and half-width of the frequency distribution ($\omega_{0}$ and $\Delta$) and the PRC offset ($q$). From now on, and without lack of generality, we set $\omega_{0}=1$, since this can always be achieved after a trivial rescaling of time and parameter $\varepsilon$.

IV.1 Rectified-Poisson (RP) pulse

Figure 3 shows the phase diagram for the RP pulse ($r=0.5$) with negative PRC offset ($q=-1$), obtained using Eq. (11) with the assistance of the MatCont toolbox of matlab. The diagram is qualitatively identical to those presented in Ref. Pazó and Montbrió (2014) for the AS pulse, indicating certain robustness of the dynamics against modifications of the pulse shape. In Fig. 4(a) we show a sketch of the phase portraits in the regions of interest. In the shaded region, labeled $2$, Eq. (11) has one attractor of limit-cycle type, meaning that $Z$ exhibits periodic oscillations. This is reflecting a state of macroscopic synchronization in which a certain part of the population is self-entrained to a common frequency. There are three different paths leading to this state, depending on which bifurcation line is crossed: Hopf, SNIC (saddle-node on the invariant circle), or Hom (homoclinic or saddle-loop). Note that the latter one is a global bifurcation that does not destabilize the steady state, and in consequence, a region of bistability between synchrony and asynchrony exists, see the dark shaded region in Figs. 3 and 4(a). Two lines of saddle-node bifurcations of fixed points, emanating from a cusp point (CP), complete the phase diagram and bound a region of bistability between two stable steady states (region 4). They correspond to two asynchronous states with a different number of quiescent oscillators. For large enough $\Delta$, namely above the CP point, the fraction of quiescent oscillators varies smoothly (i.e. non-hysteretically) with $\varepsilon$.

To verify the validity of our analytical results we carried out simulations of the full model with $N=2000$ oscillators. In Figs. 5(a) and 5(b), we present raster plots of the incoherent and the synchronized states, respectively (a dot is plotted every time an oscillator crosses a multiple of $2\pi$). Moreover, to be more systematic we swept parameter $\varepsilon$ along $\Delta=0.15$, i.e. a vertical line in Fig. 3, with the intention of testing that the bifurcations were indeed reproduced. As the rotation of the oscillators is not uniform, $|Z|$ alone is not a good order parameter to detect bifurcations. It is more convenient to use the order parameter proposed by Shinomoto and Kuramoto Shinomoto and Kuramoto (1986):

[TABLE]

where the bar means long-time average. For asynchronous dynamics, the Shinomoto-Kuramoto order parameter Eq. (15) satisfies $\zeta=0$, while $\zeta\neq 0$ indicates some degree of synchronization. Figure 5(c) shows that the results of our numerical simulations of the original Winfree model show a good agreement with the reduced ODE, Eq. (11).

We next investigate how the synchronization region changes as the pulse width varies. Figures 6(a) and 6(b) show the synchronization boundaries for $r=0$, $0.5$ and $0.95$, and opposite values of the PRC’s offset $q$ ($=\mp 1$). For the sake of clarity other bifurcation lines have been omitted. For positive offsets, we find the same result that we found in Pazó and Montbrió (2014) for AS pulses with highly unbalanced PRCs ($q\gg 0$): Narrow RP pulses ($r$ close to 1) are more efficient than broad pulses to synchronize heterogeneous populations of oscillators. Indeed, note that the synchronization boundary of the narrowest pulse ($r=0.95$) reaches the highest value of the heterogeneity parameter $\Delta$ in Fig. 6(b). However, a small discrepancy with this previous rule was already noticeable in the $q=0$ curve of Fig. 2(a) in Pazó and Montbrió (2014). Here we revisit that question and find that, as Fig. 6(a) shows, for negative PRC offsets the discrepancy is even more dramatic: the Hopf boundary is far from attaining the largest $\Delta$ value for the narrowest pulse. Hence, synchronization is not optimal for narrow pulses, but it also depends on the sign of the PRC’s offset $q$. Consequently, one is tempted to wonder if, in nature, adaptation may in some cases drive PRC offsets and pulse widths to be mutually optimized in a certain sense.

IV.2 Ariaratnam-Strogatz (AS) pulse

Thus far we found no qualitative difference between the phase diagram for RP pulses of Fig. 3, with that of Ref. Pazó and Montbrió (2014), obtained using AS pulses. Nonetheless, in this section we show that this qualitative agreement breaks down for AS pulses with PRCs with negative offset.

The AS pulse with $n=1$ is identical to the RP pulse with $r=0$, so that no differences arise in this case. Surprisingly, when we considered narrower pulses (larger values of $n$) a more complicated bifurcation scenario showed up, see Fig. 7 for $n=5$ and $q=-1$. Indeed, at a certain critical $n$, the Bogdanov-Takens point mutates its character in such a way that the Hopf line emanating from it becomes of subcritical type, while the homoclinic bifurcation now involves an unstable periodic orbit —because the sum of the eigenvalues of the saddle point, called saddle quantity, is positive. This mutated Bogdanov-Takens point is designated as BT’ hereafter. Points BT and BT’ are both equally generic zero-eigenvalue points consistent with the normal form in textbooks Guckenheimer and Holmes (1983); Kuznetsov (1998): BT is the usual representation (up to a transformation of parameters), while BT’ is also consistent upon time inversion.

In the transition from BT to BT’ two new codimension-two points appear:

A generalized Hopf (GH) point on top of the Hopf line where the bifurcation shifts from super- to sub-critical. 2. 2.

A neutral saddle (NS) point where the homoclinic connection is degenerate, since it involves a saddle point with zero saddle quantity Kuznetsov (1998). At the NS point the line of homoclinic bifurcation of the stable limit cycle terminates.

The GH and NS points are connected by a new line curve, which is the locus of a saddle-node bifurcation of limit cycles (SNLC). Figure 4(b) shows sketches of the phase portraits when a BT’ point is present in the phase diagram. Notably, the synchronization region is detached from the BT’ point, and a region with three attractors (i.e., tristability) exists. This region is the approximate triangle with vertices at GH and NS visible both in the inset of Fig. 7 and in Fig. 4(b), region 5. There, the limit cycle (corresponding to synchronization) coexist with two stable fixed points. Note that by entering into region 5 through the saddle-node bifurcation of limit cycles line (connecting the points GH and NS) a finite-sized limit cycle with a finite basin of attraction suddenly appears.

IV.3 Transition between the synchronization scenarios BT and BT’

In view of the distinct phase diagrams associated to BT and BT’, next we investigate the conditions under which each scenario shows up. Our systematic numerical investigation indicates that the RP pulse is always associated to a BT point. In the case of the AS pulse, we determined numerically the threshold value of $q$, which we designated as $q_{*}$, where the transition between BT and BT’ occurs, i.e. the $q$ value at which a degenerate (codimension-three) BT point arises. The result covering all integer values $n\leq 10$ is depicted in Fig. 8, and demonstrates that the BT’ point only arises for sufficiently negative offsets $q$. Noteworthy, when $n$ grows, BT’ can be observed for increasing small values of $|q|$. The absence of a point for $n=1$ in Fig. 8 is not an omission; in fact, we failed to numerically find a BT’ point even after considering extremely small values of $q$ —recall also that the AS pulse with $n=1$ coincides with the RP pulse with $r=0$. Finally, we carried out numerical simulations using square pulses (not shown), and found that the bifurcation scenario associated to BT’ is observed already for $q=0$, provided that $b$ is smaller than $b_{*}(q=0)=1.02\ldots$ (a significantly large value).

To get some more physical insight, we examined the asymptotic behavior of $P(\theta)$ in a neighborhood of $\theta=\pi$ for each pulse type. They are:

[TABLE]

for RP, AS, and square pulses, respectively. The marked differences of the respective asymptotics led us to conjecture a simple rule of the thumb: pulses that fall fast enough to zero at $\theta=\pi$ are prone to exhibit the synchronization scenario with five codimension-two points, i.e. BT’. On the contrary, for pulses that fall to zero more slowly (such as the RP pulse) favor the first scenario (BT), making the second scenario (BT’) impossible or only present for small enough PRC offsets $q$.

IV.4 Dirac delta pulse

All the pulses studied in this paper have the Dirac delta as limiting case. It is not difficult to obtain the expression of $h(Z)$ for the Dirac delta pulse, see Table I. Nonetheless some caution must be taken here: for obtaining the mean field $h(Z)$ the thermodynamic limit ($N\to\infty$) is assumed prior to the zero width pulse limit, and it is well known that these two limits do not commute Zillmer et al. (2007). Therefore, the results we obtain for Dirac delta pulses cannot be exactly reproduced in numerical simulations, which necessarily involve a finite number of oscillators. Accordingly, the results obtained here for Dirac delta pulses must be interpreted as a limit of the bifurcation lines for very narrow pulses. This allows us to put aside the pulse shape, and to focus solely on the influence of the PRC offset parameter $q$.

Figure 9 shows phase diagrams in the $(\Delta,\varepsilon)$-plane for Dirac delta pulses, and for several values of $q$. The curves displayed are Hopf bifurcation lines that emanate from the origin and approach the vertical axis when $\varepsilon\to\infty$. As mentioned above, there are subtle questions regarding this pulse, so the Hopf bifurcation lines have to be understood simply as the limit of the Hopf curves for very narrow pulses. In fact, the absence of the saddle-node bifurcations lines indicates a certainly singular behavior in that limit.

Yet, from Fig. 9, we can conclude that the synchronization region increases monotonically with $q$. Our physical interpretation of this feature is that, for negative PRC offsets, the bias of the PRC tends to slow down the oscillators favoring the formation of a cluster with quiescent oscillators (partial oscillation death). On the contrary, positive PRC offsets generally favor phase advances, retarding the accumulation of quiescent oscillators, and leaving room for the synchronization to occur more widely.

Let us finally point out that the bifurcation lines can be analytically obtained by transforming Eq. (11) into a complex-valued ODE for a new variable $w=(1+Z)/(1+Z)$, such that $h$ equals $\text{Re}(w)$. In the new coordinate system, and with the assistance of mathematica, we derived convoluted, but nonetheless exact equations of the Hopf boundaries in parametric form:

[TABLE]

with

[TABLE]

where $y\in(0,\infty)$. In passing, we note, that for $q=0$ a simple explicit formula can be found, see Ref. [36] in Pazó and Montbrió (2014).

V Limit of weak coupling and nearly identical oscillators

To conclude, we investigate the Winfree model in the limit of weak coupling and weak heterogeneity, i.e. $|\varepsilon|\ll 1$, $\Delta\ll 1$. This is an important limiting case, since the method of averaging can be applied and the Winfree model reduces to the well-studied Kuramoto-Sakaguchi (KS) model Kuramoto (1984); Sakaguchi and Kuramoto (1986). Our aim in this section is to investigate how the synchronization threshold of the Winfree model deviates from that of its corresponding KS model, for different pulse functions and PRC types.

V.1 Averaging approximation: Kuramoto-Sakaguchi model

In the classical analysis by Kuramoto Kuramoto (1984), weakly interacting oscillators with frequencies close to a resonance are described by means of the averaging approximation. In the case of 1:1 resonance (nearly identical frequencies), the interaction term between any two oscillators becomes a function of their phase difference.

Given that the PRCs considered here are chosen to be sine-shaped, the only resonant term is the first harmonic. Thus, the averaging calculation leads to the KS model Kuramoto (1984); Sakaguchi and Kuramoto (1986); Pazó and Montbrió (2014):

[TABLE]

The parameter $\Pi$ is a “shape factor” that depends on the pulse shape. The shape factor depends only on the first harmonic of the pulse, or more precisely, $\Pi=c_{1}$, see Eq. (12). The dependence of $\Pi$ on the parameter controlling the pulse width can be found in the last column of Table 1. In all cases, the effective interaction strength increases as the pulses become narrower. In fact, the largest $\Pi$ value is attained for the Dirac delta pulse.

Generally, for unimodal frequency distributions, the KS model Eq. (16) displays a simple transition between incoherence (asynchrony) and macroscopic synchronization at a critical finite value of $\varepsilon$ —but see Omel’chenko and Wolfrum (2012, 2016) for exceptions. The synchronous state is characterized by the appearance of a subset of oscillators that rotate with a common frequency and have their phases locked, thanks to the mutual coupling that is able to overcome the disparity of the natural frequencies. For the Lorentzian distribution of frequencies, Eq. (3), the critical coupling of the synchronization transition in the thermodynamic limit ($N\to\infty$) can be obtained analytically Sakaguchi and Kuramoto (1986):

[TABLE]

where the superscript “(av)” is used to emphasize that this is the critical coupling of the averaged model in Eq. (16). Curiously, within this approximation $\varepsilon_{c}$ does not depend on $q$. (This has to be attributed to the special properties of the Lorentzian distribution $g(\omega)$ which uses to yield particularly simple results in Kuramoto-like models.) The exact critical coupling of the Winfree model is computed numerically below and, as presumed, depends on $q$.

V.2 Synchronization threshold: Winfree vs. KS model

To test the goodness of the averaging approximation we compare the synchronization threshold of the Winfree model with the threshold of its averaged counterpart, given by Eq. (17). Our aim is to determine if certain pulses deviate more from the averaging approximation, and whether these results depend on the PRC offset $q$. To make the comparison significant we considered different pulse types with the same $\Pi$ values. In different panels of Fig. 1, pulses plotted with the same line style have identical shape factor $\Pi$, and therefore they yield the identical KS model upon averaging. In turn, the prediction of Eq. (17) is exactly the same for all pulse types (provided the same $\Pi$ value). In order to measure the deviation of the Winfree model from the the KS model we define the quantity

[TABLE]

which is proportional to the difference between the exact and the approximated critical couplings (normalized by the approximated critical coupling). For each pulse type and $q$ value, the locus of the Hopf bifurcation $\varepsilon_{H}(\Delta)$ is numerically available from the exact low-dimensional Eq. (11).

In Fig. 10 we graph $\rho$ for $\Pi=10/11$ and the three pulse types considered, adopting three values $q=-1$, [math], and $1$ in panels (a), (b) and (c), respectively. As expected $\rho(\Delta\to 0)=0$, indicating the validity of the averaging approximation in this limit. As $\Delta$ is increased from zero $\rho(\Delta)$ becomes positive for $q=-1$ and negative for $q=1$, which implies that synchronization is hindered (promoted) with respect to the averaging approximation for negative (positive) $q$. (This is also consistent with Fig. 9.) Numerical evidence shows that

[TABLE]

where $\xi(q)$ is a pulse-dependent odd function (and possibly monotonically decreasing). Note that this means that the best pulse type, in a certain sense, for $q>0$, becomes the worst for $q<0$. For instance, for $q=-1$ synchronizability is the best for the square pulse among the pulses considered, but this becomes just the opposite for $q=1$. The case $q=0$ (zero PRC offset) is the boundary between these scenarios, because $\xi(0)=0$. Accordingly, the inset of Fig. 10(b) confirms a nonlinear dependence, namely quadratic, of $\rho(\Delta\ll 1)$.

VI Conclusions

The Winfree model is broadly known, but scarcely studied in detail. The application of the OA ansatz to the Winfree model allows for the detailed investigation of its collective dynamics. Here we systematically investigated the dynamics of the Winfree model for three different pulse types and various widths, and for sinusoidal PRCs with positive, negative, and zero offsets $q$.

The case of negative offset was not considered in Pazó and Montbrió (2014), but has revealed to be interesting and nontrivial. The claim that narrow pulses are optimal for synchronizing large populations of oscillators Pazó and Montbrió (2014), does not hold for negative PRC offsets. In this case we observe that the optimal pulse, allowing to synchronize ensembles with a higher degree of heterogeneity, has an intermediate width, see Fig. 6(a). Moreover, for negative offsets (but not only) it is more likely to find a synchronization scenario with five codimension-two points (incl. BT’), in contrast to the scenario with three points reported in Pazó and Montbrió (2014). Under which conditions each scenario is found depends on the particular pulse type. From our results we inferred that pulses which are closer to zero at phases far from the peak phase are more likely to exhibit a BT’ point. In fact, the RP pulse does not exhibit a BT’ point for any $q$ value, while the square pulse already does for a balanced PRC ($q=0$). We also considered the limit of infinitely narrow pulses (Dirac delta pulses) and provided exact formulas for the synchronization boundary (a Hopf bifurcation). Additionally, we demonstrated that positive PRC offsets display larger synchronization regions and are capable of synchronizing more heterogeneous ensembles.

Finally, we have compared the synchronization threshold of the Winfree model with its averaging approximation (the Kuramoto-Sakaguchi model), and found what we may summarize as antisymmetry with the PRC offset parameter $q$.

In future studies, it would be interesting to find techniques to efficiently analyze the Winfree model with frequency distributions beyond the Lorentzian one. The study in Ariaratnam and Strogatz (2001) for a uniform distribution of natural frequencies is valuable, but it is difficult to extend it to nonvanishing $q$. Generalizing the model by considering other sources of heterogeneity is also an interesting venue for future research.

Acknowledgements.

We acknowledge support by MINECO (Spain) under Projects No. FIS2014-59462-P, No. FIS2016-74957-P, No. PSI2016-75688-P and No. PCIN-2015-127. E.M and D.P also acknowledge support by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 642563.