September 2015
EPL, 111 (2015) 50007
doi: 10.1209/0295-5075/111/50007

www.epljournal.org

Optimal synchronization of oscillatory chemical reactions
with complex pulse, square, and smooth waveforms signals
maximizes Tsallis entropy
Hisa-Aki Tanaka1 , Isao Nishikawa1 , Jürgen Kurths2 , Yifei Chen3 and István Z. Kiss3
1

The University of Electro-Communications - Tokyo 182-8585, JapanThe Potsdam Institute for Climate Impact Research - Telegraphenberg A31, 14473 Potsdam, Germany
3
Department of Chemistry, Saint Louis University - St. Louis, MO 63103, USA2

received 1 June 2015; accepted in ﬁnal form 24 August 2015
published online 18 September 201505.45.Xt – Synchronization; coupled oscillators
05.45.-a – Nonlinear dynamics and chaos
PACS 84.35.+i – Neural networks
PACS
PACS

Abstract – We show that the mathematical structure of Tsallis entropy underlies an importantand ubiquitous problem in nonlinear science related to an eﬃcient synchronization of weakly forced
nonlinear oscillators. The maximization of the locking range of oscillators with the use of phasemodels is analyzed with general constraints that encompass forcing waveform power, magnitude, orarea. The optimization problem is then recasted as a general form of Tsallis entropy maximization.
The solution of these optimization problems is shown to be a direct consequence from Hölder’sinequality. The resulting new maximization principle is conﬁrmed in numerical simulations andexperiments with chemical oscillations with nickel electrodissolution. While weakly nonlinear
oscillators have generic optimal waveforms (sinusoidal, 50% duty cycle square wave, and equally
paced bipolar pulses for power-, area-, and magnitude-constraints, respectively), strongly nonlinear
oscillators require more complex waveforms such as smooth, square, and pulse ones.
c EPLA, 2015
Copyright Introduction. – Synchronization (injection locking)of oscillators to an external forcing (injection signal) ofwidely diﬀerent waveforms is often used to provide means
to timing of essential system processes. In many branchesof science and engineering, methods for eﬃcient synchronization have been developed in recent years, because
the entrainment process plays a prominent role in thefunctioning of dynamical systems [1–3]. Diﬀerent ﬁelds
often use diﬀerent waveforms: phase locking in nerve
membrane dynamics with a sinusoidal current [4,5], circadian synchronization with light-dark cycles of a square
waveform [3,6], and reset of the timing in pacemakers andelectrical devices with pulse signals [7,8].
The various forcing waveforms called for establishinga theory for optimal synchronization of weakly forced
nonlinear oscillators; phase-model–based theories proved
to be useful because of their eﬀectiveness in describing
the synchronization process as well as their analyticaltractability [1,2,9–11]. The shape of the forcing waveform profoundly impacts the eﬃciency of synchronization for applications requiring low-power consumptions.

The width of the locking region requires waveforms related to a phase response curve (which shows the phaseshift in the oscillations due to an instantaneous pulse signal) [1]. Fast synchronization [2], or desynchronization [9]is related to the derivative of the phase response curve.
Diﬀerent forcing techniques using square or pulse injection
signals were developed for examples where the magnitude
and the area of the signals were minimized [3,6,7]. Most
previous theories resorted to variational calculus, such as
the Euler-Lagrange equation, which are intended inherently for local optimization and limited to smooth signals [1,2,9]. In contrast, recent progresses [10,11] haveshown that the optimal-synchronization problem for minimum power, magnitude, and area signals can be solved
for a global optimization even with non-smooth signals,
by means of Hölder’s inequality [12]:f g1 ≤ f α gβ ,

(1)

1


1
in which f α ≡ |f (θ)|α α denotes [ |f (θ)|α dθ] α , and−1−1
1 ≤ α, β ≤ ∞ and α + β = 1. f α amounts to the50007-p1

