PRL 111, 024102 (2013)

week ending
12 JULY 2013

PHYSICAL REVIEW LETTERS

Optimal Waveform for Fast Entrainment of Weakly Forced Nonlinear Oscillators
Anatoly Zlotnik,1,* Yifei Chen,2,† István Z. Kiss,2,‡ Hisa-Aki Tanaka,3,§ and Jr-Shin Li1,∥
1

Electrical and Systems Engineering, Washington University in St. Louis, St. Louis, Missouri 63130, USA
2
Department of Chemistry, Saint Louis University, St. Louis, Missouri 63103, USA
3
Department of Electronic Engineering, The University of Electro-Communications, Tokyo 182-8585, Japan(Received 24 January 2013; revised manuscript received 2 April 2013; published 9 July 2013)
For many biological and engineered systems, a central function or design goal is to abbreviate the timerequired to synchronize a rhythmic process to an external forcing signal. We present a theory for deriving
the input that effectively minimizes the average transient time required to entrain a phase model, whichenables a practical technique for constructing fast entrainment waveforms for general nonlinear oscillators. This result is verified in numerical simulations using the Hodgkin-Huxley neuron model, and inexperiments on an oscillatory electrochemical system.
DOI: 10.1103/PhysRevLett.111.024102

PACS numbers: 87.19.lr, 05.45.Xt, 82.45.Qr, 87.19.ln

The entrainment process is fundamental to numerous
scientific and engineering applications in which oscillating
systems are asymptotically synchronized to an external
periodic signal [1,2]. The ability to optimize entrainment
has important implications for achieving rapid cardiac
resynchronization [3] and quick adjustment from jet lag[4], maximizing the growth rate of plants [5], and implementing phase-locked loop circuits and injection-locked
microintegrated oscillators [6]. When the weak perturbation approximation is made, a rescaling of the phase
response curve (PRC) was shown to be the minimum
energy signal for spiking or entraining oscillators at a given
period [7–9], and a weighted sum of appropriately shifted
PRCs maximizes the range of frequency detunings forwhich entrainment occurs [10,11]. An alternative essential
objective is to minimize the time to entrainment at a given
forcing signal energy, in order to establish a fixed phaserelationship between the system and forcing signal assoon as possible after the forcing is applied [12]. This
notion of fast entrainment can also be used to minimize
the time required to reestablish entrainment after interruptions caused by disturbances [13].
In this Letter, we use phase model reduction to derive
an asymptotically optimal waveform that maximizes theaverage rate of entrainment for general weakly forced
nonlinear oscillators. The rate of entrainment is characterized by the coefficient of exponential decay in thephase difference between the system and forcing signal.
We present a theory by which the entrainment time scaleis minimized for a specified forcing energy, where theoptimal waveform is a sum of the PRC and its derivativewith weights that depend on the difference between
the natural and forcing frequencies. These findings canbe applied to weakly nonlinear oscillators just past the
Hopf bifurcation, as well as strongly nonlinear relaxation
oscillators. We confirm our results with numerical simulations using the Hodgkin-Huxley (HH) neuron model, as
well as in experiments on an oscillatory chemical system0031-9007=13=111(2)=024102(5)arising through the electrodissolution of nickel in sulfuric
acid.
Phase coordinate transformation is a model reduction
technique that is useful for examining nonlinear oscillating
systems [14,15], and can also be used for system identification when the dynamics are complex or unknown [2].
Such models have been studied extensively, with a particular focus on neural [14,16] and electrochemical [17–19]
systems. Consider a full state-space model described by asmooth ordinary differential equation system x_ ¼ fðx; uÞ,xð0Þ ¼ x0 , where xðtÞ 2 Rn is the state and uðtÞ 2 R isa control, with an attractive, nonconstant limit cycleðtÞ ¼ ðt þ TÞ that satisfies _ ¼ fð; 0Þ on the periodic
orbit  ¼ fy 2 Rn : y ¼ ðtÞ for 0  t < Tg  Rn . Thissystem is reduced to a scalar phase equationc_ ¼ ! þ Zð c Þu;(1)where ! is the natural frequency of oscillation, Z isa smooth PRC, and c ðtÞ is the asymptotic phase [20].
The model is valid for inputs u such that the state-spacesystem remains within a neighborhood U of  [21], and the
PRC can be computed numerically [14,22,23].
The primary objective in entrainment design is to lockthe system to an input with the desired frequency  using
a control uðtÞ ¼ vðtÞ where v is 2 periodic. We makethe weak perturbation approximation, i.e., v ¼ "v1 where
v1 has unit energy and "  1, so that given this control
the actual state of the system is guaranteed to remain in U,
and the phase model (1) remains valid. We define a slow
phase variable ðtÞ ¼ c ðtÞ  t that satisfies the dynamic
equation _ ¼ c_   ¼ ! þ Zðt þ ÞvðtÞ, where! ¼ !   denotes the detuning between the naturaland forcing frequencies. To study the asymptotic behaviorof the slow phase, we eliminate the explicit dependence ontime on the right-hand side by using formal averaging [17].
Given a periodic forcing with frequency , we denotethe forcing phase  ¼ t. We also define an averaging
operator hi: P ! R on the set of 2-periodic functions by

