PRL 105, 088301 (2010)

PHYSICAL REVIEW LETTERS

week ending
20 AUGUST 2010

Optimal Waveform for the Entrainment of a Weakly Forced Oscillator
Takahiro Harada,1 Hisa-Aki Tanaka,2 Michael J. Hankins,3 and István Z. Kiss3,*
1

Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan
Department of Electronic Engineering, The University of Electro-Communications, Tokyo 182-8585, Japan3
Department of Chemistry, Saint Louis University, St. Louis, Missouri 63103, USA(Received 19 March 2010; published 16 August 2010)

2

A theory for obtaining a waveform for the effective entrainment of a weakly forced oscillator ispresented. Phase model analysis is combined with calculus of variation to derive a waveform with which
entrainment of an oscillator is achieved with a minimum power forcing signal. Optimal waveforms arecalculated from the phase response curve and a solution to a balancing condition. The theory is tested inchemical entrainment experiments in which oscillations close to and farther away from a Hopf bifurcation
exhibited sinusoidal and higher harmonic nontrivial optimal waveforms, respectively.
DOI: 10.1103/PhysRevLett.105.088301

PACS numbers: 82.40.Bj, 02.30.Xx, 05.45.Xt, 82.40.Np

Entrainment of oscillators to an external signal in nonlinear dissipative systems is a fundamental concept ofimportance in a large variety of applications [1]. Two
prominent examples include the time-scale adjustment of
the circadian system to light [2] and the cardiac system to apacemaker [3]. The entrainment process is rigorously described theoretically by phase or amplitude equations andcircle maps for weakly and strongly perturbed nonlinear
systems, respectively [1]. The general result of the theoretical analysis is that nonlinear oscillators can adjust theirfrequencies to that of the external source above a criticalforcing amplitude. In the forcing amplitude versus forcing
frequency diagram there are long vertical entrainment
regions called Arnold tongues.
A widely accepted tool for studying entrainment is thephase response curve (PRC) [2]. Specifically, the phase
response function (infinitesimal PRC) indicates the phaseshift of an oscillator due to an infinitesimal perturbation ofa system variable [4]. A classical problem in nonlinear
dynamics uses the phase response function and forcing
waveform with which all important features of the entrainment process (e.g., locking range, defined as the width of
the Arnold tongue at a given forcing amplitude) can beobtained for weakly perturbed systems [1].
Many applications require optimization of the entrainment process. This is often achieved by adjusting theforcing waveform to achieve a target entrainment feature.
A variety of control targets were explored: optimal input
was determined for establishing fast entrainments [5],circadian phase resetting [6,7], starting or stopping of theoscillations [7,8], and maximal resonance (energy transfer)
between the system and the forcing signal [9]. Control ofdeterministic [10] and stochastic [11,12] neuronal spiking
activity was achieved with a phase modeling approachcombined with variational methods to optimize spiking
time [10] and variance of firing rates [12].
In this Letter, we give a full account for the inverse of theclassical entrainment problem: namely, What is the minimal power forcing waveform that produces efficient en0031-9007=10=105(8)=088301(4)trainment of a limit-cycle oscillator in the weak forcing
limitation? Although the quality of entrainment could involve features such as stability and basin of attraction, here
we consider efficient entrainment as the occurrence of
maximum width (or minimum slopes) of the Arnold
tongue. This ‘‘locking range’’ quality marker has beencommonly used for injection-locked microintegrated oscillators as well as phase-locked loop circuits [13]. We
propose a versatile, efficient approach to obtain an exact
functional form of the optimal waveform provided that the
response function related to the forcing action had been
established. The theoretically obtained optimal waveforms, which exhibit some unexpected symmetry relationships with the response function, are tested in a simple
numerical model that include higher harmonics in the response function typically seen in strongly nonlinear oscillators. The experimental applicability of the method isdemonstrated with optimal entrainment of a chemical system, the oscillatory electrodissolution of nickel in sulfuricacid.
The entrainment process of a limit-cycle oscillator in the
weak forcing limit can be modeled by [14]dc
¼ ! þ Zð c Þfð!tÞ;(1)
dt
where c is the phase of the oscillator, Z is the phase
response function, and ! and ! are the natural frequencyof the oscillator and the forcing frequency, respectively. In
the weak forcing limit, Eq. (1) is further simplified byaveraging [15], asd!
¼ "! þ #ð!Þ;(2)
dt
where ! and "! are given by ! ¼ c % !t and "! ¼
! % !, respectively [16]. The interaction function #ð!Þ is
obtained from the forcing waveform f and the phase
response function Z, as #ð!Þ ¼ hZð" þ !Þfð"Þi, where "represents !t and h&i denotes
the average by " over its peH
riod 2#: h&i ' ð2#Þ%1 d". Entrainment occurs when thephase difference is locked, i.e., d!=dt ¼ "! þ #ð!Þ ¼ 0