Hisa-Aki Tanaka et al.
Lα -norm of the function f if α > 0. For this inequality, thefollowing reverse form also holds (i.e., the reverse Hölder’sinequality):f g1 ≥ f α gβ ,(2)where 0 < α < 1 and α−1 + β −1 = 1. Note that eq. (2)11
is directly obtained by setting |f g| α |g|− α 1 to the l.h.s.
of Hölder’s inequality (1) [12]. In eqs. (1), (2), f g1 =f α gβ holds if and only if there exist constants A andB, not both 0, such that A|f (θ)|α = B|g(θ)|β for almost
all θ [12], which is referred to as the equality condition.
The general solution of the optimal synchronization implies that there might be a general mechanism that underlies the seemingly very diﬀerent optimization problems(for power-, magnitude-, or area-constraints, respectively).
In fact, entropy-based quantities often reﬂect the presenceof a maximization mechanism under given constraints.
Tsallis entropy [13,14] is one of such instances, which was
proposed as a versatile framework for expansion of therealm of classical (Boltzmann-Gibbs) entropies for nonlinear processes, in particular, those that exhibit power-law
behavior.
In this letter, as illustrated in ﬁg. 1, we show thatthe mathematical structure of Tsallis entropy underlies
an important and ubiquitous problem in nonlinear science
related to an eﬃcient synchronization of weakly forced
nonlinear oscillators. We formulate a uniﬁed theoretical framework for optimal synchronization and Tsallis
entropy maximization, i.e., a new maximization principle. A general analysis of this maximization solution is
performed for the characterization of optimal forcing signals for power-, magnitude-, and area-constraints. The
theory is tested with numerical modeling of an electrochemical oscillatory process. Finally, we demonstrate theeﬃciency of smooth, square, and pulse signals obtainedwith the theory for establishing optimal synchronization(with power-, magnitude-, and area-constraints, respectively) in laboratory experiments with oscillatory nickel
dissolution in sulfuric acid.

Fig. 1: Optimal forcings fopt, r maximizing synchronizability,
and the associated escort functions Popt, q maximizing Tsallis
entropy. (a) Arnold tongue determining the locking range. Δω
denotes the frequency diﬀerence between the oscillator and theforcing. (b) Phase response function Z for the Hodgkin-Huxleyneuron model [2]. (c), (d), (e): curves represent the optimal
forcings for r = 2, 1, ∞, respectively. (c) Power-reduced forcing, (d) area-reduced forcing, (e) magnitude-reduced forcing for
the Hodgkin-Huxley neuron model [2]. (f), (g), (h): escort functions Popt, q (x) of probability distributions maximizing Tsallis
entropy for q = 1, 3, 0, respectively. Note that for q = 1 Tsallis
entropy is identiﬁed with the Boltzmann-Gibbs (BG) entropy
since it converges to the BG entropy as q → 1 + 0 [13].

Phase locking with a weak forcing f (θ) occurs when theslowly evolving phase diﬀerence φ between the oscillator
and the forcing exhibits a stationary state, i.e.,dφ/dt = Δω + Γ(φ) = 0,(3)

where Δω = ω − Ω is the diﬀerence between the forcing frequency Ω and the natural frequency ω of the oscillations, and Γ(φ) is a periodic function deﬁned byπ1
Γ(φ) ≡ 2πZ(θ + φ)f (θ)dθ and Z(θ) is the phase re−πsponse (sensitivity) function (PRF) [16], obtained numerically or experimentally [17,18]. Note that this PRF is
the same notion as the phase response curve (PRC) under
weak inputs. The range of the frequency diﬀerence Δω between the oscillator and the forcing, in which the solution
for a stable steady state exists for φ, deﬁnes the locking
range (synchronizability) L[f (θ)] for a certain ﬁxed forcing
waveform f (θ). For maximizing L[f ] under the constraintμ1
2π f (θ) = 0, the functional J[f (θ)] ≡ L[f ] + 2π f (θ)
is πintroduced, where μ is a Lagrange multiplier and · =· dθ [1]. Moreover, J[f ] is rewritten as the following−π
inner product of f and g:

Tsallis entropy maximization underlies optimalsynchronization. – In order to make the essence of our
ﬁnding more transparent, we focus on the case of a single
oscillator with one external weak forcing; however, the results can be extended for designing optimal m : n mutualsynchronization [11] in networks of weakly coupled nonlinear oscillators. We consider a general class of periodic
functions f (θ) as the weak forcing, namely, those satisfying the following constraint: ||f ||r = M , in which both r
and M are positive constants. Especially for r = 1, 2, ∞,this constraint respectively deﬁnes the area-, power-, andmagnitude-reduced forcings for a certain small M . Beπ1f (θ)dθ = 0,sides f r = M , another constraint 2π−πi.e., a charge-balance constraint [1,2,15], is introduced, because total injection (injected charge) is required to be 0 wherein practical situations.
50007-p2

