week ending
22 NOVEMBER 2013

PHYSICAL REVIEW LETTERS

PRL 111, 214101 (2013)

Phase Reduction Method for Strongly Perturbed Limit Cycle Oscillators
Wataru Kurebayashi,* Sho Shirasaka, and Hiroya Nakao
Graduate School of Information Science and Engineering, Tokyo Institute of Technology, O-okayama 2-12-1, Meguro,
Tokyo 152-8552, Japan(Received 30 April 2013; published 22 November 2013)
The phase reduction method for limit cycle oscillators subjected to weak perturbations has significantly
contributed to theoretical investigations of rhythmic phenomena. We here propose a generalized phasereduction method that is also applicable to strongly perturbed limit cycle oscillators. The fundamental
assumption of our method is that the perturbations can be decomposed into a slowly varying component ascompared to the amplitude relaxation time and remaining weak fluctuations. Under this assumption, we
introduce a generalized phase parameterized by the slowly varying component and derive a closed
equation for the generalized phase describing the oscillator dynamics. The proposed method enables us to
explore a broader class of rhythmic phenomena, in which the shape and frequency of the oscillation mayvary largely because of the perturbations. We illustrate our method by analyzing the synchronization
dynamics of limit cycle oscillators driven by strong periodic signals. It is shown that the proposed
method accurately predicts the synchronization properties of the oscillators, while the conventional
method does not.
DOI: 10.1103/PhysRevLett.111.214101

PACS numbers: 05.45.Xt, 05.45.a, 84.35.+i

Rhythmic phenomena are ubiquitous in nature and ofgreat interest in many fields of science and technology,
including chemical reactions, neural networks, genetic
circuits, lasers, and structural vibrations [1–8]. Theserhythmic phenomena often result from complex interactions among individual rhythmic elements, typically modeled as limit cycle oscillators. In analyzing such systems,
the phase reduction method [1–6] has been widely usedand considered an essential tool. It systematically approximates the high-dimensional dynamical equation of a perturbed limit cycle oscillator by a one-dimensional reduced
phase equation, with just a single phase variable  representing the oscillator state.
A fundamental assumption of the conventional phasereduction method is that the applied perturbation is sufficiently weak; hence, the shape and frequency of the limit
cycle orbit remain almost unchanged. However, this assumption hinders the applications of the method tostrongly perturbed oscillators, because the shapes and frequencies of their orbits can significantly differ from thosein the unperturbed cases. Indeed, strong coupling candestabilize synchronized states of oscillators that are stablein the weak coupling limit [9]. The effect of strong coupling can further lead to nontrivial collective dynamics
such as quorum-sensing transition [8], amplitude deathand bistability [9], and collective chaos [10]. The assumption of weak perturbations can also be an obstacle tomodeling real-world systems, which are often subjected
to strong perturbations.
Although the phase reduction method has recently been
extended to stochastic [11], delay-induced [12], and collective oscillations [13], these extensions are still limited to
the weakly perturbed regime. To analyze a broader class of
0031-9007=13=111(21)=214101(5)synchronization phenomena exhibited by strongly driven
or interacting oscillators, the conventional theory should be
extended. This Letter proposes an extension of the phasereduction method to strongly perturbed limit cycle oscillators, which enables us to derive a simple generalized
phase equation that quantitatively describes their dynamics. Although not all of the above collective phenomena[8–10] are the subject of discussion in this study, ourformulation will give an insight into a certain class ofthem, e.g., bistability between phase-locked and driftingstates [9]. We use our method to analyze the synchronization dynamics of limit cycle oscillators subjected to strong
periodic forcing, which cannot be treated appropriately by
the conventional method.
We consider a limit cycle oscillator whose dynamics
depends on a time-varying parameter IðtÞ ¼½I1 ðtÞ; . . . ; Im ðtÞ> 2 Rm representing general perturbations, described by_ ¼ FðXðtÞ; IðtÞÞ;
XðtÞ

(1)