088301-1! 2010 The American Physical Society

PRL 105, 088301 (2010)

PHYSICAL REVIEW LETTERS[7]. The range of frequency difference "!, where thesolution for stable steady state exists for !, defines thelocking range R½f) for a certain forcing waveform [17].
Therefore, the locking range is the difference between the
maximum (at ! ¼ !þ ) and minimum (at ! ¼ !% ) valuesof #ð!Þ where a phase-locked solution exists [18]. R½f) is
thus given by #ð!þ Þ % #ð!% Þ.
Now we are in position to formulate the optimal entrainment problem mathematically: the optimal forcing waveform (fopt ) maximizes the locking range R under certain
constraints. A convenient practical constraint is the total
power of the waveform over its period, hfð"Þ2 i. Therefore,
the optimal forcing waveform fopt gives maximal locking
range for a given (constant) forcing power P. We consider
this as a variational problem maximizing the functionalform
S ½f) ' R½f) % $ðhf2 i % PÞ;(3)where $ is the Lagrange multiplier. The solution to thevariational problem f* , a suitable candidate for the optimal
waveform, is obtained by ensuring that the first variation%S vanishes and the second variation %2 S is negative [19]:f* ð"Þ ¼ ð2$Þ%1 fZð" þ !þ Þ % Zð" þ !% Þg:(4)

The Lagrange multiplier can be obtained by substituting
the solution of Eq. (4) in the constant power constraint
pffiffiffiffiffiffiffiffiffiffi(hf*2 i % P ¼ 0): $ ¼ ð1=2Þ Q=P with Q ' hfZð" þ!þ Þ % Zð" þ !% Þg2 i. Note that f* in Eq. (4) has zeroaverage: hf* i ¼ 0.
To obtain f* of Eq. (4) the maximum (!þ ) and theminimum (!% ) ! values of # with the forcing waveform,
pffiffiffiffiffiffiffiffiffiffi#ð!Þ ¼ P=QhZð" þ !ÞfZð" þ !þ Þ % Zð" þ !% Þgi; (5)have to be determined. The conditions for the maximum
and minimum of # are as follows:#0 ð!+ Þ ¼ 0;

#00 ð!þ Þ < 0;and