J[f (θ)] =1
f (θ)g(θ),2πg(θ) = Z(θ + Δφ) − Z(θ) + μ.

(4)

(5)

Optimal synchronization maximizes Tsallis entropy
μ and Δφ in eq. (5) are determined from Γ and the constraint on f r = M , as shown later. It should also be
noted that with the optimal waveform, Δφ becomes thediﬀerence between phase values of the maximum and minimum of Γ(φ). Thus, the optimal synchronizability can bedeﬁned by the following problem:
maximize J[f ],subject to f r = M.

(6)This optimal synchronization problem was previously
solved by using the calculus of variations, only forr = 2 [1,2]. In addition to the case of r = 2, this problem
was recently solved for general r, i.e., for 1 ≤ r ≤ ∞, via
Hölder’s inequality [10,11]. This ﬁnding hints at a new,complete analysis of Tsallis entropy maximization. Tsallis
entropy is deﬁned by the functional S[p(x)] of the probability density function (PDF) p(x),Optimization with Hölder’s inequality. – The features of the solution of the optimization problem (6)(or (11)) are now directly answered by Hölder’s inequality. Below, we outline the resulting solutions to both ofthe optimization problems (6) and (11), respectively forthe following cases: i) 1 < α < ∞, ii) α = ∞, iii) α = 1,
and iv) α < 1 of eqs. (1), (2). For illustration, ﬁg. 1 shows
the solutions for the optimization problem for a PRC for
the Hodgkin-Huxley neuron model along with the escort
functions of probability distributions maximizing the Tsallis entropy. Here, mathematical details related to obtaining solutions with Hölder’s inequality are omitted, as they
are presented in [10,11].
i) Case for 1 < α < ∞: The optimal forcing fopt, r
is obtained by setting α = r in the equality condition ofeq. (1),1(7)S[p(x)] ≡ [1 − p(x)q x ]/(q − 1),r ) r−1 ,(12)fopt, r (θ) = M sgn[g(θ)](|g(θ)|/g r−1
∞where ·x = −∞ · dx [13]. In addition, the escort function of p(x) is also deﬁned as P (x) ≡ p(x)q /p(x)q x . For where g is given as in eq. (5) and this representation is
the parameter q, we usually assume 0 < q < 3 and q = 1. proved to be unique, i.e., there is no other optimal forcing
without that in eq. (12) [10]. The two parameters Δφ andTsallis entropy maximization is then deﬁned as
μ in g(θ) of eq. (5) are determined by solving two nonlinear
equations of Δφ and μ, which are obtained by plugging
maximize S[p(x)],subject to p(x)x = 1 and h(x)P (x)x = σ 2 , (8) eq. (12) into the constraints f  = 0 and f r = M [10,11].
This solution (12) results in a relatively smooth forcing
where we assume h(x) = x2 which constrains the second signal; the shape of the signal often resembles that of the
moment of P (x) to be σ 2 .
PRC as far as r takes a moderate value. As a special case,
Note for a given q and σ, eq. (8) is equivalent to the the power-reduced forcing (in ﬁg. 1(c)) is obtained fromminimization of p(x)q x under the constraint (h(x) − eq. (12) for r = 2.
σ 2 )pq (x)x = 0, i.e., h(x)P (x)x = σ 2 in eq. (8). ThereSimilarly to the optimal forcing fopt, r , the optimal PDF
fore, by introducing the following functional T [p(x)] and to Tsallis entropy (7) is obtained by setting α = q −1 (i.e.,the multiplier λ, Tsallis entropy maximization is recast as 0 < q < 1) in the equality condition of eq. (1),the equivalent optimization:11
popt, q (x) = γ q |λh(x) + (1 − λ)σ 2 | 1−q ,(13)minimize (or maximize)
T [p(x)] ≡ σ 2 p(x)q x + λ(h(x) − σ 2 )p(x)q x
where γ and λ are determined from the constraints in
eq. (8) and a cut-oﬀ is naturally introduced, i.e., pq (x) ≡ 0= p(x)q [λh(x) + (1 − λ)σ 2 ]x
2
for q > 1 (or q < 1), subject to p(x) = 1. (9) for x satisfying λh(x) + (1 − λ)σ < 0. As a result,

x


 1
1
1 − q x2 1−qpopt, q (x) =1−Zq
3 − q σ2 +This T [p] of eq. (9) is directly mapped to Hölder’s inequality in eq. (1) or eq. (2) by inserting
p(x)q = f and λh(x) + (1 − λ)σ 2 = g,

(10)

(14)
1

2−q 1
2 2
in [13,14] is recovered, where Zq = ( 3−q1−q σ ) B( 1−q , 2 )to T [p]. Namely, Tsallis entropy maximization is equiva- and [a] ≡ max(a, 0).
+
lent to the optimization:ii) Case for α = ∞: Similarly to the case of 1 < α < ∞,Hölder’s inequality (1) results in the optimal forcing by
maximize (or minimize) f g for q < 1 (or q > 1),
setting α = r = ∞,(11)
subject to f  q1 = 1,