024102-1

Ó 2013 American Physical Society

week ending

PHYSICAL REVIEW LETTERS
12 JULY 2013
PRL 111, 024102 (2013)
R2
The constraints can be adjoined to the objective usinghxi ¼ 1=ð2Þ 0 xðÞd. The weak ergodic theorem for
multipliers  and  to yield the formulation
measure-preserving dynamical systems on the torus [24]implies that for any periodic function v, the interaction
J ½v ¼ 0v ð’ Þ þ ðhv2 i  PÞ þ ðv ð’ Þ þ !Þ
function¼ hZ0 ð þ ’ ÞvðÞi þ ðhv2 i  PÞv ðÞ ¼ hZð þ ÞvðÞiþ ðhZð þ ’ ÞvðÞi þ !Þ
1 Z 2
¼
Zð þ ÞvðÞd
1 Z 2
2 0ðvðÞ½Zð þ ’ Þ  Z0 ð þ ’ Þ þ vðÞ
¼
2 0
1 ZT
Zðt þ ÞvðtÞdt
(2)
¼ lim
 P þ !Þd:(7)
T!1 T 0is a smooth function in P . Using the weak perturbation
approximation and the formal averaging theorem [25],the time-averaged slow phase dynamics are, up to Oð"2 Þ,’_ ¼ ! þ v ð’Þ:(3)

Equation (3) is used to study the asymptotic behavior of (1)under periodic forcing. We say that the system is entrainedby a control u ¼ vðtÞ when (3) satisfies ’_ ¼ 0, which
eventually occurs if there exists a phase ’ satisfying
! þ v ð’ Þ ¼ 0. For nonzero control waveform v andnonzero PRC Z, the function v ð’Þ is not identically zero,
so when the system is entrained there exists at least
one ’ 2 ½0; 2Þ that is an attractive fixed point of (3).
We have shown that the minimum energy periodic waveform that entrains a single oscillator with natural frequency! to a target frequency  is a rescaling of the PRC, givenby vðÞ ¼ !ZðÞ=hZ2 i [9,11].
Our goal here is to entrain the system (1) to a target
frequency  as quickly as possible by using a periodic
control v of fixed power P ¼ hv2 i. Ideally, the interaction
function would be of a piecewise-constant form, so that the
averaged slow phase ’ converges to a fixed point ’ at a
uniform rate from any initial value. However, the discontinuity as ’ ! ’ would result in a singularity in the controlv, making it infeasible in practice. An alternative is tomaximize j’_  j, the rate of convergence of the averaged
slow phase in the neighborhood of its attractive fixed point’ . The calculus of variations can then be used to obtain asmooth optimal candidate solution that also performs wellin practice. When the system (3) is entrained by a controlv, there exists an attractive fixed point ’ satisfying
v ð’ Þ þ ! ¼ 0 and 0v ð’ Þ < 0, where 0 is the differentiation operator. In order to maximize the rate of entrainment in a neighborhood of ’ using a control of power P,
the value of j’j_ should be maximized for values of ’ near’ , which occurs when 0v ð’ Þ is large. This results inthe following problem formulation for fast entrainment:
maxJ ½v ¼ 0v ð’ Þ(4)s:t: hv2 i ¼ P(5)