#00 ð!% Þ > 0: (6)The first condition in Eq. (6), combined with Eq. (5), gives
hZ0 ð" þ !þ ÞZð" þ !% Þi ¼ hZ0 ð" þ "!ÞZð"Þi ¼ 0; (7)where "! ' !þ % !% . ("! is introduced to remove
phase ambiguity.) We shall refer to Eq. (7) as the balancing
condition because this equation realizes optimality bybalancing both terms in Eq. (4). The trivial "! ¼ 0 solution to Eq. (7) is discarded because it does not allow
entrainment [#ð!Þ ' 0 in Eq. (5)].
However, other solutions do exist in Eq. (7), because@hZ0 ð" þ "!ÞZð"Þi=@"!j"!¼0 ¼ %hZ0 ð"Þ2 i < 0,and0hZ ð" þ "!ÞZð"Þi is a periodic, bounded function of "!
in a large class of systems [1]. In particular, we have foundthat in models with twice differentiable, continuous Z a
solution with "! ¼ # exists; we call this solution and thecorresponding optimal waveform ‘‘generic’’ [20].
Here we test the above theoretical predictions by asimple model with the following response function:

Zð"Þ ¼ sin" þ a sinð2"Þ:week ending
20 AUGUST 2010

(8)This Z simulates the behavior of the Stuart-Landau oscillator [1] with a ¼ 0; therefore, we can consider a as ameasure of distance from Hopf bifurcation that can introduce higher harmonics in the response function. The balancing condition of Eq. (7) for this model is explicitly
written as ½1 þ 4a2 cosð"!Þ) sinð"!Þ ¼ 0. For jaj <
1=2, there is only one (nontrivial) solution: "! ¼ #.
The optimal waveform for
pffiffiffithis
ffiffiffi generic solution is independent of a: fopt ð"Þ ¼ % 2P sin". Although the response
function does contain a second order harmonic for 0 <jaj < 1=2, this term does not affect the shape of the optimal
waveform. This finding suggests that with Z containing
only weak second (or, in general, even) harmonics, thesinusoidal forcing is the optimal since the generic solution
to the balancing condition always exists. For example, forsystems close to Hopf bifurcation, which contain mostly
first and weak second (even) harmonics in Z, the optimal
waveform is sinusoidal. However, the odd harmonics do
appear in the generic optimal waveform, and thus systems
with relatively strong third (and higher odd) harmonics in
Z are expected to retain the odd harmonics in the optimal
waveform. For thepgenericsolution the locking range isffiffiffiffiffiffi
calculated as R ¼ 2P, which is again independent of a.
In the range jaj , 1=2, we have three different solutions
for the balancing condition and thus three candidates forthe optimal waveform. The generic solution with "! ¼ #exists; however, there appear two additional, ‘‘nongeneric’’solutions satisfying cos"! ¼ %1=4a2 . For these nongenericpsolutions
pffiffiffiffiffiffi ranges are identical: R ¼ ð1 þ
ffiffiffi pffiffiffi the locking
2
4a Þ P=ð2 2aÞ , 2P. This shows that in the range
jaj , 1=2, the nongeneric optimal waveforms outperformthe generic waveform; for large values of a, the improvement of locking range for the nongeneric over genericwaveforms increases approximately linearly. The nongeneric waveforms are not purely sinusoidal and depend on
the parameter a. Two nongeneric waveforms for a ¼ 0:95
are shown in Fig. 1(a).
We have also verified these theoretical predictions, byusing a standard genetic algorithm which numericallysearches for fopt . Figure 1(a) shows all fopt obtained by
the algorithm, for both cases of full [Eq. (1)] and averaged[Eq. (2)] phase models, with ! ¼ 10, a ¼ 0:95 (>1=2)
and P ¼ 0:5 or P ¼ 10:0. The numerical algorithm foundthe exact same optimal waveform as predicted by thetheory, for both cases up to around P ¼ 2:0. The numericalvalue of the locking range for these optimal solutions[maximum in R landscape in Fig. 1(b)] is also the sameas predicted by the theory; the nongeneric optimal solution
performs about 23.1% better than the simple sinusoidal(generic optimal) forcing, and it also performs about 89.8%
better than the pulse forcing in Fig. 2(c). Beyond P ¼ 2:0,
a small discrepancy appears in the optimal waveforms
obtained with the genetic algorithm. However, the shapeof the bimodal landscapes of R½f) in Fig. 1(b) is preserved
up to around P ¼ 10:0 for the case of Eq. (1), which

088301-2

PRL 105, 088301 (2010)

PHYSICAL REVIEW LETTERS

FIG. 1 (color online). Optimal forcing fopt obtained in themodel with the response function of Eq. (8). Panel (a) shows
all fopt , after rescaling for comparison, obtained by the genetic
algorithm for a ¼ 0:95: solid curves for P ¼ 0:5 with Eq. (2),dotted curves for P ¼ 10:0 with Eq. (1). In this genetic algorithm, fopt is searched among all functions of the form c1 sin " þ
c2 sinð2" þ c 2 Þ. Panel (b) shows the landscape of R for a ¼
0:95 (>1=2), with respect to (c1 , c 2 ). The peaks 1 and 2 of thelandscape correspond to the waveforms 1 and 2 in (a), respectively.