in which f g = p(x)q [λh(x) + (1 − λ)σ 2 ]x and f  q1 =
pq  q1 = p(x)x = 1. From eqs. (4), (6) and eq. (11), we
realize that the mathematical structure of Tsallis entropy
maximization underlies the optimal synchronization. We
note that Rényi entropy maximization [13] is recast into
the same form of eq. (11) if the constraints in eq. (8)are imposed, for instance.

fopt, ∞ (θ) = M sgn[g(θ)],(15)

which is consistent with the limit of eq. (12) as r → ∞and its uniqueness is proved [10]. Δφ and λ in g(θ) ofeq. (15) are determined, similarly to the above case for1 < α < ∞. Note that this forcing shown in ﬁg. 1(e)corresponds to the bang-bang control obtained under themagnitude-reduced constraint [15].

50007-p3

Hisa-Aki Tanaka et al.
As for Tsallis entropy, α = ∞ (i.e., q = α−1 = 0) implies p(x)q ≡ 1 for p(x) > 0 in eq. (7). Hence, Tsallis
entropy√
is maximized
for any p(x) with the compact sup√port [− 3σ, 3σ] although a detailed argument for thisis omitted here. Also from eq. (13) it is proved that theescort function of popt, q (x) converges to the step function
Popt, 0 (x) shown in ﬁg. 1(h) (pointwise, q → 0 + 0).
iii) Case for α = 1: In contrast to the case ii), fromeq. (1) with α = 1, fopt, 1 is proved to become a pair consisting of one arbitrarily tall negative pulse and one arbitrarily tall positive pulse (i.e., bipolar pulses) separated
by Δφmax [10] as shown in ﬁg. 1(d). Since this optimal
forcing for r = 1 includes only one parameter Δφmax , theoptimal Δφmax is determined algorithmically from Z, simply by tuning Δφmax in the range of [0, π] and observing
the resulting locking range L[f ] for each Δφmax numerically or experimentally.
In contrast to the case of q = 0, Tsallis entropy (7) is not
deﬁned at q = 1. However, it is known that it converges
to the Gaussian distribution as q → 1 + 0 [13].
iv) Case for α < 1: As opposed to the above cases of 1 ≤r (= α) ≤ ∞, it is not known that fopt, r has any physicalinterpretation for r < 1. However, for Tsallis entropy, thecase of q = α−1 > 1 is important from its applications,since it captures the power-law distributions with heavy
tails. Similarly to the case of 0 < q < 1, the optimal
PDF popt, q (x) is uniquely obtained as the same formulain eq. (14), by setting α = q −1 (i.e., q > 1) in the equality
condition of eq. (2). Note that Zq is obtained as Zq =
3−q
1
2 12
( 3−q
q−1 σ ) B( 2(q−1) , 2 ) and q < 3 is naturally required.
Simulations. – The principal predictions of the theory
are ﬁrst conﬁrmed in numerical simulations of a chemicalcorrosion model (see [19]) describing the electrodissolutionof nickel in sulfuric acid. The dimensionless model simulates the variation of the electrode potential E and surface
coverage Θ of oxide species with time t at a constant circuit potential V and cell resistance S,
ch exp(E/2)
V −E
dE
=
−
+ a exp(E) (1 − Θ) ,dt
S
1 + ch exp(E)
exp(E/2)b ch exp(2E)dΘ=(1 − Θ) −Θ. (16)Gdt
1 + ch exp(E)c ch + exp(E)The model parameters were chosen to produce limit cycle oscillations close to a supercritical Hopf bifurcation,
namely, the dimensionless circuit potential V is set at 15,
acid concentration ch is set at 1600, and kinetic parameters are a = 0.3, b = 6 × 10−5 , c = 0.001, S = 20, and
G = 0.01. The PRF Z for perturbation of the electrodepotential is shown in ﬁg. 2(a). The PRF was applied to obtain three optimal waveforms for the power-, magnitude-,and area-reduced constraints as shown in ﬁg. 2(b). For
power-constraint, the waveform coincides with the generic
optimal waveform Z(θ + π) − Z(θ) [1], and has a relatively
smooth character that retains the odd harmonic components of the PRC. The waveforms for the magnitude- andFig. 2: (Color online) Simulation results from optimal synchronization of a chemical oscillatory model (16). (a) PRF.
(b) Three waveforms: Sq: square wave (thick dashedline). P: pulse wave (thick solid line). SO: smooth (powerreduced) optimal waveform determined as Z(θ + π) − Z(θ).
(c), (e), (g): synchronizabilities Er (P ), Er (M ), Er (A) for
power-, magnitude-, and area-reduced constraint, respectively. For comparison, results of sinusoidal synchronization
simulations are also shown (waveform S). (d), (f), (h): synchronizability histogram for 1000 random waveforms for power-,magnitude-, and area-reduced constraint, respectively. The
random waveforms were generated with uniformly distributed
random numbers for the Fourier coeﬃcients of the waveforms
up to ﬁve harmonics. The arrow denotes the maximum synchronizability limit obtained by smooth optimal fopt, 2 , squarefopt, ∞ , and pulse fopt, 1 , respectively.