v ð’ Þ þ ! ¼ 0:

(6)

v2P

The associated Euler-Lagrange equation isZð þ ’ Þ  Z0 ð þ ’ Þ þ 2vðÞ ¼ 0;

(8)and solving for v yields the candidate solutionvðÞ ¼

1 0½Z ð þ ’ Þ  Zð þ ’ Þ:2

(9)The multipliers  and  can be found by substituting (9)into the constraints (5) and (6). This yields the equations1½hðZ0 Þ2 i  2hZ0 Zi þ 2 hZ2 i ¼ P;
42

(10)

1
½hZ0 Zi  hZ2 i ¼ !:(11)2Because Z is 2 periodic, one can show, e.g., using
Fourier series, that hZ0 Zi ¼ 0, so that (11) easily yields
 ¼ 2!=hZ2 i. Substituting this result into (10) leads toa quadratic equation (10) for . Now, by substituting (9)into 0 ð’ Þ ¼ hZð þ ’ ÞvðÞi we obtain 0 ð’ Þ ¼
hðZ0 Þ2 i=ð2Þ, so we choose  < 0 when solving (10) for in order to maximize the objective in (4). Thus the
optimal waveform and multiplier simplify tov
uﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Z0 ðÞ !ZðÞ
1u
hðZ0 Þ2 iuvðÞ ¼
(12)

;

¼

t
2;
2
2 P  ð!Þ
hZ2 i
2
hZ i
where we disregard the phase shift ’ , because the entrainment process is asymptotic. For zero frequency detuning,
the optimal waveform is a rescaling of the derivative Z0of the PRC. As j!j increases, v continuously transforms
towards Z, which is the minimum energy waveform forfrequency control [9]. This transition reflects the conceptual tradeoff between the fast entrainment objective (4) andfrequency control constraint (6), which can be satisfied
only when P > ð!Þ2 =hZ2 i.
Consider a system with a sinusoidal PRC, given byZðÞ ¼ a sinðÞ. Using angle sum identities and the factthat in this case hðZ0 Þ2 i ¼ hZ2 i, one can show that v isof form vðÞ ¼ P sinðÞ. Indeed, for the case of a sinusoidal PRC, a sinusoidal input optimizes the minimum energy[9] and rapid phase-locking objectives simultaneously.
However, the utility of our approach is most evident for
oscillating systems with complex dynamics, in particular
those that exhibit relaxation, and hence higher harmonics

024102-2

PHYSICAL REVIEW LETTERS

0.25
0.2
0.15
0.1
0.05
0
−0.05
−0.1

0.5

v (ψ)

0

*

Z(ψ)

PRL 111, 024102 (2013)

−0.5

0

2

4

ψ

−1

6

0

2

ψ

4

6

FIG. 1. (a) PRC of the HH model; (b) Optimal fast entrainment
waveform v for HH system at  ¼ ! and P ¼ 0:1 mW.

in the PRC. As an example, consider the HH system [26],which is a fundamental model used in the study of neuraldynamics [11]. When the baseline current Ib injected into
the axon is sufficiently high, the voltage V spikes repeatedly. Our goal is to modulate the additional injected currentIðtÞ to entrain the spiking frequency to a desired target  inas short a time as possible. We first reduce the HH system
to a phase model as in (1) where u ¼ IðtÞ, where!  0:429 rad=sec, and the PRC Z is given in Fig. 1(a).
After selecting the control power P and the target frequency , we use (12) to compute the optimal waveformv, which is shown in Fig. 1(b). Numerically, we use theFourier series coefficients of Z to evaluate expressions
derived from the PRC, such as Z0 , hZ2 i, and so on. In thiscomputational example, we focus on initial convergence
rates for fast entrainment, which can be quantified by therate k at which the phase difference between successiveinterspike intervals converges exponentially to zero,according to ’n ¼ ekn , as shown in Fig. 2. The optimal
0.3