suggests that, at least in this particular example, beyond thestrictly weak forcing limit the theoretical prediction of
optimal waveform could be used as an initial candidatethat can be further optimized with other techniques.
The theoretical method for obtaining the optimal forcing
waveform was demonstrated in chemical experiments with
Ni electrodissolution. The experiments were carried out ina standard three-electrode electrochemical apparatus witha 1 mm diameter Ni wire as working, a 1.57 mm diameterPt coated Ti rod as the counter, and Hg=Hg2 SO4 =K2 SO4 ðsatÞ as the reference electrode. The potential of
the Ni wire was set to a circuit potential V ¼ V0 þ
Af0 ð"Þ, where V0 is a base potential, A is the amplitude
of forcing, and f0 ð"Þ is the normalized forcing waveformwith hf02 i ¼ 0:5. The current, proportional to dissolution
rate, was measured by the potentiostat. The frequency ofthe current oscillations was determined with the Hilberttransform method [1]. The experiments were carried out at
10 - C in 3 mol=L sulfuric acid solution. Further experimental details are given in previous studies [21,22]. Whena resistor of 1 k! is attached to the Ni electrode, oscillations arise through Hopf bifurcation at V0 . 1:07 V due tothe hidden negative differential resistance of the electrochemical system [23].
In accordance with previous experimental results [22], at
V0 ¼ 1:100 V, slightly above the bifurcation point, theoscillations are smooth with a PRC that contains predominantly first harmonic components [see Figs. 2(a) and 2(b)].
The experiments were carried out at fixed forcing frequencies and the amplitude of the waveform was increased until
entrainment occurs where the critical amplitude Ac wasrecorded. To compare the forcing capability of different
waveforms we determined
the reduced entrainabilities asqffiffiffiffiffiffiffiffi2
Er ¼ j"!j=ðAc hf0 iÞ. The average value of the reduced
entrainability for fixed positive and negative detunings is
proportional to the R value; therefore, it is a useful quantity
to characterize the entraining capability of a waveform. In

week ending
20 AUGUST 2010

FIG. 2 (color online). Experiments: Optimal entrainment ofsmooth oscillators, V0 ¼ 1:100 V. (a) Waveform. (b) Phase response curve (PRC). The PRC was obtained with pulsing theoscillations with an amplitude of %100 mV and width of 0.1 rad
and measuring the phase advance ($) in five independent experiments. The error bars indicate standard deviations in the fiveexperiments. (c) Three waveforms, sine-sawtooth-pulse, withwhich entrainments are tested. (d) Reduced entrainabilities of
the three waveforms: P, pulse; Saw, sawtooth; Sin, sine. The
detuning "!=! was set to +5%; at each detuning two experiments were carried out.

a typical set of experiments 5% detuning was applied, i.e.,"!=! ¼ +0:05; it was confirmed that smaller detuning(2%) gives Er values within 5% relative error. Experimentswith larger detunings resulted in large changes of oscillation waveforms and thus the weak forcing assumption of
the theory was not satisfied.
We have tested three waveforms with smooth oscillators; sine, sawtooth, and pulse waves are shown in
Fig. 2(c). The optimal waveform was determined by theabove algorithm [24]; only the generic waveform exists,
and because the PRC has very small third and higher harmonics (less than 1%), we can consider the sine waveform
as optimal within experimental error. The entrainabilities
in Fig. 2(d) show that, as predicted, the sine waveform
indeed outperforms the sawtooth and pulse waves.
When V0 is increased to 1200 mV, the oscillations
become moderately relaxational; the waveform and the
PRC exhibit higher harmonics [see Figs. 3(a) and 3(b)].
The analysis has revealed that three optimal waveforms
exist: one generic and two nongeneric, given in Figs. 3(c)–3(e), respectively. In Fig. 3(e) we also show a waveformthat is obtained by adding Gaussian random numbers withstandard deviation of 0.1 to the Fourier coefficients of themost optimal waveform; this waveform is referred to as
noisy optimal 2 waveform. The reduced entrainabilities
follow the theoretical predictions: sawtooth, sine, andpulse waves are inferior to the entraining power of the
optimal waveforms. As predicted, the generic optimal
waveform performs slightly worse than the nongeneric
optimal waveforms. The best waveform was optimal 2;

088301-3

PRL 105, 088301 (2010)

PHYSICAL REVIEW LETTERS

week ending
20 AUGUST 2010

I. Z. K. acknowledges financial support from Research
Corporation. H. T. thanks Dr. Y. Ando and N. Sagayama for
their help with genetic algorithms.

FIG. 3 (color online). Experiments: Optimal entrainment ofmoderately relaxational oscillators, V0 ¼ 1:200 V. (a) Oscillation waveform. (b) Phase response curve. The PRC was obtained
similarly to that in Fig. 2. (c)–(e) Optimal waveforms:(c) Generic optimal waveform (Og). (d) Optimal wave
form 1 (O1). (e) Thick solid line: optimal waveform 2 (O2).
Thin line: noisy optimal waveform 2 (O2n). (f) Reduced entrainabilities of the tested waveforms. The Er values are averages
of experiments with +2% and +5% detunings.