area-constraint are also the simplest ones expected for
50% duty cycle square wave and two evenly paced bipolar
pulses, respectively. The relatively simple shapes of the
optimal waveforms imply that under these conditions theoscillations have relatively weak nonlinear characters.
To compare the synchronizing capability of diﬀerent
waveforms, we deﬁned the reduced synchronizabilities Er
as follows. Numerical simulations for the synchronizationof the oscillations were carried out by periodically varying the circuit potential variable as V = 15 + Af (Ωt).
The simulations were conducted at ﬁxed forcing frequencies Ω and the amplitude A of the waveform was increased
until synchronization occurs where the critical amplitude
Ac was recorded. The reduced synchronizabilities Er aredeﬁned as
Er = |Δω|/Ac R,(17)

where R is the rescaling factor that enables a comparison
1of diﬀerent signals for their power (R = f 2 = f 2  2 ),magnitude (R = f ∞ = ess. sup. (|f |)), or area (R =f 1 = |f |). The average value of the reduced synchronizability for ﬁxed positive and negative detunings is
proportional to the width of the Arnold tongue at a given
value of the constraint; therefore, it is a useful quantity
to characterize the synchronizing capability of a waveform. In the following, 2% detuning was applied, namelyΔω/ω = ±0.02.
The simulations conﬁrmed that for power-, magnitude-,and area-reduced synchronizabilities the smooth optimal,
square, and pulse waveforms generate the widest locking
range L[f ] (see ﬁgs. 2(c), (e), and (g); for comparison,

50007-p4

Optimal synchronization maximizes Tsallis entropy

Fig. 3: (Color online) Experimental results from optimal synchronization control of oscillations close to Hopf bifurcation(Vs = 1.080 V) with three chosen waveforms. (a) The experimentally measured nearly harmonic PRF of the oscillations
and the shape of oscillations (inset). (b) The three considered control waveforms: sine (red), pulse (green) and square(blue). (c), (d), (e): comparison of the reduced synchronizability Er (P ), Er (A), Er (M ) for three constraints: power-, area-,and magnitude-reduced constraint, respectively. The entrainabilities are linearly renormalized such that the value for thesinusoidal waveform is 100. The same colors are used for waveforms as in panel (b).

results for synchronization with sinusoidal waveform arealso presented). As a further conﬁrmation of the nontrivial shapes of the optimal waveforms, we also performed
corresponding numerical simulations with 1000 waveforms
with random Fourier coeﬃcients. The histogram of thesynchronizabilities (see ﬁgs. 2(d), (f), and (h)), showthat these random waveforms have bell-shape distributions whose range is below the synchronizability limitedby the corresponding optimal waveform.
Experiments. – The enhancement of synchronization
with optimized forcing waveforms was veriﬁed in experiments with an electrochemical oscillating reaction. The
experimental system was a standard three-electrode electrochemical cell (controlled with a Gamry Reference 600
potentiostat) with a 1 mm diameter Ni-wire working electrode in 3 M sulfuric acid electrolyte, a 1.5 mm Pt-coatedTi-rod counterelectrode, and a Hg/Hg2 SO4 /K2 SO4 (sat)reference electrode. (All potentials are given with respect
to the reference electrode.) The cell was kept at a temperature of 10 ◦ C. When a certain ﬁxed circuit potential
V = Vs is applied with a 1 kOhm external resistance attached in series with the Ni wire, oscillatory behavior canbe observed that occurs through a Hopf bifurcation at
Vs = 1.06 V with the period T ∼ 2.2 s. The forcing signal
waveform f (Ωt) (= f (θ)) is superimposed on the circuit
potential such that V = Vs + Af (θ), where A is the amplitude. The output signals were collected and analyzed
using a National Instruments Labview data acquisition
system at 100 Hz and 1000 Hz sample acquisition rate forthe synchronization experiments and PRF measurements,

