July 2014
EPL, 107 (2014) 10009
doi: 10.1209/0295-5075/107/10009

www.epljournal.org

Design and control of noise-induced synchronization patterns
W. Kurebayashi1(a) , T. Ishii1 , M. Hasegawa2 and H. Nakao1
1Graduate School of Information Science and Engineering, Tokyo Institute of Technology - 2-12-1 O-okayama,
Meguro-ku, Tokyo 152-8552, Japan2
Department of Electrical Engineering, Tokyo University of Science - 6-3-1 Niijuku, Katsushika-ku,
Tokyo 125-8585, Japan
received 4 March 2014; accepted in final form 10 June 2014
published online 7 July 201405.45.Xt – Synchronization; coupled oscillators
02.50.Ey – Stochastic processes
PACS 05.40.Ca – Noise
PACS

PACS

Abstract – We propose a method for controlling synchronization patterns of limit-cycle oscillators
by common noisy inputs, i.e., by utilizing noise-induced synchronization. Various synchronization
patterns, including fully synchronized and clustered states, can be realized by using linear filtersthat generate appropriate common noisy signals from given noise. The optimal linear filter canbe determined from the linear phase response property of the oscillators and the power spectrumof the given noise. The validity of the proposed method is confirmed by numerical simulations.
c EPLA, 2014
Copyright ⃝Introduction. – Various nonlinear dynamical systems
tend to synchronize when driven by a common noisy
input. This phenomenon, called noise-induced synchronization, is observed in many systems, for example, inneurons [1–4], electric circuits [5], electronic devices [6],microbial cells [7], lasers [8], and chaotic dynamical systems [9,10]. It has been clarified that noise-induced synchronization has quite a diﬀerent mechanism from phaselocking to periodic forcing, i.e., the oscillators are not entrained by the input but still exhibit mutual synchronization, characterized by coherent distributions of the phase
diﬀerences. Analytical investigations of this phenomenon
for limit-cycle oscillators can be performed by using thephase reduction method [11], and its properties have been
widely studied in the last decade [12–24].
In contrast to phase locking that requires periodic forcing whose frequency is close to rational multiples of thenatural frequency of the oscillators, noise-induced synchronization can occur even for white noise [12–18,20–24].
Thus, noise-induced synchronization may more easily berealized than phase locking, because common environmental noise is ubiquitous in nature. Indeed, it is conjectured
that some plants utilize common environmental noise forsynchronization to realize biological functions [25,26]. Insuch biological systems, some kind of filtering mechanisms
for the environmental noise may also exist to improvenoise-induced synchronization. Such filtering mechanisms,(a) E-mail: kurebayashi.w.aa@m.titech.ac.jp

if any, would also be useful in practical applications suchas noise-induced synchronization of sensor networks [27].
The aim of this letter is to provide a theoretical basis for
optimizing noise-induced synchronization by filtering the
given noise.
In our previous works [16,21], we developed quantitative theories that predict global statistical properties of
the noise-induced synchronization, such as formation ofvarious synchronization patterns and fluctuations aroundthese patterns, not only whether the oscillators synchronize with the others. The synchronization patterns arecharacterized by the probability density function (PDF)of the phase diﬀerences between the oscillators, which canbe calculated from the phase response property of the oscillators and the statistical properties of the noisy inputs.
In this letter, on the basis of ref. [21], we propose
a method for designing and controlling various noiseinduced synchronization patterns of limit-cycle oscillators, including the fully synchronized and clustered states.
Since the synchronization pattern of the oscillators depends on the statistical properties of the noisy input,
we can design the synchronization pattern by optimizing
the noisy input so that some objective function, e.g., thedegree of synchronization, is maximized. We develop an
optimization method for the noisy input, namely, for a linear filter that transforms given noise into an appropriate
noisy input so that the desired synchronization pattern isrealized. The validity of the proposed method is confirmed
by numerical simulations.

10009-p1