where XðtÞ ¼ ½X1 ðtÞ; . . . ; Xn ðtÞ> 2 Rn is the oscillator
state and FðX; IÞ ¼ ½F1 ðX; IÞ; . . . ; Fn ðX; IÞ> 2 Rn is an
I-dependent vector field representing the oscillator dynamics. For example, X and I can represent the state of aperiodically firing neuron and the injected current, respectively [4,6]. In this Letter, we introduce a generalized phase, which depends on the parameter IðtÞ, of the oscillator. In
defining the phase , we require that the oscillator stateXðtÞ can be accurately approximated by using ðtÞ with asufficiently small error and that ðtÞ increases at a constantfrequency when the parameter IðtÞ remains constant.
The former requirement is a necessary condition for thephase reduction, i.e., for deriving a closed equation for the214101-1

Ó 2013 American Physical Society

PRL 111, 214101 (2013)

PHYSICAL REVIEW LETTERS

generalized phase, and the latter enables us to derive ananalytically tractable phase equation.
To define such , we suppose that I is constant until
further notice. We assume that Eq. (1) possesses a family of
stable limit cycle solutions with period TðIÞ and frequency!ðIÞ :¼ 2=TðIÞ for I 2 A, where A is an open subset ofRm (e.g., an interval between two bifurcation points). An
oscillator state on the limit cycle with parameter I can beparameterized by a phase  2 ½0; 2Þ as X0 ð; IÞ ¼
½X0;1 ð; IÞ; . . . ; X0;n ð; IÞ> . Generalizing the conventional
phase reduction method [1–5], we define the phase  suchthat, as the oscillator state XðtÞ ¼ X0 ððtÞ; IÞ evolves along
the limit cycle, the corresponding phase ðtÞ increases at a
_ ¼ !ðIÞ for each I 2 A.
constant frequency !ðIÞ as ðtÞWe assume that X0 ð; IÞ is continuously differentiable withrespect to  2 ½0; 2Þ and I 2 A.
We consider an extended phase space Rn  A, as depicted
schematically in Fig. 1(a). We define C  Rn  A as acylinder formed by the family of limit cycles [X0 ð; IÞ, I]for  2 ½0; 2Þ and I 2 A and define U  Rn  A as aneighborhood of C. For each I, we assume that any orbit
starting from an arbitrary point ðX; IÞ in U asymptoticallyconverges to the limit cycle X0 ð; IÞ on C. We can thenextend the definition of the phase into U, as in the conventional method [1–5], by introducing the asymptotic phaseand isochrons around the limit cycle for each I. Namely, we
can define a generalized phase function ðX; IÞ 2 ½0; 2Þof ðX; IÞ 2 U such that ðX; IÞ is continuously differentiable with respect to X and I, and ð@ðX; IÞ=@XÞ FðX; IÞ ¼ !ðIÞ holds everywhere in U, where ð@=@XÞ ¼½ð@=@X1 Þ; . . . ; ð@=@Xn Þ> 2 Rn is the gradient ofweek ending
22 NOVEMBER 2013

ðX; IÞ with respect to X and the dot ðÞ denotes an innerproduct. This ðX; IÞ is a straightforward generalizationof the conventional asymptotic phase [1–5] and guarantees
that the phase of any orbit XðtÞ in U always increases_
with a constant frequency as ðXðtÞ;IÞ ¼ !ðIÞ at each I.
For any oscillator state on C, ðX0 ð; IÞ; IÞ ¼  holds. In
general, the origin of the phase can be arbitrarily defined foreach I as long as it is continuously differentiable withrespect to I. The assumptions that X0 ð; IÞ and ðX; IÞare continuously differentiable can be further relaxed for acertain class of oscillators, such as those considered inRef. [14].
Now suppose that the parameter IðtÞ varies with time. To
define  that approximates the oscillator state with a sufficiently small error, we assume that IðtÞ can be decomposed
into a slowly varying component qðtÞ 2 A and remaining
weak fluctuations pðtÞ 2 Rm as IðtÞ ¼ qðtÞ þ pðtÞ.
Here, the parameters  and  are assumed to be sufficiently
small so that qðtÞ varies slowly as compared to the relaxation time of a perturbed orbit to the cylinder C of thelimit cycles, which we assume to be Oð1Þ without lossof generality, and the oscillator state XðtÞ always remains
in a close neighborhood of X0 ð; qðtÞÞ on C, i.e.,XðtÞ ¼ X0 ððtÞ; qðtÞÞ þ Oð; Þ holds (see Supplemental
Material [15]). We also assume that qðtÞ is continuously
differentiable with respect to t 2 R. Note that the slow
component qðtÞ itself does not need to be small.
Using the phase function ðX; IÞ, we introduce a generalized phase ðtÞ of the limit cycle oscillator (1) asðtÞ ¼ ðXðtÞ; qðtÞÞ. This definition guarantees that ðtÞ
increases at a constant frequency when IðtÞ remains