Fig. 4: (Color online) Optimal synchronization control of experimental relaxation oscillations far away from Hopf bifurcation (Vs = 1.180 V). (a) The experimentally measured PRF
and the waveform (inset). (b) The three optimal control waveforms corresponding to the three examples: optimal wave (red)for power, optimal thin pulse (green, pulse width: 0.011T )for area, and double square (blue) for magnitude constraints.
(c), (d), (e): comparison of the reduced synchronizability (Er )for all the six controls in three constraints: power-, area-, andmagnitude-reduced constraint, respectively. The entrainabilities are linearly renormalized such that the value for the sinusoidal waveform is always 100. The same colors are used forwaveforms as in (b), with other waveforms in black. (Symbols: S for sine, O for power optimal, S1 for 0.5T -separated
square shown in ﬁg. 3(b), S2 for optimal double square, P1 for
0.5T -separated thin pulse, P2 for optimal thin pulse).

respectively. The phase of the oscillations was determinedwith the Hilbert transform method [18].
At a ﬁxed circuit potential Vs = 1.080 V, the system
is close to Hopf bifurcation and the oscillation waveform
and the phase response function are nearly harmonic asshown in ﬁg. 3(a). For such nearly harmonic waveform
PRFs, the optimal waveforms are smooth sinusoidal, twoevenly paced pulses, or a 50% duty cycle square wave forthe power-, area-, and magnitude-reduced constraint, respectively. (Here, the pulse wave has a width of 0.1T .)We measured reduced synchronizabilities Er of these three
waveforms at each constraint, and compared them separately (ﬁgs. 3(c), (d) and (e)). As predicted by the theory,the smooth sinusoidal, pulse, and square waves are superior to the others for power-, area-, or magnitude-reduced
constraint, respectively. Even for this relatively harmonic
PRF, the optimal wave can outperform the others by 6%
to 93% depending on the speciﬁc constraint and waveform.
The chemical oscillation waveform and PRF include
strong higher harmonics when the electrodes are polarized at an elevated potentials of Vs = 1.180 V as shown inﬁg. 4(a). For this system, we also used the theory to determine the optimal waveforms for the three constraints asshown in ﬁg. 4(b). For the minimum-power constraint, theoptimal waveform is a smooth signal whose shape is connected to that of the PRF in a highly nontrivial manner.
For the area constraint, the two pulses are separated by a

50007-p5