W. Kurebayashi et al.
Model. – We consider an ensemble of N uncoupled to a one-dimensional phase equation for small ϵ,
identical limit-cycle oscillators subjected to correlatedθ̇j = ω + ϵZ(θj )[f ∗ (ξj (t) + ηj (t)) + ζj (t)]noise and independent noise, described by the following
Langevin equations:(2)+ ϵ2 ν(θj ) + O(ϵ3 ),(1)Ẋj (t) = F (Xj ) + ϵG(Xj )[Ij (t) + ζj (t)],where θj (t) ∈ [0, 2π) is the phase of the oscillator j,
Z(θj ) is a sensitivity function that characterizes the refor j = 1, . . . , N . Here, Xj (t) ∈ Rn is the state ofsponseof the oscillator phase to noisy inputs, and ν(θj )the oscillator j at time t, F (Xj ) ∈ Rn is a vector field
represents
the eﬀect of amplitude relaxation dynamics ofrepresenting the oscillator dynamics, G(Xj ) ∈ Rn represtochastic limit-cycle oscillators [20,28] (this term evensents the coupling of the oscillator to the noisy inputs,
tually vanishes and does not play a role in the followIj (t) ∈ R is the correlated noise, ζj (t) ∈ R is the indepening argument). The sensitivity function Z(θ) is given as
dent noise that drives each oscillator independently, andZ(θ) = G⊤ (X)∇X θ(X)|X=X0 (θ/ω) , where θ(X) is theϵ (≪ 1) is a small parameter that controls the strength
isochron of the limit cycle and ∇X θ(X)|X=X0 (θ/ω) repreof the noisy inputs. We assume that eq. (1) possesses
sents its gradient on the limit-cycle orbit at phase θ [11].
a stable limit-cycle orbit X0 (t) with period T and frequency ω := 2π/T when ϵ = 0. The correlated noise isCharacterization of synchronization patterns. –generated from two given noisy signals by a linear filter As discussed in refs. [16,21], the phase diﬀerence between
as Ij (t) = f ∗ (ξ(t) + ηj (t)), where the noise ξ(t) ∈ R is two oscillators φj,k := θj − θk characterizes the noisecommon to all oscillators, ηj (t) ∈ R is independently ap- induced synchronized state. Since statistical properties
plied to each oscillator, f (τ ) ∈ R is a filter function that of ηj (t) and ζj (t) do not depend on the oscillator index j,transforms the given noise to appropriate noise for real- the PDF of the phase diﬀerence φj,k does not depend on
izing desired synchronization patterns, and the star (∗) the indices j and k. Thus, in the following, we denote the! +∞denotes convolution f ∗ α(t) = −∞ f (τ )α(t − τ )dτ . We phase diﬀerence by φ without the oscillator indices.
In our previous work [21], we obtained the stationary
introduced two independent noise terms ηj (t) and ζj (t) to
take into account the eﬀect of external disturbances before PDF U (φ) of the phase diﬀerence φ by employing eﬀective white-noise approximation of the phase equations (2)and after filtering.
For example, the oscillators described by eq. (1) can be subjected to correlated colored noise and by deriving an
regarded as spiking neurons receiving artificial injection averaged Fokker-Planck equation for φ from the multivaricurrents as in ref. [3]. In this case, the filtered noise Ij (t) ate Fokker-Planck equation for the phase variables {θj }.
represents the injected current to each neuron, and we may It turns out that the correlation functions of the noise play
suppose ηj (t) = 0. The filter f (τ ) is used for generating an important role, and U (φ) is explicitly given by
an appropriate injection current, and ζj (t) is independent
1
1
,
(3)
U (φ) = ·noise inherent in each neuron, e.g., channel or synapticū g(0) − g(φ) + h(0)noise. In the situation of ref. [27] where noise-induced
synchronization of wireless sensor networks is considered, where! +π ū ∈ R is a normalization constant determined byeach oscillator described by eq. (1) corresponds to each −π U (φ)dφ = 1, and g(φ) ∈ R and h(φ) ∈ R are correof the noise terms in eq. (2), defined as
sensor node. The sensor node measures a noisy environ- lation functions! +∞mental signal ξ(t) + ηj (t), and the filter f (τ ) implemented g(φ) = −∞ ⟨Z(θ(t))f ∗ ξ(t)Z(θ(t − τ ) + φ)f ∗ ξ(t − τ )⟩dτ!
on each sensor node transforms the signal into an appro- and h(φ) = +∞ ⟨Z(θ(t))[f ∗ ηj (t) + ζj (t)]Z(θ(t − τ ) +−∞priate noisy input Ij (t) that induces synchronization of φ)[f ∗ ηj (t) + ζj (t − τ )]⟩dτ . In Fourier representation,
the sensor nodes. We may suppose ζj (t) = 0 in this case. these functions can be written as g(φ) = "+∞ g eiℓφℓ=−∞ ℓ"+∞
We assume that ξ(t), ηj (t), and ζj (t) are mutually indeand h(φ) = ℓ=−∞ hℓ eiℓφ , where the Fourier coeﬃcientspendent zero-mean Gaussian noise, i.e., ⟨ξ(t)⟩ = ⟨ηj (t)⟩ =g , h ∈ R are given by⟨ζj (t)⟩ = 0 and ⟨ξ(t)ηj (t − τ )⟩ = ⟨ξ(t)ζj (t − τ )⟩ = ℓ ℓ⟨ηj (t)ηk (t−τ )⟩ = ⟨ζj (t)ζk (t−τ )⟩ = ⟨ηk (t)ζℓ (t−τ )⟩ = 0 forgℓ = |zℓ |2 |A(ℓω)|2 Pξ (ℓω),
any j, k, and ℓ (j ̸= k), where ⟨·⟩ denotes the ensemble avhℓ = |zℓ |2 |A(ℓω)|2 Pη (ℓω) + |zℓ |2 Pζ (ℓω),(4)
erage over realizations of ξ(t), ηj (t), and ζj (t). For simplic!
+π1ity, we assume that the statistical properties of ηj (t) and and zℓ := 2πe−iℓθ Z(θ)dθ is the Fourier coeﬃcient of−πζj (t) do not depend on the oscillator index j. Their power Z(θ). From eq. (3), we see that the PDF U (φ) is sym! +∞spectra are given by Pξ (Ω) := −∞ e−iΩτ ⟨ξ(t)ξ(t − τ )⟩dτ , metric about φ = 0 and has a maximum at φ = 0. For! +∞Pη (Ω) := −∞ e−iΩτ ⟨ηj (t)ηj (t − τ )⟩dτ , and Pζ (Ω) := example, a PDF U (φ) with a single peak at φ = 0 represents the synchronized state of the oscillators, and U (φ)! +∞ −iΩτe⟨ζj (t)ζj (t−τ )⟩dτ . We also define the amplitude with k peaks represents the k-clustered state. Thus, when−∞! +∞ −iΩτ
response of the filter f (τ ) as A(Ω) := | −∞ ef (τ )dτ |. the amplitude response A(Ω) of the filter f (τ ) enhances
By the phase reduction method [11], we can reduce the the k-th mode of the correlation function g(φ) in eq. (4),high-dimensional oscillator dynamics described by eq. (1) k-clustered distribution is emphasized in the PDF U (φ).
10009-p2