FIG. 1 (color online). Phase dynamics of a modified Stuart-Landau oscillator. (a) A schematic diagram of the extended phase space
Rn  A with n ¼ 2 and m ¼ 1. (b) Frequency !ðIÞ. (c) I-dependent stable limit cycle solutions X0 ð; IÞ. (d),(e) Sensitivity functionsð; IÞ and ð; IÞ. (f),(g) Time series of the phase ðtÞ of the oscillator driven by (f) a periodically varying parameter Ið1Þ ðtÞ or (g) achaotically varying parameter I ð2Þ ðtÞ. For each of these cases, results of the conventional (top panel) and proposed (middle panel)~ ¼ ðXðtÞ; qc Þ and the generalized phase ðtÞ ¼ ðXðtÞ; qðtÞÞ measured
methods are shown. Evolution of the conventional phase ðtÞfrom the original system (lines) is compared with that of the conventional and generalized phase equations (circles). Time series of thestate variable xðtÞ (red solid line) and time-varying parameter IðtÞ (blue dashed line) are also depicted (bottom panel). The periodicallyvarying parameter is given by I ð1Þ ðtÞ ¼ qð1Þ ðtÞ þ pð1Þ ðtÞ with qð1Þ ðtÞ ¼ 0:05 sinð0:5tÞ þ 0:02 sinðtÞ and pð1Þ ðtÞ ¼ 0:02 sinð3tÞ, andthe chaotically varying parameter is given by I ð2Þ ðtÞ ¼ qð2Þ ðtÞ þ pð2Þ ðtÞ with qð2Þ ðtÞ ¼ 0:007L1 ð0:3tÞ and pð2Þ ðtÞ ¼ 0:001L2 ðtÞ,where L1 ðtÞ and L2 ðtÞ are independently generated time series of the variable x of the chaotic Lorenz equation [3], x_ ¼ 10ðy  xÞ,
y_ ¼ xð28  zÞ  y, and z_ ¼ xy  8z=3.

214101-2

PRL 111, 214101 (2013)

week ending
22 NOVEMBER 2013

PHYSICAL REVIEW LETTERS

constant and leads to a closed equation for ðtÞ. Expanding_ ¼ FðX; qðtÞÞ þ GðX; qðtÞÞpðtÞ þ
Eq. (1) in  as XðtÞ_ ¼
Oð2 Þ and using the chain rule, we can derive ðtÞ!ðqðtÞÞ þ ð@ðX;IÞ=@XÞjðXðtÞ;qðtÞÞ  GðX;qðtÞÞpðtÞ þ_þ Oð2 Þ, where GðX; IÞ 2
ð@ðX; IÞ=@IÞjðXðtÞ;qðtÞÞ  qðtÞnmR
is a matrix whose ði; jÞth element is given
by @Fi ðX; IÞ=@Ij , @=@I ¼ ½ð@=@I1 Þ; . . . ; ð@=@Im Þ>
2 Rm is the gradient of ðX; IÞ with respect to I, and_qðtÞ
denotes dqðtÞ=dðtÞ.
To obtain a closed equation for , we use the lowest-orderapproximation in  and , i.e., XðtÞ ¼ X0 ððtÞ; qðtÞÞ þOð; Þ. Then, by defining a phase sensitivity functionZð; IÞ ¼ ð@ðX; IÞ=@XÞjðX0 ð;IÞ;IÞ 2 Rn and two othersensitivity functions  ð; IÞ ¼ G> ðX0 ð; IÞ; IÞZð; IÞ 2
Rm and ð; IÞ ¼ ð@ðX; IÞ=@IÞjðX0 ð;IÞ;IÞ 2 Rm , we canobtain a closed equation for the oscillator phase ðtÞ as_ ¼ !ðqðtÞÞ þ  ð; qðtÞÞ  pðtÞðtÞ_þ ð; qðtÞÞ  qðtÞþ Oð2 ; 2 ; Þ;(2)