10
sin
prc
opt

∆φ

n

0.1
0

−0.1
−0.2

sin
prc
opt

8

log|∆φn|

0.2

6
4
2

1

2

3

4

0

5

1

2

n

3

4

5

n

1.02

week ending
12 JULY 2013

waveform (12) achieves a significantly greater average rate
k for all values of  and initial states on .
The experimental utility of the phase model techniquefor fast entrainment is demonstrated by manipulating
an oscillatory chemical process [27]. A standard threeelectrode setup was used with a 1 mm diameter nickel
working, a Hg=Hg2 SO4 =ðsatÞK2 SO4 reference, and a Pt
coated Ti rod counterelectrode immersed in 3 mol/L sulfuric acid solution at 10 C. The nickel working electrodewas polarized with a potentiostat (Gamry Instruments,
Reference 600) at a circuit potential V ¼ V0 þ AFðÞ,where A and F are the forcing amplitude and waveform,
respectively, and V0 is the base potential. Each forcing
waveform F has power P ¼ 0:5. The current, proportional
to the dissolution rate, was measured by the potentiostat
at a rate of 200 Hz. When 1 k resistance was attached to
the nickel wire, nonlinear current oscillations with a periodof 2.11 s were obtained at V0 ¼ 1:15 V as shown in theinset of Fig. 3. In each instance of the experiment, the PRC,
such as the example in Fig. 3(a), was obtained using thepulse perturbation method [10,28]. The phase of the oscillation was obtained using the linear interpolation technique[2] by setting the phase of the nth current peak to 2n. The
transformation of the PRC as the circuit potential increases
has been previously analyzed in detail [28].
Using (12), an optimal fast entrainment waveform
was constructed for equal forcing and natural frequencies,
 ¼ !, in order to remove the effect of the frequency
control constraint (6). When the optimal waveform withamplitude of A ¼ 12:5 mV was applied to entrain the freerunning chemical oscillator, the phase difference between
the current oscillations and the forcing signal, shown in
Fig. 4(a), monotonically decreased until a final phasedifference of f ¼ 5:51 rad was attained after 30 sec.
The behavior of the phase difference ðtÞ¼ðtÞfafter 20 sec can be closely described by an exponential
decay ln½ðtÞ=ð0Þ ¼ kt, which is shown in
Fig. 4(b), with a rate of entrainment k ¼ 0:243 s1 for

2.5
2

0.6

1.5
1

1

k

Ω/ω

1.01

0.4

0.5
0.99

0

0.2

−0.5
0.98
1

2

3

θ

4

5

6

−1

0

sin

prc

optimalFIG. 2 (color online). Simulations with the HH model:(a) Convergence of phase difference ’n between interspikeintervals; (b) Exponential fit for k when  ¼ 1:01! and P ¼
0:5. The k values are 0.5415, 0.9510, and 1.4139 for sine, PRC,
and optimal waveforms, respectively. (c) Initial convergence rates
k (in color) for 5 cycles with  2 ½0:98!; 1:02! and ð0Þ 2½0; 2 when using the optimal waveform. (d) Average initial k onð; ð0ÞÞ 2 ½0:98!; 1:02! ½0; 2 for sine, PRC, and optimal
waveforms is 0.3933, 0.5365, and 0.7691, respectively. Initial
divergence takes place in 7.69%, 12.73%, and 6.04% of initialconditions for sine, PRC, and optimal waveforms, respectively.

FIG. 3 (color online). Electrodissolution experiments: (a) PRC
and current waveform (inset) of the electrochemical oscillations.
The PRC is measured by stimulating the system using a sequence of pulses (A ¼ 200 mV magnitude and  ¼ 0:05 s pulse
width) and measuring the corresponding phase shift () as afunction of the phase; Z ¼ =ðAÞ rad=mV=s measurements(dots) and Fourier fit with five harmonics (curve). (b) Optimal
waveform using (12) with  ¼ !, P ¼ 0:5, and the PRC in (a).