Design and control of noise-induced synchronization patterns
Design of synchronization patterns. – Equation (3) indicates that we can design the stationary PDF
of the phase diﬀerence U (φ), i.e., the synchronization pattern, by varying the correlation function g(φ). Therefore,
given the power spectra Pξ (Ω), Pη (Ω), Pζ (Ω), and thesensitivity function Z(θ), we can try to find an optimalfilter f (τ ) that gives g(φ). In this study, rather than explicitly specifying the precise PDF of the oscillators as
the target, we aim to maximize its statistical property,
e.g., the degree of synchronization or clustering. This isbecause such macroscopic properties, rather than precise
functional forms of U (φ), are relevant in practical applications. Note also that we cannot generate arbitrary PDFs
but only optimize the PDF given in the form of eq. (3).
To design the optimal filter f (τ ), we introduce the following objective functional R{f } of f (τ ) characterizing
the statistical property of U (φ) as a measure for choosing
an optimal synchronization pattern:
R{f } =# +π

U (φ; f )q(φ)dφ,

(5)

−π

where we explicitly show the dependence of U (φ; f ) on f .
The function q(φ) determines what statistical property wefocus on. We try to design synchronization patterns withdesired statistical properties by choosing appropriate q(φ).
In the numerical simulations given below, we will use
the following functions for designing the synchronization
patterns: q1 (φ) = cos φ, q2 (φ) = δ(φ), q3 (φ) = cos 3φ,
and q4 (φ) = cos 2φ, where δ(·) denotes the Dirac delta
function. When we use q1 (φ), the objective functional
R{f } corresponds to the order parameter introduced inref. [18], which characterizes the degree of noise-induced
synchronization. When we use q2 (φ), the objective functional R{f } corresponds to the maximum of the PDF U (φ)
at φ = 0, which also characterizes the degree of synchronization, but in a more strict way, i.e., it counts only theoscillator pairs with exactly zero phase diﬀerence. When
we use q3 (φ), the objective functional R{f } characterizes
three-clustered states, in which three synchronized subgroups of oscillators are formed. Similarly, q4 (φ) characterizes two-clustered states.
In the following,"∞ we represent q(φ) as a Fourier series,q(φ) = q̃0 + 2 ℓ=1 q̃ℓ cos ℓφ, where the coeﬃcient q̃ℓ ∈ R
represents the weight"of the ℓ-th Fourier mode. Expanding∞1+
U (φ) as U (φ)=ℓ=1 ũℓ cos ℓφ, eq. (5) can be written
2π"∞
as R{f } = ℓ=1 q̃ℓ ũℓ . By finding optimal ũℓ for given q̃ℓ ,we can obtain a PDF U (φ) and a filter f (τ ) that maximizes
the objective functional R{f }.
Optimization of the filter. – By maximizing theobjective functional R{f }, we seek for the optimal filterf (τ ). However, unconstrained maximization of R{f } often leads to divergent f (τ ). We also need to take into account that our present theory is not valid for strong noisy
inputs, because the phase reduction method requires the
input given to the oscillators to be suﬃciently weak [11].