which is a generalized phase equation that we propose inthis study. The first three terms in the right-hand side ofEq. (2) represent the instantaneous frequency of the oscillator, the phase response to the weak fluctuations pðtÞ, andthe phase response to deformation of the limit cycle orbit
caused by the slow variation in qðtÞ, respectively, all ofwhich depend on the slowly varying component qðtÞ.
To address the validity of Eq. (2) more precisely, let ðIÞ( > 0) denote the absolute value of the second largest
Floquet exponent of the oscillator for a fixed I, whichcharacterizes the amplitude relaxation time scale ofthe oscillator (  1=ðIÞ). As argued in Supplemental
Material [15], we can show that the error terms in Eq. (2)remain sufficiently small when =ðqðtÞÞ  1 and=ðqðtÞÞ2  1, namely, when the orbit of the oscillator
relaxes to the cylinder C sufficiently faster than the variations in qðtÞ.
~ ¼
Note that if we define the phase variable as ðtÞðXðtÞ; qc Þ with some constant qc instead of ðtÞ ¼~ gives the conventional phase. Then,ðXðtÞ; qðtÞÞ, ðtÞ~_ ¼ !c þ
we obtain the conventional phase equation ðtÞ~  pðtÞ þ Oð2 Þ with qðtÞ ¼ qc and pðtÞ ¼
 c ðÞIðtÞ  qc . Here, !c :¼ !ðqc Þ is a natural frequency,~ ¼  ð;~ qc Þ ¼ GðX0 ð;~ qc Þ; qc Þ> Zð;~ qc Þ, and Zð;~ qc Þ
 c ðÞis the conventional phase sensitivity function at I ¼ qc [2].
This equation is valid only when =ðqc Þ  1 [i.e.,kIðtÞ  qc k=ðqc Þ  1]. By using the near-identity transformation [16], we can show that the conventional equationis actually a low-order approximation of the generalized
equation (2) (see Sec. III of Supplemental Material [15]).
In practice, we need to calculate  ð; IÞ and ð; IÞnumerically from mathematical models or estimate themthrough experiments. We can show that the following
relations hold (see Supplemental Material [15] for thederivation): ð; IÞ ¼  ð; IÞ ¼ ð0 ; IÞ 1 ðIÞ :¼
2

@X0 ð; IÞ>Zð; IÞ;
@I

(3)

1 Z
½ ð0 ; IÞ   ðIÞd0 ; (4)!ðIÞ 0Z 2
0 ð; IÞd ¼

d!ðIÞ
;
dI

(5)