024102-3

PHYSICAL REVIEW LETTERS

PRL 111, 024102 (2013)1.50ln(∆φ−∆φf )1∆φ−∆φf (rad)

2

−1

1

−2

0.5

−3
0
0

10

20

30

−4

Time (s)2030

3Predicted
Experiment

0.3

2

κ

k (1/s)

10

Time (s)

0.4

0.2

1

0.1
0
0

0

10

Forcing Amplitude (mV)20

0week ending
12 JULY 2013

addition to the PRC itself. The methodology is promisingfor fast reestablishment of entrainment in oscillators thatintermittently break phase locking due to environmental orinternal effects, such as biological systems with fluctuations in chemical reaction rates due to the small number ofmolecules in a cell [29]. Finally, observe that our methodology is suitable for weak phase resetting, while strong
resetting requires control approaches that do not depend
on averaging but involve substantial changes to the state of
the oscillator.
This work was supported by the NSF Career Award
No. 0747877, the AFOSR YIP FA9550-10-1-0146, the
NSF Grant No. CHE-0955555, and the SLU PRF.

sine PRC optimal
Waveforms

FIG. 4 (color online). Electrodissolution experiments:(a) Phase difference ðtÞ for sine (dashed), PRC (thin), andoptimal (thick) forcing at A ¼ 12:5 mV. (b) Semilog plot ofðtÞ from (a). For t > 20 sec after forcing is applied, ðtÞ
decays exponentially to zero. Rates of entrainment (slope of thelinear fits): kðsinÞ ¼ 0:1093s1 , kðZÞ ¼ 0:167 s1 , kðoptimalÞ ¼
0:243 s1 . (c) Rate of entrainment as a function of forcing
amplitude for sin ( ), PRC (4), and optimal (h) forcing.
The slopes  in s1 mV1 of the fitted lines characterize theperformance of the waveforms for fast entrainment. (d)Normalized entrainment rates predicted from PRC estimates[left: ðsinÞ ¼ 0:68, ðPRCÞ ¼ 0:89, ðoptimalÞ ¼ 2:00] andmeasured experimentally [right: ðsinÞ ¼ 1:19, ðPRCÞ ¼
1:39, ðoptimalÞ ¼ 2:38] are highest for the optimal waveform.