Thus, we should introduce some constraint on the filter
function f (τ ).
In this study, we formulate the constrained optimization
problem of the objective functional R{f } as follows:
maximize

R{f },(6)subject toσ 2 := ⟨Ij (t)2 ⟩ = C,(7)fwhere the condition eq. (7) constrains the variance σ 2 ofthe filtered noise Ij (t) to be a constant C. Using the power
spectra Pξ (Ω) and Pη (Ω) and the amplitude response A(Ω)
of the filter f (τ ), the variance σ 2 can be written as
σ2 =

# +∞
−∞

|A(Ω)|2 [Pξ (Ω) + Pη (Ω)]dΩ.

(8)

By solving the optimization problem described by eqs. (6)and (7), we can, in principle, obtain the optimal filter f (τ )for maximizing the objective functional R{f }.
Actually, we should also take into account that the optimal solution of eqs. (6) and (7) may not be implemented
in practice. The optimal amplitude response A(Ω) obtained as above often has"delta peaks at Ω = ℓω (ℓ ∈ Z),∞
i.e., |A(Ω)|2 = a0 δ(Ω) + ℓ=1 aℓ [δ(Ω − ℓω) + δ(Ω + ℓω)](aℓ is some coeﬃcient and ω is the natural frequencyof the oscillator), because the PDF U (φ) depends only
on the harmonic components Pξ (0), Pξ (ω), Pξ (2ω), . . . andPη (0), Pη (ω), Pη (2ω), . . . of the noise (see eqs. (3) and (4)).
Such a delta-peaked amplitude response A(Ω) corresponds
to a physically unrealistic filter that extracts only purely
harmonic components from the noise, which leads to phaselocking rather than noise-induced synchronization of theoscillators. Besides, such singular A(Ω) cannot be realized
in practical implementation of the linear filter f (τ ).
To overcome this problem, we restrict the class of A(Ω)
and further assume that the square of the amplitude response is expressed as a finite sum of narrow-band basis
functions as|Ã(Ω)|2 :=