where ð@X0 ð; IÞ=@IÞ 2 Rnm is a matrix whose ði; jÞthelement is given by ð@X0;i ð; IÞ=@Ij Þ, 0 2 ½0; 2Þ is aconstant, and  ðIÞ is the average of ð; IÞ with respectto  over one period of oscillation. From mathematicalmodels of limit cycle oscillators, Zð; IÞ can be obtained
numerically by the adjoint method for each I [5,6],
and then  ð; IÞ and ð; IÞ can be computed from ð; IÞ ¼ G> ðX0 ð; IÞ; IÞZð; IÞ and Eqs. (3) and (4).
Experimentally, Zð; IÞ and  ð; IÞ can be measured byapplying small impulsive perturbations to I, while ð; IÞ
can be obtained by applying small stepwise perturbations
to I.
To test the validity of the generalized phase equation (2),
we introduce an analytically tractable model, a modified
Stuart-Landau (MSL) oscillator (see Ref. [17] and Fig. 1
for the definition and details). We numerically predict thephase ðtÞ of a strongly perturbed MSL oscillator by bothconventional and generalized phase equations and comparethem with direct numerical simulations of the original
system. In applying the conventional phase reduction,
we set qc ¼ hIðtÞit , where hit denotes the time average.
In Fig. 1, we can confirm that the generalized phaseequation (2) accurately predicts the generalized phaseðXðtÞ; qðtÞÞ of the original system, while the conventional phase equation does not well predict the conventional phase ðXðtÞ; qc Þ because of large variations in IðtÞ.
As an application of the generalized phase equation (2),
we analyze k:l phase locking [18] of the system (1) to aperiodically varying parameter IðtÞ with period TI and_ t¼
frequency !I , in which the frequency tuning (lhik!I ) occurs. Although the averaging approximation [19]for the phase difference c~ ðtÞ ¼ lðtÞ  k!I is generally
used to analyze the phase locking [2,18], we cannot
directly apply it in the present case because the frequency!ðqðtÞÞ can vary largely with time. Thus, generalizing theconventional definition, we introduce the phase difference
as c ðtÞ ¼ lðtÞ  k!I t  lhðtÞ with an additional termlhðtÞ to remove the large periodic variations in c ðtÞdue to !ðqðtÞÞ, wherehðtÞ is a TI R
-periodic function
Rdefined as hðtÞ ¼ t0 ½!ðqðt0 ÞÞ  TI1 T0 I !ðqðtÞÞdtdt0 .
By virtue of this term, temporal variations in c_ remainof the order Oð; Þ, i.e., j c_ j  1, which enables us toapply the averaging approximation to c .
Introducing a small parameter  representing the magnitude
of variations in c , one can derive a dynamical equation for c

214101-3

PRL 111, 214101 (2013)

PHYSICAL REVIEW LETTERS

as c_ ðtÞ ¼ fð c ; tÞ, where fð c ; tÞ ¼ lgð c =l þ k!I t=l þ_ and gð; tÞ denotes the right-hand sidehðtÞ; tÞ  k!I  lhðtÞ
of Eq. (2). Using first- and second-order averaging [19], we
can introduce slightly deformed phase differences c 1;2 satisfying c 1;2 ðtÞ ¼ c ðtÞ þ OðÞ and obtain the first- andsecond-order averaged equations c_ 1 ðtÞ ¼ f1 ð c 1 Þ þOð2 Þ and c_ 2 ðtÞ ¼ f1 ð c 2 Þ þ 2 f2 ð c 2 Þ þ Oð3 Þ, wheref1 ð c Þ and f2 ð c Þ are given by f1 ð c Þ ¼ ðlTI Þ1 
RlTI
R I
and
f2 ð c Þ ¼ ðlTI Þ1 lT
0 fð c ; tÞdt0 ½uð c ; tÞ 
1 ð c Þð@uð c ; tÞ=@ c Þdt,
ð@fð c ; tÞ=@ c Þ frespectively,Rt00and uð c ; tÞ ¼ 0 ½fð c ; t Þ  f1 ð c Þdt . These averaged
equations can be considered autonomous by neglecting theOð2 Þ and Oð3 Þ terms, respectively. Averaged equations forthe conventional phase equation can be derived similarly.
Thus, if the averaged equation has a stable fixed point, k:lphase locking is expected to occur. As demonstrated below, thefirst-order averaging of the generalized phase equation alreadypredicts qualitative features of the phase-locking dynamics,
while the second-order averaging gives more precise results
when the parameter IðtÞ varies significantly.
As an example, we use the MSL oscillator and investigate their phase locking to periodic forcing (SeeSupplemental Material [15] for other examples). Figure 2
shows the results of the numerical simulations. We apply

week ending
22 NOVEMBER 2013