the optimal waveform. This rate was found to be lower
for other waveforms such as sine and the PRC Z itself, as
shown in Figs. 4(a) and 4(b). To compensate for measurement errors and data processing inaccuracies, wemeasured the rate k at 7 amplitudes A between 2.5 and15 mV. The slopes  of the k vs A plots in Fig. 4(c)correspond to 0v ð’ Þ for normalized PRC and v, andare compared, along with values predicted using theestimated PRC for each experiment, in Fig. 4(d). The
optimal waveform performs significantly better.
The proposed technique for constructing optimal fast
entrainment waveforms can be applied to any nonlinear
oscillator, and requires no knowledge about its initial state.
Entrainment is achieved over the minimum number ofcycles possible for a given control energy such that phasemodel approximation and averaging remain valid. The
conditions required for such approximations to be appropriate for entrainment have been explored in previous work
[11,21]. When the initial slow phase of the system is far
from a stable fixed point, several cycles may be requiredfor convergence to the phase-locked state to be realized,
and this occurs least on average for the optimal waveform.
In contrast to previous studies on the control of oscillators
using phase models [7–10,20,21], the derivative of thephase response curve (PRC) plays an important role in*azlotnik@ese.wustl.edu†
ychen18@slu.edu‡
izkiss@slu.edu§
htan@synchro2.ee.uec.ac.jp∥
Corresponding author.
jsli@ese.wustl.edu[1] S. H. Strogatz, Nonlinear Dynamics and Chaos (Perseus
Books Group, Cambridge, MA, 2001), 1st ed.
[2] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: A Universal Concept in Nonlinear
Science (Cambridge University Press, New York, 2001).
[3] T. Naqvi and D. Winter, U.S. Patent No. 8 112 150 (2012).
[4] J. Waterhouse, T. Reilly, G. Atkinson, and B. Edwards,
Lancet 369, 1117 (2007).
[5] C. R. McClung, Adv. Genet. 74, 105 (2011).
[6] X. Feng, C. White, A. Hajimiri, and M. Roukes, Nat.
Nanotechnol. 3, 342 (2008).
[7] H. Moehlis, E. Brown, and H. Rabitz, J. Comp. Nonlin.
Dyn. 1, 358 (2006).
[8] I. Dasanayake and J.-S. Li, Phys. Rev. E 83, 061916 (2011).
[9] A. Zlotnik and J.-S. Li, in Proceedings of the 2011 ASME
Dynamic Systems and Control Conference, Arlington, VA,
2011 (ASME, New York, 2011), Vol. 1, pp. 479–484 [http://
proceedings.asmedigitalcollection.asme.org/proceeding.aspx?articleid=1638452].
[10] T. Harada, H.-A. Tanaka, M. J. Hankins, and I. Z. Kiss,
Phys. Rev. Lett. 105, 088301 (2010).
[11] A. Zlotnik and J.-S. Li, J. Neural Eng. 9, 046015 (2012).
[12] A. Granada and H. Herzel, PLoS One 4, e7057 (2009).
[13] S. Hata, K. Arai, R. F. Galán, and H. Nakao, Phys. Rev. E
84, 016229 (2011).
[14] B. Ermentrout, Neural Comput. 8, 979 (1996).
[15] S. H. Strogatz, Physica (Amsterdam) 143D, 1 (2000).
[16] F. Hoppensteadt and E. Izhikevich, Phys. Rev. Lett. 82,
2983 (1999).
[17] Y. Kuramoto, Chemical Oscillations, Waves, andTurbulence (Springer, New York, 1984).
[18] I. Z. Kiss, Y. Zhai, and J. L. Hudson, Science 296, 1676
(2002).
[19] S. Nakata, K. Miyazaki, S. Izuhara, H. Yamaoka, and
D. Tanaka, J. Phys. Chem. A 113, 6876 (2009).
[20] E. Brown, J. Moehlis, and P. Holmes, Neural Comput. 16,
673 (2004).

024102-4

PRL 111, 024102 (2013)

PHYSICAL REVIEW LETTERS[21] D. Efimov and T. Raissi, in Proceedings of the 8th IFAC
Symposium on Nonlinear Control Systems Bologna, 2010
(IFAC, Laxenberg, Austria, 2010), pp. 332–337 [http://
www.ifac-papersonline.net/Detailed/43863.html].
[22] W. Govaerts and B. Sautois, Neural Comput. 18, 817 (2006).
[23] B. Ermentrout, A Guide to XPPAUT for Researchers andStudents (SIAM, Philadelphia, 2002).
[24] I. Kornfeld, S. Fomin, and Y. Sinai, Ergodic Theory:
Differentiable Dynamical Systems, Grundlehren derMathematischen Wissenschaften Vol. 245 (SpringerVerlag, Berlin, 1982).

week ending
12 JULY 2013[25] F. Hoppensteadt and E. Izhikevich, Weakly Connected
Neural Networks (Springer-Verlag, New Jersey,1997).
[26] A. Hodgkin and A. Huxley, J. Physiol. 117, 500(1952) [http://www.ncbi.nlm.nih.gov/pubmed/12991237].
[27] O. Lev, A. Wolfberg, L. M. Pismen, and M. Sheintuch,
J. Phys. Chem. 93, 1661 (1989).
[28] I. Z. Kiss, Y. M. Zhai, and J. L. Hudson, Phys. Rev. Lett.
94, 248301 (2005).
[29] D. Andrieux and P. Gaspard, J. Chem. Phys. 128, 154506(2008).

024102-5