Colored noise induces synchronization of limit cycle oscillators

Colored noise induces synchronization of limit cycle oscillators

Wataru Kurebayashi\inst1    Kantaro Fujiwara\inst1    Tohru Ikeguchi\inst1,2
Abstract

Driven by various kinds of noise, ensembles of limit cycle oscillators can synchronize. In this Letter, we propose a general formulation of synchronization of the oscillator ensembles driven by common colored noise with an arbitrary power spectrum. To explore statistical properties of such colored noise-induced synchronization, we derive the stationary distribution of the phase difference between two oscillators in the ensemble. This analytical result theoretically predicts various synchronized and clustered states induced by colored noise and also clarifies that these phenomena have a different synchronization mechanism from the case of white noise.

Introduction

Driven by common noise, many nonlinear dynamical systems can synchronize. This phenomenon is called noise-induced synchronization, which is observed in various kinds of the nonlinear dynamical systems, for example, neural networks [1, 2], electric circuits [3], electronic devices [4], microbial cells [5], lasers [6] and chaotic dynamical systems [7, 8]. It has been theoretically proven that limit cycle oscillators can synchronize driven by common noise [9]. Many studies have investigated the synchronization property in case of various types of drive noises, for example, Gaussian white noise [10, 11, 12] and Poisson impulses [13]. In ref. [10], using a formulation of limit cycle oscillators driven by common and independent Gaussian white noises, Nakao et al. analytically obtained the probability density function (PDF) of phase differences between two oscillators, which enables us to effectively characterize the synchronization property. However, although there are some numerical studies[8, 14, 15], analytical conventional studies are limited to the case that drive signals are white noise (temporally uncorrelated noise). If we can assume that the drive signal is white noise, we can use the Fokker-Planck approximation [16] to explore statistical properties of oscillator ensembles. However, such an ideal condition is rare in the real world. For example, in neural circuits, it is known that colored noise with negative autocorrelation plays a key role to propagate synchronous activities [17]. However, it still remains unclear how the oscillators behave if they are driven by common colored noise.

Recently, it has been clarified how a limit cycle oscillator behaves if it is driven by colored non-Gaussian noise [18, 19, 20]. In this Letter, utilizing effective white-noise Langevin description proposed in ref. [19], we extend the formulation in ref. [10] to colored noise that has an arbitrary power spectrum. We then analytically derive the PDF of the phase difference between the oscillators if these oscillators are driven by common colored noise. We also conducted numerical simulations to verify our analytical results. The results show that the PDF of the phase difference explicitly depends on the power spectrum of the drive noise.

Model

We used the following system that consists of identical limit cycle oscillators subject to common and independent multiplicative colored noises. The dynamics of the th oscillator is described by

 ˙X(j) = F(X(j))+√DG(X(j))ξ(t) (0) +√ϵH(X(j))η(j)(t),

for , where is the -dimensional state variable of the th oscillator; is an unperturbed vector field that has a stable -periodic limit cycle orbit ; is the common noise, which drives all of the oscillators; () is the independent noise, which is received independently by each oscillator; and represent how the oscillators are coupled to the common and independent noises; and are parameters to control the intensities of the common and independent noises. We introduced the following three assumptions: (i) and are independent, identically distributed zero-mean colored noises, namely, , , , and (), where denotes the transpose and represents the temporal average; (ii) and can be approximated as the convolution of an arbitrary filter function and white noise; and (iii) and have correlation times shorter than the time scale of the phase diffusion ().

To characterize the statistical properties of the drive noises and , we define correlation matrices and as and (). For the sake of simplicity, we assumed that all independent noises have the same statistical property characterized by . The (, )th element of is the cross correlation function of the th and th elements of the common noise . The diagonal elements of are autocorrelation functions. In the same way, we can characterize the statistical property of by using .

Phase reduction

Under the assumption that the noise intensity is sufficiently weak ( and ), we can apply the phase reduction method [21, 20] to eq. (Model). By introducing a phase variable , eq. (Model) is reduced to the following phase equation:

 ˙ϕ(j) = ω+√DZG(ϕ(j))⋅ξ(t) (0) +√ϵZH(ϕ(j))⋅η(j)(t)+O(D,ϵ),