four types of periodically varying parameters and predict if
the oscillator exhibits either 1:1 or 1:2 phase locking to theperiodically varying parameter qðtÞ [a small fluctuationpðtÞ is also added for completeness]. We derive averaged
equations for the phase differences c 1;2 by using the proposed and conventional methods and compare the resultswith direct numerical simulations of the MSL oscillator. We
find that our new method correctly predicts the stablephase-locking point already at first-order averaging, while
the conventional method does not. In particular, the conventional method can fail to predict whether phase locking
takes place or not, as shown in Figs. 2(g) and 2(h), evenafter the second-order averaging. In this case, the exponential dependence of the frequency !ðIÞ on the parameter I is
the main cause of the breakdown of the conventional
method (see Sec. III of Supplemental Material [15] for a
discussion). Typical trajectories of ½xðtÞ; yðtÞ; qðtÞ> areplotted on the cylinder C of limit cycles in the extended
phase space ½x; y; I> , which shows that the oscillator statemigrates over C synchronously with the periodic forcing.
The trajectories are closed when phase locking occurs.
In summary, we proposed a generalized phase reduction
method that enables us to theoretically explore a broader
class of strongly perturbed limit cycle oscillators. Althoughstill limited to slowly varying perturbations with weak

FIG. 2 (color online). Phase locking of the modified Stuart-Landau oscillator. Four types of periodically varying parameters IðjÞ(j ¼ 3, 4, 5, 6) are applied, which lead to 1:1 phase locking to Ið3Þ ðtÞ [(a), (e), and (i)], 1:1 phase locking to Ið4Þ ðtÞ [(b), (f), and (j)], 1:2
phase locking to Ið5Þ ðtÞ [(c), (g), and (k)], and failure of phase locking to Ið6Þ ðtÞ [(d), (h), and (l)]. (a)–(d) Time series of the statevariable xðtÞ of a periodically driven oscillator (red solid line) and periodic external forcing (blue dashed line). (e)–(h) Dynamics of thephase difference c 1;2 with an arrow representing a stable fixed point (top panel) and time series of c 1;2 with 20 different initial states(bottom panel). (i)–(l) Orbits of a periodically driven oscillator (blue line) on the cylinder of the limit cycles (light blue line) plotted inðjÞthe extended phase space. The parameter IðjÞ ðtÞ is given by IðjÞ ðtÞ ¼ qðjÞ ðtÞ þ pðjÞ ðtÞ, qðjÞ ðtÞ ¼ ðjÞ sinð!ðjÞ
I tÞ, and p ðtÞ ¼ðjÞ
ð3;4;5;6Þ
ð3;4;5;6Þ
¼ 0:1, 0.3, 0.4, 0.4 and !I¼ 1:05, 1.10, 0.57, 0.51.
0:02 sinð5!I tÞ with

214101-4

PRL 111, 214101 (2013)

PHYSICAL REVIEW LETTERS

fluctuations, our method avoids the assumption of weakperturbations, which has been a major obstacle in applying
the conventional phase reduction method to real-world
phenomena. It will therefore facilitate further theoreticalinvestigations of nontrivial synchronization phenomena ofstrongly perturbed limit cycle oscillators [9,10]. As a final
remark, we point out that a phase equation similar toEq. (2) has been postulated in a completely different context, to analyze the geometric phase in dissipative dynamical systems [20]. This formal similarity may provide aninteresting possibility of understanding synchronization
dynamics of strongly perturbed oscillators from a geometrical viewpoint.
Financial support by JSPS KAKENHI (25540108 and
22684020), CREST Kokubu project of JST, and FIRST
Aihara project of JSPS are gratefully acknowledged.

*kurebayashi.w.aa@m.titech.ac.jp[1] A. T. Winfree, The Geometry of Biological Time (Springer,
New York, 2001).
[2] Y. Kuramoto, Chemical Oscillations, Waves andTurbulence (Dover, New York, 2003).
[3] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: A Universal Concept in Nonlinear
Sciences (Cambridge University Press, Cambridge,
England, 2001).
[4] F. C. Hoppensteadt and E. M. Izhikevich, Weakly
Connected Neural Networks (Springer, New York, 1997).
[5] G. B. Ermentrout and D. H. Terman, Mathematical
Foundations of Neuroscience (Springer, New York,
2010).
[6] E. Brown, J. Moehlis, and P. Holmes, Neural Comput. 16,
673 (2004).
[7] K. Wiesenfeld, C. Bracikowski, G. James, and R. Roy,
Phys. Rev. Lett. 65, 1749 (1990); S. H. Strogatz, D. M.
Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature(London) 438, 43 (2005); I. Z. Kiss, C. G. Rusin, H. Kori,
and J. L. Hudson, Science 316, 1886 (2007).
[8] J. Garcia-Ojalvo, M. B. Elowitz, and S. H. Strogatz, Proc.
Natl. Acad. Sci. U.S.A. 101, 10955 (2004); S. De Monte,
F. d’Ovidio, S. Danø, and P. G. Sørensen, Proc. Natl.
Acad. Sci. U.S.A. 104, 18377 (2007); A. F. Taylor, M. R.
Tinsley, F. Wang, Z. Huang, and K. Showalter, Scienceweek ending
22 NOVEMBER 2013