m
$

ℓ=−m

c|ℓ| W (Ω − ℓω),(9)where cℓ (ℓ = 0, 1, 2, . . .) is a weight coeﬃcient, W (Ω)
represents a narrow-band basis function prespecified before the optimization process, e.g., a Gaussian function,and m is the maximum wave number of the filter. We assume that the basis function W (Ω) is localized in the range|Ω| < ω, i.e., W (Ω) ≈ 0 holds for |Ω| ≥ ω. The parameterm should be suﬃciently large to obtain a good filter. The
restricted amplitude response Ã(Ω) is experimentally feasible, because W (Ω−ℓω) in eq. (9) can be implemented bya band-pass filter that passes frequencies around Ω = ℓω.
We introduce%a new parameter%
β = [β0 , β1 , . . . , βm ]⊤ ∈m+1
as β0 = |b0 c0 | and βℓ = 2|bℓ cℓ | for ℓ = 1, . . . , m,R
! +∞where bℓ (ℓ = 0, 1, 2, . . .) is defined as bℓ = −∞ W (Ω −ℓω)[Pξ (Ω) + Pη (Ω)]dΩ. The optimization problem (6)10009-p3

W. Kurebayashi et al.
and (7) with the above restriction can then be expressed as# +πR̃(β) :=Ũ (φ; β)q(φ)dφ,(10)
maximizeβsubject to−π

σ 2 ≈ ||β||2 =

m−1
$
l=0

|βℓ |2 = C,

(11)