Hisa-Aki Tanaka et al.
nontrivial time of 0.25T instead of the usual 0.5T for the
harmonic case. For the magnitude constraint, the optimalsquare waveform has a double-square shape; each squareshape has a diﬀerent size and the two square waves are
separated by 0.38T . We measured the synchronizability
for the three optimal waves as well as those obtained forthe harmonic PRFs. This comparison is expected to reveal the importance of higher harmonics in the PRF forthe experimental system. The experimental result for the
power constraint (ﬁg. 4(c)) shows that the optimal smooth
waveform clearly outperforms the other tested waveforms;
for example, the synchronizability is enhanced by 50%
compared to the smooth sinusoidual waveform. For thearea-reduced constraint (ﬁg. 4(e)) the situation is similar;
the optimal pulse waveform is 19% better than the 0.5T
separated pulse, and over two times better than a simple
sine wave. These results imply that the sinusoidal waveform can be very poor for the power- and area-reduced
constraint cases. Finally, for the magnitude-reduced forcing (ﬁg. 4(d)), the optimal wave achieves about the same
performance as a 50% duty cycle square wave does. However, in this example we can still see some eﬀect of higher
harmonics; the diﬀerence between the square wave andthe sine wave is greatly enhanced (60%) in contrast to therelatively smaller enhancement (25%) observed with thenearly harmonic PRF.
Conclusion. – We showed that the optimal synchronization problem can be translated into Tsallis entropy
maximization, which can be solved by Hölder’s inequality.
This mathematical ﬁnding provides an important physical realization of the widely known and commonly usedinequality. The numerical and experimental results in achemical system conﬁrmed the eﬃciency of Tsallis entropy
maximization analogy design of forcing waveforms for optimal synchronization. Because the theory in the present
paper is independent of the details of a system, it could
also apply to other optimal designs such as those in practical electrical circuits [20], in micro- and nano-integrated
circuit oscillators [8] as well as in biology for PRF-based
synchronization of neuronal spiking or circadian rhythm.
Further theoretical treatment is required for noisy environment, oscillator networks, or for strongly forced oscillators [21]; such attempts could be possible with thedevelopment of the phase model theory.
∗∗∗
This work was supported by MEXT (No. 26286086) andby the Telecommunications Advancement Foundation in

Japan. IZK acknowledges support from National Science
Foundation CHE-0955555 grant. The authors are gratefulto Prof. Constantino Tsallis for fruitful discussions
and to Prof. Hiroki Suyari for his critical reading of
our manuscript. The authors would also like to thank thereviewers for providing valuable comments.
REFERENCES[1] Harada T., Tanaka H.-A., Hankins M. J. and Kiss
I. Z., Phys. Rev. Lett., 105 (2010) 088301.
[2] Zlotnik A., Chen Y., Kiss I. Z., Tanaka H.-A. and
Li J. S., Phys. Rev. Lett., 111 (2013) 024102.
[3] Bagheri N., Stelling J. and Doyle F. J., PLoS
Comput. Biol., 4 (2008) e1000104.
[4] Aihara K., Matsumoto G. and Ikegaya Y., J. Theor.
Biol., 109 (1984) 249.
[5] Jackson J. C., Windmill J. F. C., Pook V. G. andRobert D., Proc. Natl. Acad. Sci. U.S.A., 106 (2009)
10177.
[6] Fukuda H., Murase H. and Tokuda I. T., Sci. Rep.,
3 (2013) 1533.
[7] Luther S. et al., Nature, 475 (2011) 235.
[8] Takano K., Motoyoshi M. and Fujishima M., Proceedings of the IEEE Asian Solid-State Circuits Conference,
Korea., 2007 (IEEE) 2007, p. 336.
[9] Wilson D. and Moehlis J., SIAM J. Appl. Dyn. Syst.,
13 (2014) 276.
[10] Tanaka H.-A., Physica D: Nonlinear Phenom., 288(2014) 1.
[11] Tanaka H.-A., J. Phys. A: Math. Theor., 47 (2014)
402002.
[12] Hardy G., Littlewood J. E. and Pólya G.,
Inequalities, second edition (Cambridge Mathematical
Library) 1988.
[13] Tsallis C., Introduction to Nonextensive Statistical
Mechanics (Springer, New York) 2009.
[14] Prato D. and Tsallis C., Phys. Rev. E, 60 (1999) 2398.
[15] Nabi A. and Moehlis J., J. Math. Biol., 64 (2012) 981.
[16] Kuramoto Y., Chemical Oscillations, Waves andTurbulence (Springer, Berlin) 1984.
[17] Winfree A. T., The Geometry of Biological Time, 2ndedition (Springer, New York) 2001.
[18] Pikovsky A. S., Rosenblum M. G. and KurthsJ., Synchronization: A Universal Concept in Nonlinear
Sciences (Cambridge University Press, Cambridge) 2001.
[19] Pyragas K., Pyragas V., Kiss I. Z. and Hudson J. L.,
Phys. Rev. Lett., 89 (2002) 244103.
[20] Nagashima T., Wei X., Tanaka H.-A. and Sekiya H.,
IEEE Trans. Circuits Syst., 61 (2014) 2904.
[21] Kurebayashi W., Shirasaka S. and Nakao H., Phys.
Rev. Lett., 111 (2013) 214101.

50007-p6