where is a phase variable that corresponds to the state of the th oscillator , () is the natural frequency, and and are the phase sensitivity functions that represent the linear response of the phase variable to the drive noises[21, 20]. The phase sensitivity functions and are defined as and . As discussed in Ref. [20], the term is necessary to describe the exact phase dynamics, while the phase diffusion is not affected by the term. As we will focus on the phase diffusion in the following sections, we do not take this term into account.

Effective Langevin description

To quantify the synchronization property without loss of generality, we consider the relationship of only two oscillators, that is, the two-body problem of and , and define the phase difference (). As we focus on the stochastic dynamics of , we define as the PDF of the phase difference . Utilizing the effective white-noise Langevin description [19], the evolution of is described by the following effective Fokker-Planck equation:

 ∂f∂t+∂∂θv(1)(θ)f−12⋅∂2∂θ2v(2)(θ)f=0, (0)

where and are effective drift and diffusion coefficients. We have the drift coefficient because . Meanwhile, the diffusion coefficient is obtained as

 v(2)(θ) = ∫+∞−∞dτ⟨[˙θ(t)−⟨˙θ⟩][˙θ(t−τ)−⟨˙θ⟩]⟩ (0) = ∫+∞−∞dτ⟨[˙ϕ(1)(t)−˙ϕ(2)(t)] [˙ϕ(1)(t−τ)−˙ϕ(2)(t−τ)]⟩

where represents the temporal average. For simplicity of notation, we define as . Then, we obtain

 v(2)(θ)=d11+d22−d12−d21=2d11−2d12. (0)

The phase variable can be expanded as by using and as expansion parameters, where , and () are approximate perturbed solutions of . We have , and . Using these perturbed solutions, eq. (Phase reduction) can be written as . Using this approximation and the fact that and , we obtain

 djk = D∫+∞−∞dτ⟨˙ϕ(j)D,1(t)˙ϕ(k)D,1(t−τ)⟩ (0) +ϵ∫+∞−∞dτ⟨˙ϕ(j)ϵ,1(t)˙ϕ(k)ϵ,1(t−τ)⟩ +O(D32,ϵ32).

Thus, using eq. (Effective Langevin description), we can calculate as follows:

 d11 = D2π∫+∞−∞dτ∫+π−πdϕ (0) ZG(ϕ)⊤Cξ(τ)ZG(ϕ−ωτ) +ϵ2π∫+∞−∞dτ∫+π−πdϕ ZH(ϕ)⊤Cη(τ)ZH(ϕ−ωτ)+O(D32,ϵ32).

In the same way, is given by

 d12 = D2π∫+∞−∞dτ∫+π−πdϕ (0) ZG(ϕ)⊤Cξ(τ)ZG(ϕ−θ−ωτ) +O(D32,ϵ32).

The detailed derivations of eqs. (Effective Langevin description) and (Effective Langevin description) are shown in Appendix A.

Finally, from eqs. (Effective Langevin description), (Effective Langevin description) and (Effective Langevin description), we have the efficient diffusion coefficient :

 v(2)(θ)=2D[g(0)−g(θ)]+2ϵh(0), (0)

where and are correlation functions defined as

 g(θ) = 12π∫+∞−∞dτ∫+π−πdϕ (0) ZG(ϕ)⊤Cξ(τ)ZG(ϕ−θ−ωτ), h(θ) = 12π∫+∞−∞dτ∫+π−πdϕ (0) ZH(ϕ)⊤Cη(τ)ZH(ϕ−θ−ωτ).

If we assume that the drive noise is white, namely, , eqs. (Effective Langevin description) and (Effective Langevin description) are exactly equivalent to eq. (6) in Ref. [10], where is an identity matrix. The results show that eqs. (Effective Langevin description) and (Effective Langevin description) are a natural generalization of eq. (6) in Ref. [10].

We obtain the explicit form of the Fokker-Planck equation of eq. (Effective Langevin description) from eqs. (Effective Langevin description)–(Effective Langevin description). The stationary distribution of the phase difference is given as the stationary solution of eq. (Effective Langevin description). Then, if we put in eq. (Effective Langevin description), we obtain

 f0(θ)=νv(2)(θ)=ν′D[g(0)−g(θ)]+ϵh(0), (0)

where and () are normalization constants.

Fourier representation