where Ũ (φ; β) is the PDF of the phase diﬀerence φ obtained by plugging eq. (9) into eq. (4), and eq. (11) follows
from eqs. (8) and (9) and the definition of β. Thus, we
can employ β as a design parameter of the optimization
problem described by eqs. (10) and (11).
Optimization algorithm. – To solve the optimization
problem given by eqs. (10) and (11), we use the gradientdescent algorithm. We randomly choose an initial value
β (0) and iteratively calculate β (j) for j ≥ 1, where β (j) isthe design parameter β at j-th iteration. At each iteration, we update β (j) asβ̃ (j+1) = β (j) + α∇β R̃(β (j) )(12)√and normalize it as β (j+1) = C β̃ (j+1) /||β̃ (j+1) ||, so that
β (j) satisfies the constraint (11). Here, α is a constant that
controls the step size, and ∇β R̃(β) ∈ Rm+1 represents
the gradient of R̃(β) with respect to β, i.e., ∇β R̃(β) =R̃(β)R̃(β) ⊤[ ∂ ∂β, . . . , ∂∂β] .
0m
For simplicity of notation, we define u(φ) := 1/[g(0) −g(φ) + h(0)]. The normalization constant is given by ū =
! +π
(j)
)
u(φ)dφ. Then, the gradient ∂ R̃(βin eq. (12) can∂βℓ−π! +π 1 ∂u(φ)∂ R̃(β)∂ ū
be expressed as ∂βℓ = −π ū2 [ ∂βℓ ū − u(φ) ∂β]q(φ)dφ,ℓ
and the gradients (∂u(φ)/∂βℓ ) and (∂ ū/∂βℓ ) can be calculated from eqs. (3), (4) and (9) and the definition of β as
"m−1
∂u(φ)
βℓ |zk |2
2[Pξ (kω)(cos kφ −k=0 W (kω − ℓω) bℓ∂βℓ = 2u (φ)!
+π ∂u(φ)∂ ū1) − Pη (kω)], ∂βℓ = −π ∂βℓ dφ.
Because W (Ω) is localized in the range |Ω| < ω,
the above expression can be simplified as follows:
2
∂u(φ)
ℓ|
= 2u2 (φ)W (0) βℓ |z[Pξ (ℓω)(cos ℓφ − 1) − Pη (ℓω)],∂βℓbℓwhich reduces the computational cost of the optimization
process. Using the optimized β, we can obtain the optimal amplitude response Ã(Ω) as the square root of eq. (9),2whose coeﬃcients c0 , . . . , cm are given by c0 = |βb00| , cℓ =|βℓ |22bℓ ,for ℓ = 1, 2, . . . , m. Thus, using, e.g., the leastsquares method [29], we can calculate the optimal filterf (τ ) from the amplitude response Ã(Ω).
Numerical simulations. – To confirm the validityof our method, we performed numerical simulations using several examples of the objective functionals. In
the first example, we use the FitzHugh-Nagumo (FHN)
model [30,31] of a periodically firing neuron. This model
has a two-dimensional state variable (u, v), which obeysv̇(t) = v − v 3 /3 − u + I0 + I(t) and u̇(t) = µ(v + c − du),µ = 0.08, c = 0.7, d = 0.8, and I0 = 0.875. The frequencyof the oscillation is approximately ω = 0.173 and the noisy

Fig. 1: (Color online) FitzHugh-Nagumo model. (a) Sensitivity function Z(θ), (b) Fourier coeﬃcients |zℓ |, and (c) power
spectra of the noisy inputs Pξ (Ω), Pη (Ω) and Pζ (Ω).

input Ij (t) is given to v(t). The sensitivity function Z(θ)
to I(t) and its Fourier coeﬃcients are shown in figs. 1(a)and (b).
As the noisy inputs, we use the Ornstein-Uhlenbeck
noise, whose power spectra are given by Pξ (Ω) =P̃OU (Ω; 0.5), Pη (Ω) = Pζ (Ω) = 0.1 · P̃OU (Ω; 0.5), withP̃OU (Ω; γ) := γ 2 /(γ 2 + Ω2 ) (fig. 1(c)). As the basis function W (Ω), we employ a rectangular function, W (Ω) =1 (|Ω| < ω/2), 0 (otherwise). The other parameters are set
as follows: the coupling strength to the noise is ϵ = 0.01,
the variance of the filtered noise is C = 10, and the maximum wave number of the filter is m = 5 (because |zℓ | issmall when ℓ ≥ 6). The parameter α used for the gradientdescent is α = 0.5. Note that the gradient descent algorithm finds only a local optimum and does not guarantee
global optimality. In order to obtain the global optimum,
the algorithm should be repeated from suﬃciently many
initial states β (0) .
Figure 2 shows the numerical results for the FHN model.
For the functions q1 (φ) and q2 (φ) defined previously,synchronized states are successfully formed (figs. 2(b)and (d)). For the sinusoidal q1 (φ), we obtain a filter
that emphasizes only the first Fourier mode (fig. 2(a))(see eq. (4)), which results in a bell-shaped PDF U (φ)(fig. 2(b)). For the delta-shaped q2 (φ), in contrast, we
obtain a nontrivial filter that consists of multiple modes(fig. 2(c)). In this case, the PDF U (φ) has a sharper peakthan that for q1 (φ) (fig. 2(d)) and a more precisely synchronized state is realized. Note that the high-frequency
components of Ã(Ω) in fig. 2(c), which are stronger thanthe low-frequency components, do not significantly affect the statistical property of Ij (t), because P (ℓω) is
suﬃciently small for large ℓ. Therefore, we can safely
neglect the high-frequency components of Ã(Ω) whose

10009-p4

Design and control of noise-induced synchronization patterns

Fig. 2: (Color online) FitzHugh-Nagumo model. (a), (c) and(e): amplitude response A(Ω) of the optimal filter designed by
the proposed method and (b), (d) and (f): probability density
function U (φ) of the phase diﬀerence φ for q1 (φ) ((a) and (b)),
q2 (φ) ((c) and (d)), and q3 (φ) ((e) and (f)). The insets display
snapshots of the oscillators in the v-u plane.

Fig. 4: (Color online) Hodgkin-Huxley model. (a), (c) and(e): amplitude response A(Ω) of the optimal filter designed by
the proposed method and (b), (d) and (f): probability density
function U (φ) of the phase diﬀerence φ for q1 (φ) ((a) and (b)),
q2 (φ) ((c) and (d)) and q4 (φ) ((e) and (f)). The insets display
snapshots of the oscillators in the V -m plane, where m is achannel variable [32].

Fig. 3: (Color online) Hodgkin-Huxley model. (a) Sensitivity
function Z(θ) and (b) Fourier coeﬃcients.

In the second experiment, we use the Hodgkin-Huxley(HH) model [32]. It also models periodic firing of a neuron, but it has more realistic, higher-dimensional dynamics without odd symmetry, in contrast to the FHN model.
We apply the noisy input Ij (t) as well as a constant input I0 = 10 to the V variable (i.e., membrane potential) of
the HH model. The oscillation frequency is approximately
ω = 0.438, and the sensitivity function and its Fourier coeﬃcients are shown in figs. 3(a) and (b). In addition to
q1 (φ) and q2 (φ), we use q4 (φ) for optimization with theaim of forming two-clustered states. The power spectraPξ (Ω), Pη (Ω) and Pζ (Ω), the parameters C and α, andthe basis function W (Ω) are the same as before. The noise
intensity is ϵ = 0.1 and the maximum wave number of thefilter is m = 4 (because |zℓ | almost vanishes at ℓ = 5).
Figure 4 shows the numerical results for the HH model.
Synchronized states are successfully formed for q1 (φ) andq2 (φ) as shown in figs. 4(b) and (d). When we use thedelta-shaped q2 (φ), we obtain a nontrivial filter consistingof multiple modes (fig. 4(c)) and the PDF U (φ) has awave numbers are larger than m. For the function q3 (φ),
we obtain a filter that emphasizes only the third Fouriermode (fig. 2(e)), which yields a three-clustered state as
expected (fig. 2(f)). Note that we cannot form a twoclustered state in the FHN model, because the phase response property of this model has odd symmetry, i.e., thesecond Fourier coeﬃcient |z2 | of Z(θ) is vanishingly small
as shown in fig. 1(b).

10009-p5

W. Kurebayashi et al.
sharper peak than the case with the sinusoidal q1 (φ)(fig. 4(d)). In contrast to the FHN model, we can realize
a two-clustered state as shown in fig. 4(f), because Z(θ) of
the HH model has a suﬃciently large second Fourier component as shown in fig. 3(b). Note that the realizabilityof a particular state is determined by the sensitivity function Z(θ) that characterizes the phase response propertyof the driven oscillator, rather than by the dimensionality
or complexity of the oscillator model.
Summary and discussion. – We have proposeda method for designing and controlling various noiseinduced synchronization patterns by filtering the input
noise, including the synchronized and clustered states. By
numerical simulations, the validity of the method has been
confirmed for two types of limit-cycle oscillators. These results will provide a theoretical basis for optimizing noiseinduced synchronization by filtering the input noise.
Though some previous works [33,34] proposed optimization methods for the phase response property of the oscillator to enhance noise-induced synchronization, thoseworks considered only the Lyapunov exponent of the phase(i.e., the exponential decay rate of the small phase diﬀerence between two oscillators), so that they could not fully
characterize the synchronized states and could not be used
to design various synchronization patterns as describedin this letter. More importantly, in contrast to previous
works [33,34] that gave the optimal phase response property of the oscillator, our present study provides a method
to generate optimal noisy inputs to the oscillator, which
can be implemented much more easily than designing theoscillator response. Thus, our method can be useful invarious real-world applications, e.g., energy-eﬃcient synchronization control in wireless sensor networks [27].
Finally, though we have not considered the eﬀect
of diﬀerences in the natural frequency of the oscillators [17,22,24] in this letter, it is often significant in practical applications. Extension of the present method tononidentical oscillators will be an important future work.
∗∗∗
The authors are grateful to H. Yasuda,
M. Harashima, Y. Honda, Y. Horio, and K. Aihara
for fruitful discussions. Financial support by KAKENHI(25540108, 26103510, 26120513), CREST Kokubu projectof JST, and FIRST Aihara project of JSPS are gratefully
acknowledged.
REFERENCES[1] Mainen Z. F. and Sejnowski T. J., Science, 268 (1995)
1503.
[2] Neiman A. B. and Russell D. F., Phys. Rev. Lett., 88(2002) 138103.
[3] Galán R. F., Fourcaud-Trocmé N., Ermentrout
G. B. and Urban N. N., J. Neurosci., 26 (2006) 3646.

[4] Ermentrout G. B., Galán R. F. and Urban N. N.,
Trends Neurosci., 31 (2008) 428.
[5] Yoshida K., Sato K. and Sugamata A., J. Sound Vib.,
290 (2006) 34.
[6] Utagawa A., Asai T., Hirose T. and Amemiya Y.,
IEICE Trans. Fundam. Electron. Commun. Comput. Sci.,
91 (2008) 2475.
[7] Zhou T., Chen L. and Aihara K., Phys. Rev. Lett., 95(2005) 178103.
[8] Uchida A., McAllister R. and Roy R., Phys. Rev.
Lett., 93 (2004) 244102.
[9] Toral R., Mirasso C. R., Hernández-Garcı́a E. andPiro O., Chaos, 11 (2001) 665.
[10] Zhou C. and Kurths J., Phys. Rev. Lett., 88 (2002)
230602.
[11] Kuramoto Y., Chemical Oscillatoins, Waves andTurbulence (Springer, Berlin) 1984.
[12] Teramae J.-N. and Tanaka D., Phys. Rev. Lett., 93(2004) 204103.
[13] Goldobin D. S. and Pikovsky A., Physica A, 351(2005) 126.
[14] Goldobin D. S. and Pikovsky A., Phys. Rev. E, 71(2005) 26220.
[15] Nakao H., Arai K., Nagai K., Tsubo Y. andKuramoto Y., Phys. Rev. E, 72 (2005) 26220.
[16] Nakao H., Arai K. and Kawamura Y., Phys. Rev.
Lett., 98 (2007) 184101.
[17] Yoshimura K., Davis P. and Uchida A., Prog. Theor.
Phys., 120 (2008) 621.
[18] Marella S. and Ermentrout G. B., Phys. Rev. E, 77(2008) 041918.
[19] Hata S., Shimokawa T., Arai K. and Nakao H., Phys.
Rev. E, 82 (2010) 036206.
[20] Goldobin D. S., Teramae J.-N., Nakao H. and
Ermentrout G. B., Phys. Rev. Lett., 105 (2010) 154101.
[21] Kurebayashi W., Fujiwara K. and Ikeguchi T., EPL,
97 (2012) 50009.
[22] Burton S. D., Ermentrout G. B. and Urban N. N.,
J. Neurophysiol., 108 (2012) 2115.
[23] Lai Y. M. and Porter M. A., Phys. Rev. E, 88 (2013)
012905.
[24] Zhou P., Burton S. D., Urban N. N. and
Ermentrout G. B., to be published in Front. Comput.
Neurosci., 113 (2013), doi: 10.3389/fncom.2013.00113.
[25] Satake A. and Iwasa Y., J. Ecol., 90 (2002) 830.
[26] Lyles D., Rosenstock T. S., Hastings A. and Brown
P. H., J. Theor. Biol., 259 (2009) 701.
[27] Yasuda H. and Hasegawa M., IEICE Trans. Commun.,96 (2013) 2749.
[28] Yoshimura Y. and Arai K., Phys. Rev. Lett., 101 (2008)
154101.
[29] Jackson L. B., Digital Filters and Signal Processing(Kluwer Academic Publishers, Norwell) 1996.
[30] FitzHugh R., Biophys. J., 1 (1961) 445.
[31] Nagumo J., Arimoto S. and Yoshizawa S., Proc. IRE,
50 (1962) 2061.
[32] Tuckwell H. C.,
Introduction to Theoretical
Neurobiology (Cambridge University Press) 1988.
[33] Abouzeid A. and Ermentrout G. B., Phys. Rev. E, 80(2009) 011911.
[34] Hata S., Arai K., Galán R. F. and Nakao H., Phys.
Rev. E, 84 (2011) 016229.

10009-p6