323, 614 (2009); T. Danino, O. Mondragón-Palomino, L.
Tsimring, and J. Hasty, Nature (London) 463, 326 (2010).
[9] D. G. Aronson, G. B. Ermentrout, and N. Kopell, Physica(Amsterdam) 41D, 403 (1990); R. E. Mirollo and S. H.
Strogatz, J. Stat. Phys. 60, 245 (1990); D. Hansel, G.
Mato, and C. Meunier, Neural Comput. 7, 307 (1995);
I. Z. Kiss, W. Wang, and J. L. Hudson, J. Phys. Chem. B
103, 11433 (1999); P. C. Bressloff and S. Coombes, NeuralComput. 12, 91 (2000); Y. Zhai, I. Z. Kiss, and J. L.
Hudson, Phys. Rev. E 69, 026208 (2004).
[10] V. Hakim and W. J. Rappel, Phys. Rev. A 46, R7347(1992); N. Nakagawa and Y. Kuramoto, Prog. Theor.
Phys. 89, 313 (1993); H. Nakao and A. S. Mikhailov,
Phys. Rev. E 79, 036214 (2009).
[11] K. Yoshimura and K. Arai, Phys. Rev. Lett. 101, 154101(2008); J.-N. Teramae, H. Nakao, and G. B. Ermentrout,Phys. Rev. Lett. 102, 194102 (2009); D. S. Goldobin, J.-N.
Teramae, H. Nakao, and G. B. Ermentrout, Phys. Rev.
Lett. 105, 154101 (2010).
[12] V. Novičenko and K. Pyragas, Physica (Amsterdam)
241D, 1090 (2012); K. Kotani, I. Yamaguchi, Y. Ogawa,Y. Jimbo, H. Nakao, and G. B. Ermentrout, Phys. Rev.
Lett. 109, 044101 (2012).
[13] Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y.
Kuramoto, Phys. Rev. Lett. 101, 024101 (2008); Y.
Kawamura, H. Nakao, and Y. Kuramoto, Phys. Rev. E
84, 046211 (2011).
[14] E. M. Izhikevich, SIAM J. Appl. Math. 60, 1789 (2000).
[15] See Supplemental Material at http://link.aps.org/
supplemental/10.1103/PhysRevLett.111.214101 for a
detailed derivation and discussion.
[16] J. P. Keener, Principles of Applied Mathematics:
Transformation and Approximation (Addison-Wesley,
Boston, 1988).
[17] The modified Stuart-Landau oscillator has a twodimensional state variable XðtÞ ¼ ½xðtÞ; yðtÞ> and a vector field FðX; IÞ ¼ ½e2I ðx  y  IÞ  ððx  IÞ2 þ y2 Þðx IÞ; e2I ðx þ y  IÞ  ððx  IÞ2 þ y2 Þy> with ðX; IÞ ¼
!ðIÞ ¼ e2I ,X0 ð; IÞ ¼ ½I þtan1 ½y=ðx  IÞ,I
I>
e cos; e sin , ð; IÞ ¼ eI sin, and ð; IÞ ¼
2e2I  eI cos.
[18] G. B. Ermentrout, J. Math. Biol. 12, 327 (1981).
[19] J. A. Sanders and F. Verhulst, Averaging Methods in
Nonlinear Dynamical Systems (Springer-Verlag, New
York, 1985).
[20] T. B. Kepler and M. L. Kagan, Phys. Rev. Lett. 66, 847(1991).

214101-5