note that the shape of this waveform is strongly nontrivial
and by the addition of a small amount of noise to the
waveform the Er value decreases. Overall, the optimal
waveform obtained from the theory increased entrainabilities by about 50% to 90% relative to those obtained withsine-pulse-sawtooth waves.
Construction of an optimal waveform for entrainment
was proposed and tested here with a single oscillator. The
method, however, can be extended to a group of interacting
oscillators where effects related to the collective phase response function [25] shall be considered. The optimal signal can also be applied in closed-loop feedback systems
along with synchronization engineering [21] for seeking
optimal target dynamics. A limitation of the methodology
is the requirement for weak forcing so that phase models
can be applied. This limitation leads to an extension of themethod where the forcing signal is limited to small values.
These extensions, along with other targets that consider
stability and basin of attraction, will be considered in aforthcoming publication. The proposed methodology provides a framework for efficient design of entrainment applications in electrical circuit technology (e.g., forinjection-locked oscillators) as well as in biologicalpacemakers.

*Corresponding author.
izkiss@slu.edu[1] A. S. Pikovsky, M. G. Rosenblum, and J. Kurths, Synchronization—A Universal Concept in Nonlinear Sciences(Cambridge University Press, Cambridge, England, 2001).
[2] A. T. Winfree, The Geometry of Biological Time
(Springer-Verlag, New York, 1980).
[3] L. Glass, Chaos 1, 13 (1991); R. A. Gray, Chaos 12, 941(2002).
[4] G. B. Ermentrout, R. F. Galan, and N. N. Urban, Phys. Rev.
Lett. 99, 248103 (2007).
[5] A. E. Granada and H. Herzel, PLoS ONE 4, e7057 (2009).
[6] N. Bagheri, J. Stelling, and F. J. Doyle, PLoS Comput.
Biol. 4, e1000104 (2008).
[7] D. Forger and D. Paydarfar, J. Theor. Biol. 230, 521(2004).
[8] D. Lebiedz et al., Phys. Rev. Lett. 95, 108303 (2005).
[9] V. Gintautas and A. W. Hübler, Chaos 18, 033118 (2008).
[10] J. Moehlis, E. Shea-Brown, and H. Rabitz, J. Comp.
Nonlin. Dyn. 1, 358 (2006).
[11] J. Ritt, Phys. Rev. E 68, 041915 (2003).
[12] J. Feng and H. C. Tuckwell, Phys. Rev. Lett. 91, 018101(2003).
[13] R. E. Best, Phase-Locked Loops: Design, Simulation, andApplications (McGraw-Hill, New York, 1997).
[14] A. T. Winfree, J. Theor. Biol. 16, 15 (1967).
[15] Y. Kuramoto, Chemical Oscillators Waves and Turbulence(Dover, Mineola, New York, 2003).
[16] In this Letter we focus on the 1:1 resonant case. However,the presented theory is also valid for general m:n resonant
cases.
[17] S. H. Strogatz, Nonlinear Dynamics and Chaos (AddisonWesley, Reading, MA, 1994), p. 105.
[18] There may be multiple locking ranges for f. For graphicalexamples, see p. 64 of [15]. In our setting, the largest
range should be chosen for maximizing R in these situations.
[19] Here f are continuous and have continuous first derivatives, which belong to a function space with the norm k & kdefined by k f k¼ maxjfð"Þj þ maxjf0 ð"Þj. For details,
see I. M. Gelfand and S. V. Fomin, Calculus of
Variations (Dover, Mineola, New York, 2000), Chap. 1.
[20] The existence of the generic solution to balancing condition can be proven by partially integrating Z0 ð" þ #ÞZð"Þ,resulting in hZ0 ð" þ #ÞZð"Þi ¼ 0. Note that this solution
has Z2 symmetry.
[21] I. Z. Kiss et al., Science 316, 1886 (2007).
[22] I. Z. Kiss, Y. Zhai, and J. L. Hudson, Phys. Rev. Lett. 94,
248301 (2005).
[23] M. T. M. Koper, Adv. Chem. Phys. 92, 161 (1996).
[24] Up to three harmonics are considered in the PRC.
Addition of higher harmonics resulted in very minorchanges of the waveform.
[25] Y. Kawamura et al., Phys. Rev. Lett. 101, 024101 (2008).

088301-4