To understand the results obtained in the previous section, we rewrite the correlation functions defined in eqs. (Effective Langevin description) and (Effective Langevin description) by using the Fourier representation. We introduced the Fourier series expansion of the phase sensitivity functions and as and , where denotes the imaginary unit and () and () are Fourier coefficients ().

Subsequently, we define and as the Fourier transforms of and , that is, and . Let us note that and are Hermitian matrices, namely, and because and from their definitions, where denotes the adjoint. The th elements of and represent the cross spectra of the th and th elements of and . In particular, the diagonal elements of and represent the power spectra.

Using the Fourier representations defined above, we can obtain the Fourier representations of the correlation functions and :

 g(θ)=+∞∑l=−∞gleilθ, h(θ)=+∞∑l=−∞hleilθ, (0)

where () and () are Fourier coefficients (). The derivations of and will be shown in Appendix B.

These expressions clearly suggest that the correlation functions and only depend on and (), that is, the other frequency components can be neglected. In the next section, we will demonstrate that colored noise induces various synchronized and clustered states, which are clearly explained by eq. (Fourier representation).

Numerical simulations

To demonstrate the validity of our results, we perform numerical experiments for two types of limit cycle oscillators. The first example is the Stuart-Landau oscillator, which takes the normal form of the supercritical Hopf bifurcation [21]: , , where is a state variable and and are parameters. In the simulation, we fixed , , , and , where denotes an diagonal matrix that has the diagonal elements . This model is reduced to the phase equation that has the natural frequency and the phase sensitivity function .

In the simulation, we use a two-dimensional drive noise that has the correlation matrix defined as and , where and are parameters that represent the peak frequency and the characteristic decay time. We define , the Fourier transform of , as . A drive noise characterized by can be generated by the damped noisy harmonic oscillator (See eqs. (43)–(49) in Ref. [19] for details).

We use the common noises with , and and the independent noise with . The power spectra of these common noises are shown in fig. 1 (a). From eq. (Fourier representation), the correlation functions and are given by and , for 1, 3 and 5, which correspond to the three types of the common noise. The derivations of and will be shown in Appendix C.

The correlation function calculated above indicate that the effective intensity of the common noise depends on the peak frequency and is maximal at . It means that the synchronous degree is maximized at . In fig. 1 (b)–(d), we compared the results of the direct numerical simulation using the Stuart-Landau oscillator and its corresponding phase oscillator with the analytical results. All PDFs are well fitted by the theoretical curves. Our theory clearly predicts that the highest synchronous degree is realized at .

The second example is the FitzHugh-Nagumo oscillator [22, 23]: , , where is a state variable and , , and are parameters. In the simulation, we fixed , , , , , and . For these parameters, this oscillator has the natural frequency . This oscillator models bursting behavior of a neuron, and only the first variable , which corresponds to the membrane potential of a neuron, is subject to noise.

In the simulation, we use the one-dimensional noise that has the correlation function . Different from the first example, we use the same parameters for both the common and independent noises. We used two parameter sets and . The power spectra of these drive noises are shown in fig. 2 (a). We obtain the correlation function () numerically as shown in fig. 2 (b).

In fig. 2 (c) and (d), we compared the results of the direct numerical simulation with the analytical results. The numerical results are in good agreement with the theoretical results. As theoretically predicted, a 3-cluster state is realized as shown in fig. 2 (d). If oscillators are driven by white noise, clustered states are induced only by multiplicative noise [10]. However, in case of colored noise, clustered states are induced not only by multiplicative noise but also by additive noise.

In the third example, we used the Hodgkin-Huxley oscillator[25], which enables us to demonstrate whether the theory is applicable to higher-dimensional limit cycle systems. We use green noise used in ref. [8], which is generated by applying a high-pass filter to white noise. The power spectrum is shown in fig. 3 (a). Different from the periodic noise characterized by , the green noise has a vanishing spectrum as . In the simulation, for the sake of simplicity, we used the same type of drive noise for the common and independent noises, and we set the noise intensities . Fig. 3 (b)–(d) compare the theoretical and numerical results, which show that our theory is also valid for these cases.

Summary and discussions

In this Letter, we extended a formulation to analyze various synchronized and clustered states of uncoupled limit cycle oscillators driven by common and independent colored noises. Using this formulation, we derived the probability density function of the phase difference and rewrote it by the Fourier representation. The obtained expressions clearly show that the synchronization property depends on the power spectrum of the drive noises. Such dependency has already been reported experimentally. For example, in ref. [24], the reliability, or synchronization across trials, is explored in neuronal responses to periodic drive inputs with various frequencies. The reliability is maximized at a certain frequency, which is similar to our results shown in fig. 1. Our results in this Letter supports the results in ref. [24] theoretically, because a neuron in a oscillatory state can be regarded as a noisy limit cycle oscillator.

Generally, noise in the real world often has a non-flat and characteristic power spectrum. In this sense, our formulation is a useful tool to estimate the synchronization property for both theoretical and practical aspects. Namely, the results obtained in this Letter can be applied to a wide range of purposes from mathematical modelings to technological problems.

The authors would like to thank S. Ogawa and AGS Corp. for their encouragement on this research project.

References

• [1] Mainen Z. F. Sejnowski T. J., Science, 268 (1995) 1503.
• [2] Galán R. F., Fourcaud-Trocmé N., Ermentrout G. B. Urban N. N., J. Neurosci., 26 (2006) 3646.
• [3] Yoshida K., Sato K. Sugamata A., J. Sound Vib., 290 (2006) 34.
• [4] Utagawa A., Asai T., Hirose T. Amemiya Y., IEICE Trans. Fundam., 91 (2008) 2475.
• [5] Zhou T., Chen L. Aihara K., Phys. Rev. Lett., 95 (2005) 178103.
• [6] Uchida A., McAllister R. Roy R., Phys. Rev. Lett., 93 (2004) 244102.
• [7] Zhou C. Kurths J., Phys. Rev. Lett., 88 (2002) 230602.
• [8] Wang Y., Lai Y.-C., Zheng Z.,Phys. Rev. E, 79 (2009) 056210.
• [9] Teramae J.-N. Tanaka D., Phys. Rev. Lett., 93 (2004) 204103.
• [10] Nakao H., Arai K. Kawamura Y., Phys. Rev. Lett., 98 (2007) 184101.
• [11] Yoshimura K., Davis P. Uchida A., Prog. Theor. Phys., 120 (2008) 621.
• [12] Nagai K. H. Kori H., Phys. Rev. E, 81 (2010) 065202.
• [13] Nakao H., Arai K., Nagai K., Tsubo Y. Kuramoto Y., Physical Review E, 72 (2005) 26220.
• [14] Yoshimura K., Valiusaityte I. Davis P., Phys. Rev. E, 75 (2007) 026208.
• [15] Hata S., Shimokawa T., Arai K. Nakao H., Phys. Rev. E, 82 (2010) 036206.
• [16] Risken H., The Fokker-Planck equation: Methods of solution and applications (Springer Verlag) 1996.
• [17] Câteau H. Reyes A. D., Phys. Rev. Lett., 96 (2006) 058101.
• [18] Teramae J. Tanaka D., Prog. Theor. Phys., 161 (2006) 360.
• [19] Nakao H., Teramae J.-N., Goldobin D. S. Kuramoto Y., Chaos, 20 (2010) 3126.
• [20] Goldobin D. S., Teramae J.-N., Nakao H. Ermentrout G. B., Phys. Rev. Lett., 105 (2010) 154101.
• [21] Kuramoto Y., Chemical oscillations, waves, and turbulence (Dover Publications) 2003.
• [22] FitzHugh R., Biophys. J., 1 (1961) 445.
• [23] Nagumo J., Arimoto S. Yoshizawa S., Proc. IRE, 50 (1962) 2061.
• [24] Fellous J. M., Houweling A. R., Modi R. H., Rao R. P. N., Tiesinga P. H. E. Sejnowski T. J., J. Neurophysiol., 85 (2001) 1782.
• [25] Hodgkin A. Huxley A., J. Physiol., 117 (1952) 500.

Appendix A: Derivations of eqs. ( ‣ Effective Langevin description) and ( ‣ Effective Langevin description)

Substituting and into eq. (Effective Langevin description), we obtain

 d11 = D∫+∞−∞dτ⟨[ZG(ϕ(1)0(t))⊤ξ(t)] (A.2 ) [ZG(ϕ(1)0(t−τ))⊤ξ(t−τ)]⟩ +ϵ∫+∞−∞dτ⟨[ZH(ϕ(1)0(t))⊤η(1)(t)] [ZH(ϕ(1)0(t−τ))⊤η(1)(t−τ)]⟩ +O(D32,ϵ32) = D∫+∞−∞dτ⟨ZG(ϕ(1)0(t))⊤ξ(t) +ϵ∫+∞−∞dτ⟨ZH(ϕ(1)0(t))⊤η(1)(t) η(1)(t−τ)⊤ZH(ϕ(1)0(t−τ))⟩ +O(D32,ϵ32).

We rewrite , , and by using their elements and obtain

 d11 = (A.3 ) ξk(t)ξl(t−τ)ZH,l(ϕ(1)0(t−τ))⟩ +ϵ∫+∞−∞dτ m∑k=1m∑l=1⟨ZH,k(ϕ(1)0(t)) η(1)k(t)η(1)l(t−τ)ZH,l(ϕ(1)0(t−τ))⟩ +O(D32,ϵ32),

where and are the th elements of and , and and are the th elements of and .

We assume that the phase variable and the drive noises and are approximately independent. Under this assumption, the temporal average can be divided into two parts; () and (). Thus, we obtain

 d11 = (A.4 ) +ϵ∫+∞−∞dτ m∑k=1m∑l=1⟨ZH,k(ϕ(1)0(t)) +O(D32,ϵ32) = D2π∫+∞−∞dτ∫+π−πdϕ m∑k=1m∑l=1ZG,k(ϕ)ZG,l(ϕ−ωτ)Cξ,kl(τ) +ϵ2π∫+∞−∞dτ ∫+π−πdϕ m∑k=1m∑l=1ZH,k(ϕ)ZH,l(ϕ−ωτ)Cη,kl(τ) +O(D32,ϵ32),

where and are the th elements of and . Finally, we rewrite eq. (A.4) by using , , and and obtain

 d11 = D2π∫+∞−∞dτ∫+π−πdϕ (A.5 ) ZG(ϕ)⊤Cξ(τ)ZG(ϕ−ωτ) +ϵ2π∫+∞−∞dτ∫+π−πdϕ ZH(ϕ)⊤Cη(τ)ZH(ϕ−ωτ) +O(D32,ϵ32).

In the same way, one can calculate as follows. We use the fact that and eliminate the phase variable of the second oscillator by substituting into , and then, we obtain

 d12 = D∫+∞−∞dτ⟨[ZG(ϕ(1)0(t))⊤ξ(t)] (A.6 ) [ZG(ϕ(2)0(t−τ))⊤ξ(t−τ)]⟩ +O(D32,ϵ32) = D∫+∞−∞dτ⟨ZG(ϕ(1)0(t))⊤ξ(t) +O(D32,ϵ32) = D2π∫+∞−∞dτ∫+π−πdϕ ZG(ϕ)⊤Cξ(τ)ZG(ϕ−θ−ωτ) +O(D32,ϵ32).

Appendix B: Derivation of eq. ( ‣ Fourier representation)

From eq. (Effective Langevin description), one can calculate the Fourier coefficient as follows. We introduce a new variable () and use the fact that is a Hermitian matrix. Then, we obtain

 gl = 12π∫+π−πdθg(θ)e−ilθ (B.1 ) = 12π∫+π−πdθ12π∫+∞−∞dτ∫+π−πdϕZG(ϕ)⊤ Cξ(τ)ZG(ϕ−θ−ωτ)e−ilθ = (12π∫+π−πdϕZG(ϕ)⊤e−ilϕ) (∫+∞−∞dτCξ(τ)eilωτ)(12π∫+π−πdχZG(χ)eilχ) = Y⊤G,l¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Pξ(lω)¯¯¯¯¯¯¯¯¯YG,l=Y†G,lPξ(lω)†YG,l = Y†G,lPξ(lω)YG,l,

where denotes the complex conjugate. From eq. (Effective Langevin description), can be derived likewise.

Appendix C: Derivations of the correlation functions g(θ) and h(θ)

For the Stuart-Landau oscillator we used in the simulations, we can calculate the Fourier coefficients and as and . Thus, from eq. (Fourier representation), the Fourier coefficient is given by and , where is a parameter. In the same way, the Fourier coefficient is given by and . Substituting and to eq. (Fourier representation), we can obtain the explicit forms of and .

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters