- Email: [email protected]
Noise-induced critical breakdown of phase lockings in a forced van der Pol oscillator
Physics Letters A 310 (2003) 407–414 www.elsevier.com/locate/pla
Noise-induced critical breakdown of phase lockings in a forced van der Pol oscillator Shinji Doi ∗ , Yohei Isotani, Ken-ichiro Sugimoto, Sadatoshi Kumagai Department of Electrical Engineering, Graduate School of Engineering, Osaka University, Osaka 565-0871, Japan Received 10 October 2002; accepted 25 February 2003 Communicated by C.R. Doering
Abstract The effect of noise on the period of the van der Pol oscillator forced by a sinusoidal signal is analyzed in detail. The dynamics of the forced nonlinear oscillator with noise is described by the evolution of probability densities by a Markov operator. The oscillator shows a 1 : 1 phase-locked oscillation whose period is equal to the period of the sinusoidal forcing when the forcing is sufficiently strong. Noise modulates the oscillator’s period as expected, but we show that this modulation is not uniform; the period of the oscillator changes drastically at a certain strength of noise. Thus this noise-induced breakdown of the phase locking can be considered as a stochastic bifurcation. The breakdown of stochastic phase-locking or the stochastic bifurcation is also analyzed from the viewpoint of the spectra of the Markov operator. 2003 Elsevier Science B.V. All rights reserved. PACS: 02.30.Oz; 02.50.Fz; 05.40.Ca Keywords: Nonlinear oscillator; Additive noise; Stochastic phase-locking; Stochastic kernel; Markov operator
1. Introduction The van der Pol oscillator is the representative model of a nonlinear oscillator and is widely utilized to represent various rhythmic phenomena in biology. Bonhoeffer–van der Pol or FitzHugh–Nagumo model [1] is the slight modification of the van der Pol oscillator and mimics qualitative behavior of a neuron. It is fundamentally important problem to examine how the period of the oscillator is modulated by the other oscillator or an external environment such as external noise [2]. * Corresponding author.
E-mail address: [email protected] (S. Doi).
Grasman and Roerdink [3] studied the van der Pol relaxation oscillator in the presence of noise and reduced the problem of examining the period of the oscillator to the analysis of the time necessary for a one-dimensional stochastic process to reach a boundary for the first time. Tateno et al. [4] extended this first-passage-time approach to the presence case of a periodic forcing and reduced the dynamics of the system to the evolution of probability densities by a Markov operator. Doi et al. [5] proposed the method to analyze phase-locking and bifurcation phenomena in the presence of noise using the spectra of the Markov operator. The present Letter advances this line of researches and investigates how the period of the sinusoidally
0375-9601/03/$ – see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(03)00386-4
408
S. Doi et al. / Physics Letters A 310 (2003) 407–414
forced van der Pol oscillator is changed by noise in detail. As the noise intensity increases, the oscillator phase-locked to the sinusoidal forcing is shown to change its period drastically from the period of the sinusoid at a certain strength of noise which can be considered as a breakdown point of a stochastic phase locking. We show that this breakdown point of a stochastic phase locking can be well characterized by the spectra of the Markov operator. The present Letter is organized as follows: Section 2 introduces the (piecewise-linearized) van der Pol oscillator forced by sinusoidal input in the presence of noise. Section 3 explains the mathematical framework to treat the noisy oscillator. Our method is based on the Markov operator which is a natural extension of a so-called Poincaré map or a return map in the noisefree (deterministic) case. Using this operator, the behavior of the system is described by the evolution of probability density functions of the phase of the sinusoidal input. Section 4 analyzes the mean period of the forced noisy oscillator in detail using the Markov operator and stochastic Arnold’s tongues (stochastic phaselocked regions) are presented in Section 5. The Letter concludes with some brief comments in Section 6.
2. van der Pol relaxation oscillator We consider a sinusoidally-forced (and piecewiselinearized) van der Pol relaxation oscillator in the presence of additive noise: dx = F (x, y), dt dy dW (t) = −x + v(t) + σ , dt dt 5 F (x, y) = y − x + |x + 1| − |x − 1| , 6 v(t) = A sin 2π(t/T + θ0 ) ,
(1)
Fig. 1. The x–y phase plane of (1), (2) in the limit of = 0 and the N-shaped x-nullcline (x˙ = 0). The y-nullcline (y˙ = 0) coincides with the y-axis. The closed orbit ABCD is a limit cycle. See text for details.
this case, the van der Pol oscillator shows a relaxation oscillation. Fig. 1 shows the x–y phase plane of Eqs. (1), (2). In the limit of = 0, the orbit with an initial value which is not on the N-shaped x-nullcline (x˙ = 0) instantaneously jumps to the x-nullcline in the horizontal direction since the horizontal velocity x˙ = F (x, y)/ becomes infinity unless F (x, y) = 0. Thus, all orbits are considered to move on the x-nullcline and both the sinusoid and noise terms of (2) modulate the velocity of the orbit along the x-nullcline. On the right (left) branch of the N-shaped x-nullcline, an orbit moves toward the point B (D, respectively). At the point B (D), the orbit instantaneously jumps to the point C (A, respectively), since we consider the limit of = 0. This limit is called singular since orbits are not differentiable at such jump points. The parallelogram ABCD is a limit cycle to which every orbit of (1), (2) approaches asymptotically.
(2)
where A, T and θ0 are the amplitude, the period and the initial phase of the sinusoidal input, respectively. W (t) is the standard Wiener process and σ dW (t)/dt denotes a Gaussian white noise with a noise intensity σ . Note that the cubic polynomial of the original van der Pol equation is piecewise linearized for the computational purpose [4]. Throughout, we consider a so-called singular (or relaxation) limit of = 0. In
3. Stochastic kernel and Markov operator Consider the noise-free (σ = 0) case and a Poincaré map or a return map as follows. Suppose that a state point starts at the point A of the x–y phase plane with an initial ‘phase’ θ0 of the sinusoidal input and that the point returns to the point A again after a time t. Then the phase of the sinusoidal input is changed from θ0 to θ1 : def
θ1 = p(θ0 ) = t/T + θ0 (mod 1).
S. Doi et al. / Physics Letters A 310 (2003) 407–414
Phase lockings of the noise-free case are usually analyzed by this return map p(θ ) and by the onedimensional mapping: θn+1 = p(θn ),
n = 0, 1, 2, . . . .
(3)
Note the sequence {θn } generated by this one-dimensional discrete-time dynamical system will be referred to as the orbit of (3) in the same way as the orbit of (1), (2) unless any confusions occur. Generally, in the case of m : n locking, the orbit or sequence {θn } asymptotically approaches to an n-periodic orbit (sequence): θ (0) , θ (1) , . . . , θ (n−1) , θ (0) , θ (1), . . . . We extend the deterministic map to the case with noise in the following way. In the noisy case, both variables θ0 and θ1 fluctuate owing to noise and thus are random variables Θ0 and Θ1 , respectively. Define a kernel function g(θ0 , θ1 ) using a conditional probability density function:
where the summation is taken for all ti such that θ1 = ti /T + θ0 (modulo 1). The first-passage-time pdf f0 (t, θ ) can be numerically calculated without simulations of the stochastic differential equations (1) and (2) [7]. Thus the kernel function can be obtained numerically also [4]. The function g takes relatively high values along the graph of the map p(θ0 ) of noisefree case; the values of g(θ0 , p(θ0 )) is much greater than other values of g(θ0 , θ1 ) if the noise intensity σ is not so large. Thus g can be considered as the stochastic extension of the (deterministic) return map. Using the kernel function g, we extend the system (3) to the noisy case. Let S denote a unit interval [0, 1] and D the set of absolutely-integrable non-negative functions with a unit L1 norm on S. A function which belongs to D is called a probability density function (pdf) or simply a density function. An (Markov [6]) operator P on D is defined by: Ph(θ ) =
g(θ0 , θ1 ) dθ1 = Pr{θ1 Θ1 θ1 + dθ1 | Θ0 = θ0 }.
f (t, θ ) =
f0 (s, θ )f0 t − s, θ ds,
0
θ = s/T + θ (mod 1),
(5)
since the phase of the sinusoid at the point C becomes θ after the elapse of time s on the branch AB of Fig. 1. The kernel function g(θ0 , θ1 ) of Eq. (4) is obtained by the function f (θ, t) as follows: fp (ti , θ0 ) × T , g(θ0 , θ1 ) = (6) θ1 =ti /T +θ0 (mod 1)
g(θ0 , θ )h(θ0 ) dθ0 ,
h ∈ D.
(7)
S
(4) Let the state point of the stochastic differential equations of (1), (2) start at the point A of Fig. 1 with the initial phase θ of the sinusoidal signal. We consider the time interval when the state point reaches the point B for the first time after leaving the point A. Let f0 (θ, t) be the probability density function (pdf) of this first-passage time from the point A to B. Owing to the symmetry of the trajectories of Fig. 1, the pdf of the first-passage time from the point C to D is same as f0 (θ, t). Thus the pdf f (θ, t) of the one-cycle time when the state point returns to the point A for the first time after leaving the point A with the initial phase θ of the sinusoid can be denoted by t
409
Let h0 (θ ) ∈ D denote the probability density function of the initial phase Θ0 when an orbit starts at the point A of the x–y phase plane. Then the density function of Θ1 when the orbit returns to the point A again is obtained by h1 (θ ) = Ph0 (θ ). Thus, the deterministic mapping (3) is extended to the system hn+1 (θ ) = Phn (θ ),
n = 0, 1, 2, . . . .
(8)
Note we investigate the asymptotic behavior of the sequence {hn (θ )} of probability density functions rather than the sequence {θn }. Let’s list several preliminary definitions [6]. A function h∗ (θ ) is called as the invariant density function of an operator P if the relation Ph∗ = h∗ holds. The invariant density is asymptotically stable if for any initial density function h0 ∈ D lim P n h0 − h∗ = 0.
n→∞
(9)
As is easily seen from its definition, the function g has the property g(θ0 , θ1 ) 0,
g(θ0 , θ1 ) dθ1 = 1 S
410
S. Doi et al. / Physics Letters A 310 (2003) 407–414
and is called a stochastic kernel. The inequality (10) inf g(θ0 , θ1 ) dθ1 > 0 θ0
S
is assumed to hold [4]. The operator with this property is known to have a unique asymptotically stable invariant density (see Corollary 5.7.1 of Ref. [6]). Thus the sequence {hn (θ )} produced by the operator P always approaches to a unique invariant density asymptotically as n → ∞. Note that the operator P has always a unique eigenvalue unity; the corresponding eigenfunction is the invariant density h∗ . The moduli of the other eigenvalues are less than unity. In particular, the modulus of the second (in the order of modulus) eigenvalue dominates the convergence speed of the initial pdf h0 to the invariant density h∗ in (9); the smaller the modulus of the second eigenvalue is the faster the convergence is.
4. Effects of noise on the period of the forced oscillator In this section, we investigate how the the forced oscillator (1), (2) varies when intensity σ is changed. The mean period oscillator is calculated by: T =
∞1
f (t, θ )h∗ (θ ) dθ dt,
period of the noise T of the
(11)
0 0
where the functions f and h∗ are the pdf of onecycle time and the invariant density (the stationary or asymptotic pdf of the phase) of the operator P, respectively. Fig. 2(a-i), (b-i), . . . , (f-i) show the mean periods numerically computed using (11) as a function of noise intensity σ . Panels (a-ii), (b-ii), . . ., (f-ii) show the moduli (absolute values) of the second and third eigenvalues of the operator P in the order of modulus (the first eigenvalue is always unity). Panels (a-i) and (a-ii) correspond to the case of A = 0: the absence case of the sinusoidal input. In this case, the mean period slightly decreases as the noise intensity σ is increased; noise accelerates the oscillator in the absence of the sinusoidal forcing. The inherent period of the van der Pol relaxation oscillator without both noise and sinusoidal input is
2(ln 7 − ln 3) ≈ 1.695 to which value the graph of the mean period approaches as the noise decreases. The moduli of second and third eigenvalues also decreases monotonically as noise increases. The plots of both second and third eigenvalues are completely overlapped; they are complex conjugate. Panels (b-i) and (b-ii) are the case that the amplitude and period of the sinusoid are A = 0.6 and T = 1.3, respectively. In this case, the sinusoidal input is not so strong that the oscillator could be phase-locked to the input when the noise is absent. The forced oscillator exhibits a so-called quasi-periodic oscillation without noise. The behavior of eigenvalues against noise is almost same as the previous case (a) while the mean period change is a little bit different from the previous case; the mean period does not change much as the noise intensity is increased. Although the oscillator is not phase-locked to the input, the input has much effect to the oscillator; the period of the oscillator is much shorter than the oscillator’s inherent period. As the noise intensity increases, the noise effect dominates over the effect of the sinusoidal input and the modulated period is finally growing up toward the inherent period (in the case of noise presence) shown in panel (a-i). Panel (c) corresponds to the case when the amplitude of the input is larger (A = 0.8). In this case, a 1 : 1 phase-locking appears when noise is absent. The mean period of (c-i) is nearly same as the period of the sinusoidal input in the range of small noise. Even if the noise intensity is increased, the mean period still stays around the input period. If the noise intensity is further increased, the mean period suddenly begins to increase. (Note that the inherent period of the van der Pol oscillator is about 1.7.) A striking fact is that the increase of the mean period by noise is not uniform; the mean period remains around the input period and then suddenly begins to increase. From this observation, we may say that the phase locking survives even when noise is added and the stochastic phaselocking breakdowns not gradually but critically at a certain noise intensity. The corresponding eigenvalues also present a critical behavior (see panel (c-ii)). The second and third eigenvalues are real numbers in the range of weak noise (this is a sign of stochastic phaselockings, see Ref. [5] for detail) and change their values from real to complex conjugate near the noise intensity σ = 0.065.
S. Doi et al. / Physics Letters A 310 (2003) 407–414
411
Fig. 2. Panel (i): mean periods of the forced van der Pol oscillator with noise are plotted as a function of the noise intensity σ . Panel (ii): moduli of 2nd and 3rd eigenvalues of the operator P are plotted as a function of noise intensity σ . (a) The absence case of the sinusoidal input: A = 0. (b) T = 1.3, A = 0.6. (c) T = 1.3, A = 0.8. (d) T = 1.3, A = 1.0. (e) T = 1.8, A = 0.2. (f) T = 1.8, A = 0.4.
Comparing the panels (c-i) and (c-ii), we can see that the noise intensity where the mean period starts increasing coincides nearly with the noise intensity where the eigenvalues change their values from real to complex. Thus we may say that the point of the eigenvalue change is the breakdown point of a stochastic phase-locking or the stochastic bifurcation point. Note that the second largest eigenvalue of the operator P dominates the convergence speed of the pdf sequence {hn (θ )} to the invariant density h∗ (θ ); the smaller the modulus of the eigenvalue is, the
faster the convergence becomes. Usually, if the noise intensity increases, the convergence to the stationary distribution becomes faster and the modulus of the second eigenvalue decreases toward zero. In the both panels (a-ii) and (b-ii), the moduli of the second and third eigenvalues are very close to unity and thus the convergence speed to the invariant density is very small in the range of small noise. This corresponds to the fact that the phase dynamics of the forced oscillator (in the case of noise-free or small noise) is similar to the simple rotation on the unit circle and the effect of phase spread-out to the invariant density is weak.
412
S. Doi et al. / Physics Letters A 310 (2003) 407–414
Differently from the both cases of (a-ii) and (b-ii), the modulus of the second eigenvalue in the small noise range is much less than unity in the case of (c-ii). This is because the nonlinearity of the dynamics in the case (c) is bigger than the case (a) and (b), thus the effect of phase spread-out to the invariant density caused by the nonlinearity is stronger than that of (a) and (b). Panels (d-i) and (d-ii) show the case of the stronger input. The mean period retains its value around the input period even if a very large noise is added. The transition point of both the mean period and the eigenvalue is shifted to a large noise direction. Also note that the moduli of the second eigenvalue are reduced totally than the case of (c) since the phase spread-out effect caused by the nonlinearity of the dynamics is increased by the increase of the amplitude of the sinusoidal input. Panels of (e) and (f) show the case of a different input period: T = 1.8 that is longer than the inherent
period of the van der Pol oscillator. Thus a large noise destroys the 1 : 1 phase-locking and decreases the mean period of the oscillator toward its inherent period (≈ 1.7). We can also see that the mean period changes its value drastically at a certain noise intensity and that the change point of eigenvalues is close to the change point of the mean period in both cases of (e) and (f).
5. Arnold tongue in the presence of noise In the above section, we have studied the noise effects on the forced nonlinear oscillator for several fixed values of the input amplitude and period. In this section we study the more global behavior on the input amplitude and period. Fig. 3(a) shows the Arnold tongues (boundaries of a phase-locked region) of noise-free (deterministic) forced van der Pol oscillator. Region of three phase-
Fig. 3. Arnold’s tongues. (a) Noise-free case (σ = 0). Regions of deterministic phase lockings of three types: 1 : 1, 2 : 1, and 5 : 3 lockings are shown. (b) Effect of noise on the 1 : 1 phase-locking. Noise-free and three cases of different noise intensities: σ = 0.02, 0.05, 0.09 are shown. (c) 2 : 1 phase-locking. Noise-free and the case of σ = 0.02 are compared. (d) 5 : 3 phase-locking. Noise-free and the case of σ = 0.02 are compared. See text for details.
S. Doi et al. / Physics Letters A 310 (2003) 407–414
locked patterns (1 : 1, 2 : 1, 5 : 3) are shown in the parameter space (plane) of the amplitude A and period T of the sinusoidal input. The 1 : 1 locking region occupies the most part of the parameter plane and the regions of 2 : 1 and 5 : 3 lockings are very narrow. Panel (b) shows both deterministic and stochastic Arnold’s tongues of the 1 : 1 phase-locking. The stochastic Arnold tongues are calculated as the points where the second and third eigenvalues change their values between real and complex values as mentioned in the above section. The widest region (the lowest plots) corresponds to the noise-free case and the stochastic phase-locking region shrinks slightly as the noise is increased. (The Arnold’s tongues of the σ = 0.02 case almost coincides with the noise-free case and are indistinguishable in this figure.) These calculations of the deterministic and stochastic Arnold’s tongues show the validity of our definition of the break-down point (boundaries) of stochastic phaselockings as the point of the eigenvalue change. Panel (c) is the case of 2 : 1 phase-locking: onecycle of the van der Pol oscillator is synchronized to the two cycles of the sinusoidal input. In this case, only noise-free and σ = 0.02 cases are shown since the stochastic 2 : 1 phase-locking is very fragile and the stochastic 2 : 1 locking region for stronger noise is very small or disappears. The method of computation of breakdown points (Arnold tongues) of 2 : 1 locking is same as the case of 1 : 1 locking. The last panel (d) shows both the deterministic and stochastic Arnold tongue of 5 : 3 lockings. By the same reason as the case (c), only noise-free and σ = 0.02 cases are computed. In this case, the points where the sixth and the seventh eigenvalues rather than the second and third ones change their values between real and complex are utilized to obtain the Arnold tongues, since the 5 : 3 phase-locking corresponds to the ‘threeperiodic’ behavior of the Markov operator while the 1 : 1 locking corresponds to the ‘one-periodic’ behavior of the operator (see Ref. [5] for details). From these figures we can see that the phaselocking regions (Arnold tongues) shrink by noise although this can be easily expected. We, however, stress that our numerical method which uses the Markov operator and its eigenvalues does not use numerical simulations of the stochastic differential equations (1), (2) and thus we can quantitatively study such a delicate behavior induced by noise in detail.
413
6. Discussion We have studied how the mean period of phaselocked oscillations in a sinusoidally-forced nonlinear oscillator is modulated by an additive noise using a piecewise linearized van der Pol oscillator. The dynamics of the noisy forced oscillator was expressed in terms of the evolution of probability density functions by a Markov operator [4]. This Markov operator is a natural extension of the deterministic (noise-free) Poincaré mapping or a return mapping which is the standard methodology to analyze forced oscillators. It should be note that the operator can be obtained only numerically but heavy numerical (Monte Carlo) simulations of the stochastic differential equations are not necessary. The operator is obtained by solving some integral equations numerically and thus can be computed with a high accuracy. A similar (Markov) approach has been taken by Plesser and Geisel [8] for the analysis of the noisy integrate-and-fire (IF) neuron model. Note that the IF neuron model is essentially same as the piecewise linearized van der Pol oscillator in the relaxation limit and the application of our method to the analysis of IF neuron models is a direct and simple application. Using the Markov operator, it is shown that the mean period of the noisy forced oscillator is modulated by noise not uniformly but critically; the mean period remains near the phase-locked period (the period of the sinusoidal input) against noise and drastically starts changing its value at a critical noise intensity. This behavior of the mean period modulation infers the existence of phase-lockings even in the case of noise presence (although the phase is not locked but varied by the presence of noise in a strict sense). In fact, near the critical change point of the mean period, the eigenvalues of the operator also showed the critical change: alternation of eigenvalue between real and complex values. The breakdown of deterministic phase-lockings is in itself the (saddle-node) bifurcation phenomenon. Although we have already proposed the method to analyze stochastic saddle-node and period-doubling bifurcations of the forced van der Pol oscillator [5] and of a one-dimensional mapping [9] and this method has already been utilized for the analysis of various noisy systems [10], there are still some place for the debate on the validity of our definition of stochastic
414
S. Doi et al. / Physics Letters A 310 (2003) 407–414
bifurcations. In this sense, the result of the present Letter that showed the abrupt change of the mean period at a critical noise level and the apparent coincidence of eigenvalue change, may give some basis on the usefulness and validity of our analysis method. A more rigorous and detailed discussion on the reason of the apparent closeness of both the critical change point of mean periods and the eigenvalues is necessary for the future research.
[3]
Acknowledgement
[7]
This research was partially supported by the Japanese Ministry of Education, Science, Sports and Culture, Grant-in-Aid (#10650429).
[8] [9]
References [1] R. FitzHugh, Biophys. J. 1 (1961) 445; S. Doi, S. Sato, Math. Biosci. 125 (1995) 229. [2] W. Horsthemke, R. Lefever, Noise-Induced Transitions, Theory and Applications in Physics, Chemistry, and Biology, Springer-Verlag, Berlin, 1984;
[4] [5] [6]
[10]
F. Moss, P.V.E. McClintock (Eds.), Noise in Nonlinear Dynamical Systems Vols. 1–3, Cambridge Univ. Press, Cambridge, 1987; T. Kapitaniak, Chaos in Systems with Noise, World Scientific, Singapore, 1988. J. Grasman, Asymptotic Methods for Relaxation Oscillations and Applications, Springer-Verlag, New York, 1987; J. Grasman, J.B.T.M. Roerdink, J. Stat. Phys. 54 (1989) 949. T. Tateno, S. Doi, S. Sato, L.M. Ricciardi, J. Stat. Phys. 78 (1995) 917. S. Doi, J. Inoue, S. Kumagai, J. Stat. Phys. 90 (1998) 1107. A. Lasota, M.C. Mackey, Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics, Springer, Berlin, 1994. A. Buonocore, A.G. Nobile, L.M. Ricciardi, Adv. in Appl. Probab. 19 (1987) 784. H.E. Plesser, T. Geisel, Phys. Rev. E 59 (1999) 7008. J. Inoue, S. Doi, S. Kumagai, in: Proceedings of NOLTA’97, International Symposium on Nonlinear Theory and its Applications, Hilton Hawaiian Village, Hawaii, 1997, pp. 13–16; J. Inoue, S. Doi, S. Kumagai, in: Proceedings of NOLTA’98, Le Regent, Crans-Montana, Switzerland, 1998, pp. 871–874; J. Inoue, S. Doi, S. Kumagai, Phys. Rev. E 64 (2001) 056219. T. Tateno, J. Stat. Phys. 92 (1998) 675; T. Yamanobe, K. Pakdaman, Biol. Cybern. 86 (2002) 155.
Noise-induced critical breakdown of phase lockings in a forced van der Pol oscillator Shinji Doi ∗ , Yohei Isotani, Ken-ichiro Sugimoto, Sadatoshi Kumagai Department of Electrical Engineering, Graduate School of Engineering, Osaka University, Osaka 565-0871, Japan Received 10 October 2002; accepted 25 February 2003 Communicated by C.R. Doering
Abstract The effect of noise on the period of the van der Pol oscillator forced by a sinusoidal signal is analyzed in detail. The dynamics of the forced nonlinear oscillator with noise is described by the evolution of probability densities by a Markov operator. The oscillator shows a 1 : 1 phase-locked oscillation whose period is equal to the period of the sinusoidal forcing when the forcing is sufficiently strong. Noise modulates the oscillator’s period as expected, but we show that this modulation is not uniform; the period of the oscillator changes drastically at a certain strength of noise. Thus this noise-induced breakdown of the phase locking can be considered as a stochastic bifurcation. The breakdown of stochastic phase-locking or the stochastic bifurcation is also analyzed from the viewpoint of the spectra of the Markov operator. 2003 Elsevier Science B.V. All rights reserved. PACS: 02.30.Oz; 02.50.Fz; 05.40.Ca Keywords: Nonlinear oscillator; Additive noise; Stochastic phase-locking; Stochastic kernel; Markov operator
1. Introduction The van der Pol oscillator is the representative model of a nonlinear oscillator and is widely utilized to represent various rhythmic phenomena in biology. Bonhoeffer–van der Pol or FitzHugh–Nagumo model [1] is the slight modification of the van der Pol oscillator and mimics qualitative behavior of a neuron. It is fundamentally important problem to examine how the period of the oscillator is modulated by the other oscillator or an external environment such as external noise [2]. * Corresponding author.
E-mail address: [email protected] (S. Doi).
Grasman and Roerdink [3] studied the van der Pol relaxation oscillator in the presence of noise and reduced the problem of examining the period of the oscillator to the analysis of the time necessary for a one-dimensional stochastic process to reach a boundary for the first time. Tateno et al. [4] extended this first-passage-time approach to the presence case of a periodic forcing and reduced the dynamics of the system to the evolution of probability densities by a Markov operator. Doi et al. [5] proposed the method to analyze phase-locking and bifurcation phenomena in the presence of noise using the spectra of the Markov operator. The present Letter advances this line of researches and investigates how the period of the sinusoidally
0375-9601/03/$ – see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(03)00386-4
408
S. Doi et al. / Physics Letters A 310 (2003) 407–414
forced van der Pol oscillator is changed by noise in detail. As the noise intensity increases, the oscillator phase-locked to the sinusoidal forcing is shown to change its period drastically from the period of the sinusoid at a certain strength of noise which can be considered as a breakdown point of a stochastic phase locking. We show that this breakdown point of a stochastic phase locking can be well characterized by the spectra of the Markov operator. The present Letter is organized as follows: Section 2 introduces the (piecewise-linearized) van der Pol oscillator forced by sinusoidal input in the presence of noise. Section 3 explains the mathematical framework to treat the noisy oscillator. Our method is based on the Markov operator which is a natural extension of a so-called Poincaré map or a return map in the noisefree (deterministic) case. Using this operator, the behavior of the system is described by the evolution of probability density functions of the phase of the sinusoidal input. Section 4 analyzes the mean period of the forced noisy oscillator in detail using the Markov operator and stochastic Arnold’s tongues (stochastic phaselocked regions) are presented in Section 5. The Letter concludes with some brief comments in Section 6.
2. van der Pol relaxation oscillator We consider a sinusoidally-forced (and piecewiselinearized) van der Pol relaxation oscillator in the presence of additive noise: dx = F (x, y), dt dy dW (t) = −x + v(t) + σ , dt dt 5 F (x, y) = y − x + |x + 1| − |x − 1| , 6 v(t) = A sin 2π(t/T + θ0 ) ,
(1)
Fig. 1. The x–y phase plane of (1), (2) in the limit of = 0 and the N-shaped x-nullcline (x˙ = 0). The y-nullcline (y˙ = 0) coincides with the y-axis. The closed orbit ABCD is a limit cycle. See text for details.
this case, the van der Pol oscillator shows a relaxation oscillation. Fig. 1 shows the x–y phase plane of Eqs. (1), (2). In the limit of = 0, the orbit with an initial value which is not on the N-shaped x-nullcline (x˙ = 0) instantaneously jumps to the x-nullcline in the horizontal direction since the horizontal velocity x˙ = F (x, y)/ becomes infinity unless F (x, y) = 0. Thus, all orbits are considered to move on the x-nullcline and both the sinusoid and noise terms of (2) modulate the velocity of the orbit along the x-nullcline. On the right (left) branch of the N-shaped x-nullcline, an orbit moves toward the point B (D, respectively). At the point B (D), the orbit instantaneously jumps to the point C (A, respectively), since we consider the limit of = 0. This limit is called singular since orbits are not differentiable at such jump points. The parallelogram ABCD is a limit cycle to which every orbit of (1), (2) approaches asymptotically.
(2)
where A, T and θ0 are the amplitude, the period and the initial phase of the sinusoidal input, respectively. W (t) is the standard Wiener process and σ dW (t)/dt denotes a Gaussian white noise with a noise intensity σ . Note that the cubic polynomial of the original van der Pol equation is piecewise linearized for the computational purpose [4]. Throughout, we consider a so-called singular (or relaxation) limit of = 0. In
3. Stochastic kernel and Markov operator Consider the noise-free (σ = 0) case and a Poincaré map or a return map as follows. Suppose that a state point starts at the point A of the x–y phase plane with an initial ‘phase’ θ0 of the sinusoidal input and that the point returns to the point A again after a time t. Then the phase of the sinusoidal input is changed from θ0 to θ1 : def
θ1 = p(θ0 ) = t/T + θ0 (mod 1).
S. Doi et al. / Physics Letters A 310 (2003) 407–414
Phase lockings of the noise-free case are usually analyzed by this return map p(θ ) and by the onedimensional mapping: θn+1 = p(θn ),
n = 0, 1, 2, . . . .
(3)
Note the sequence {θn } generated by this one-dimensional discrete-time dynamical system will be referred to as the orbit of (3) in the same way as the orbit of (1), (2) unless any confusions occur. Generally, in the case of m : n locking, the orbit or sequence {θn } asymptotically approaches to an n-periodic orbit (sequence): θ (0) , θ (1) , . . . , θ (n−1) , θ (0) , θ (1), . . . . We extend the deterministic map to the case with noise in the following way. In the noisy case, both variables θ0 and θ1 fluctuate owing to noise and thus are random variables Θ0 and Θ1 , respectively. Define a kernel function g(θ0 , θ1 ) using a conditional probability density function:
where the summation is taken for all ti such that θ1 = ti /T + θ0 (modulo 1). The first-passage-time pdf f0 (t, θ ) can be numerically calculated without simulations of the stochastic differential equations (1) and (2) [7]. Thus the kernel function can be obtained numerically also [4]. The function g takes relatively high values along the graph of the map p(θ0 ) of noisefree case; the values of g(θ0 , p(θ0 )) is much greater than other values of g(θ0 , θ1 ) if the noise intensity σ is not so large. Thus g can be considered as the stochastic extension of the (deterministic) return map. Using the kernel function g, we extend the system (3) to the noisy case. Let S denote a unit interval [0, 1] and D the set of absolutely-integrable non-negative functions with a unit L1 norm on S. A function which belongs to D is called a probability density function (pdf) or simply a density function. An (Markov [6]) operator P on D is defined by: Ph(θ ) =
g(θ0 , θ1 ) dθ1 = Pr{θ1 Θ1 θ1 + dθ1 | Θ0 = θ0 }.
f (t, θ ) =
f0 (s, θ )f0 t − s, θ ds,
0
θ = s/T + θ (mod 1),
(5)
since the phase of the sinusoid at the point C becomes θ after the elapse of time s on the branch AB of Fig. 1. The kernel function g(θ0 , θ1 ) of Eq. (4) is obtained by the function f (θ, t) as follows: fp (ti , θ0 ) × T , g(θ0 , θ1 ) = (6) θ1 =ti /T +θ0 (mod 1)
g(θ0 , θ )h(θ0 ) dθ0 ,
h ∈ D.
(7)
S
(4) Let the state point of the stochastic differential equations of (1), (2) start at the point A of Fig. 1 with the initial phase θ of the sinusoidal signal. We consider the time interval when the state point reaches the point B for the first time after leaving the point A. Let f0 (θ, t) be the probability density function (pdf) of this first-passage time from the point A to B. Owing to the symmetry of the trajectories of Fig. 1, the pdf of the first-passage time from the point C to D is same as f0 (θ, t). Thus the pdf f (θ, t) of the one-cycle time when the state point returns to the point A for the first time after leaving the point A with the initial phase θ of the sinusoid can be denoted by t
409
Let h0 (θ ) ∈ D denote the probability density function of the initial phase Θ0 when an orbit starts at the point A of the x–y phase plane. Then the density function of Θ1 when the orbit returns to the point A again is obtained by h1 (θ ) = Ph0 (θ ). Thus, the deterministic mapping (3) is extended to the system hn+1 (θ ) = Phn (θ ),
n = 0, 1, 2, . . . .
(8)
Note we investigate the asymptotic behavior of the sequence {hn (θ )} of probability density functions rather than the sequence {θn }. Let’s list several preliminary definitions [6]. A function h∗ (θ ) is called as the invariant density function of an operator P if the relation Ph∗ = h∗ holds. The invariant density is asymptotically stable if for any initial density function h0 ∈ D lim P n h0 − h∗ = 0.
n→∞
(9)
As is easily seen from its definition, the function g has the property g(θ0 , θ1 ) 0,
g(θ0 , θ1 ) dθ1 = 1 S
410
S. Doi et al. / Physics Letters A 310 (2003) 407–414
and is called a stochastic kernel. The inequality (10) inf g(θ0 , θ1 ) dθ1 > 0 θ0
S
is assumed to hold [4]. The operator with this property is known to have a unique asymptotically stable invariant density (see Corollary 5.7.1 of Ref. [6]). Thus the sequence {hn (θ )} produced by the operator P always approaches to a unique invariant density asymptotically as n → ∞. Note that the operator P has always a unique eigenvalue unity; the corresponding eigenfunction is the invariant density h∗ . The moduli of the other eigenvalues are less than unity. In particular, the modulus of the second (in the order of modulus) eigenvalue dominates the convergence speed of the initial pdf h0 to the invariant density h∗ in (9); the smaller the modulus of the second eigenvalue is the faster the convergence is.
4. Effects of noise on the period of the forced oscillator In this section, we investigate how the the forced oscillator (1), (2) varies when intensity σ is changed. The mean period oscillator is calculated by: T =
∞1
f (t, θ )h∗ (θ ) dθ dt,
period of the noise T of the
(11)
0 0
where the functions f and h∗ are the pdf of onecycle time and the invariant density (the stationary or asymptotic pdf of the phase) of the operator P, respectively. Fig. 2(a-i), (b-i), . . . , (f-i) show the mean periods numerically computed using (11) as a function of noise intensity σ . Panels (a-ii), (b-ii), . . ., (f-ii) show the moduli (absolute values) of the second and third eigenvalues of the operator P in the order of modulus (the first eigenvalue is always unity). Panels (a-i) and (a-ii) correspond to the case of A = 0: the absence case of the sinusoidal input. In this case, the mean period slightly decreases as the noise intensity σ is increased; noise accelerates the oscillator in the absence of the sinusoidal forcing. The inherent period of the van der Pol relaxation oscillator without both noise and sinusoidal input is
2(ln 7 − ln 3) ≈ 1.695 to which value the graph of the mean period approaches as the noise decreases. The moduli of second and third eigenvalues also decreases monotonically as noise increases. The plots of both second and third eigenvalues are completely overlapped; they are complex conjugate. Panels (b-i) and (b-ii) are the case that the amplitude and period of the sinusoid are A = 0.6 and T = 1.3, respectively. In this case, the sinusoidal input is not so strong that the oscillator could be phase-locked to the input when the noise is absent. The forced oscillator exhibits a so-called quasi-periodic oscillation without noise. The behavior of eigenvalues against noise is almost same as the previous case (a) while the mean period change is a little bit different from the previous case; the mean period does not change much as the noise intensity is increased. Although the oscillator is not phase-locked to the input, the input has much effect to the oscillator; the period of the oscillator is much shorter than the oscillator’s inherent period. As the noise intensity increases, the noise effect dominates over the effect of the sinusoidal input and the modulated period is finally growing up toward the inherent period (in the case of noise presence) shown in panel (a-i). Panel (c) corresponds to the case when the amplitude of the input is larger (A = 0.8). In this case, a 1 : 1 phase-locking appears when noise is absent. The mean period of (c-i) is nearly same as the period of the sinusoidal input in the range of small noise. Even if the noise intensity is increased, the mean period still stays around the input period. If the noise intensity is further increased, the mean period suddenly begins to increase. (Note that the inherent period of the van der Pol oscillator is about 1.7.) A striking fact is that the increase of the mean period by noise is not uniform; the mean period remains around the input period and then suddenly begins to increase. From this observation, we may say that the phase locking survives even when noise is added and the stochastic phaselocking breakdowns not gradually but critically at a certain noise intensity. The corresponding eigenvalues also present a critical behavior (see panel (c-ii)). The second and third eigenvalues are real numbers in the range of weak noise (this is a sign of stochastic phaselockings, see Ref. [5] for detail) and change their values from real to complex conjugate near the noise intensity σ = 0.065.
S. Doi et al. / Physics Letters A 310 (2003) 407–414
411
Fig. 2. Panel (i): mean periods of the forced van der Pol oscillator with noise are plotted as a function of the noise intensity σ . Panel (ii): moduli of 2nd and 3rd eigenvalues of the operator P are plotted as a function of noise intensity σ . (a) The absence case of the sinusoidal input: A = 0. (b) T = 1.3, A = 0.6. (c) T = 1.3, A = 0.8. (d) T = 1.3, A = 1.0. (e) T = 1.8, A = 0.2. (f) T = 1.8, A = 0.4.
Comparing the panels (c-i) and (c-ii), we can see that the noise intensity where the mean period starts increasing coincides nearly with the noise intensity where the eigenvalues change their values from real to complex. Thus we may say that the point of the eigenvalue change is the breakdown point of a stochastic phase-locking or the stochastic bifurcation point. Note that the second largest eigenvalue of the operator P dominates the convergence speed of the pdf sequence {hn (θ )} to the invariant density h∗ (θ ); the smaller the modulus of the eigenvalue is, the
faster the convergence becomes. Usually, if the noise intensity increases, the convergence to the stationary distribution becomes faster and the modulus of the second eigenvalue decreases toward zero. In the both panels (a-ii) and (b-ii), the moduli of the second and third eigenvalues are very close to unity and thus the convergence speed to the invariant density is very small in the range of small noise. This corresponds to the fact that the phase dynamics of the forced oscillator (in the case of noise-free or small noise) is similar to the simple rotation on the unit circle and the effect of phase spread-out to the invariant density is weak.
412
S. Doi et al. / Physics Letters A 310 (2003) 407–414
Differently from the both cases of (a-ii) and (b-ii), the modulus of the second eigenvalue in the small noise range is much less than unity in the case of (c-ii). This is because the nonlinearity of the dynamics in the case (c) is bigger than the case (a) and (b), thus the effect of phase spread-out to the invariant density caused by the nonlinearity is stronger than that of (a) and (b). Panels (d-i) and (d-ii) show the case of the stronger input. The mean period retains its value around the input period even if a very large noise is added. The transition point of both the mean period and the eigenvalue is shifted to a large noise direction. Also note that the moduli of the second eigenvalue are reduced totally than the case of (c) since the phase spread-out effect caused by the nonlinearity of the dynamics is increased by the increase of the amplitude of the sinusoidal input. Panels of (e) and (f) show the case of a different input period: T = 1.8 that is longer than the inherent
period of the van der Pol oscillator. Thus a large noise destroys the 1 : 1 phase-locking and decreases the mean period of the oscillator toward its inherent period (≈ 1.7). We can also see that the mean period changes its value drastically at a certain noise intensity and that the change point of eigenvalues is close to the change point of the mean period in both cases of (e) and (f).
5. Arnold tongue in the presence of noise In the above section, we have studied the noise effects on the forced nonlinear oscillator for several fixed values of the input amplitude and period. In this section we study the more global behavior on the input amplitude and period. Fig. 3(a) shows the Arnold tongues (boundaries of a phase-locked region) of noise-free (deterministic) forced van der Pol oscillator. Region of three phase-
Fig. 3. Arnold’s tongues. (a) Noise-free case (σ = 0). Regions of deterministic phase lockings of three types: 1 : 1, 2 : 1, and 5 : 3 lockings are shown. (b) Effect of noise on the 1 : 1 phase-locking. Noise-free and three cases of different noise intensities: σ = 0.02, 0.05, 0.09 are shown. (c) 2 : 1 phase-locking. Noise-free and the case of σ = 0.02 are compared. (d) 5 : 3 phase-locking. Noise-free and the case of σ = 0.02 are compared. See text for details.
S. Doi et al. / Physics Letters A 310 (2003) 407–414
locked patterns (1 : 1, 2 : 1, 5 : 3) are shown in the parameter space (plane) of the amplitude A and period T of the sinusoidal input. The 1 : 1 locking region occupies the most part of the parameter plane and the regions of 2 : 1 and 5 : 3 lockings are very narrow. Panel (b) shows both deterministic and stochastic Arnold’s tongues of the 1 : 1 phase-locking. The stochastic Arnold tongues are calculated as the points where the second and third eigenvalues change their values between real and complex values as mentioned in the above section. The widest region (the lowest plots) corresponds to the noise-free case and the stochastic phase-locking region shrinks slightly as the noise is increased. (The Arnold’s tongues of the σ = 0.02 case almost coincides with the noise-free case and are indistinguishable in this figure.) These calculations of the deterministic and stochastic Arnold’s tongues show the validity of our definition of the break-down point (boundaries) of stochastic phaselockings as the point of the eigenvalue change. Panel (c) is the case of 2 : 1 phase-locking: onecycle of the van der Pol oscillator is synchronized to the two cycles of the sinusoidal input. In this case, only noise-free and σ = 0.02 cases are shown since the stochastic 2 : 1 phase-locking is very fragile and the stochastic 2 : 1 locking region for stronger noise is very small or disappears. The method of computation of breakdown points (Arnold tongues) of 2 : 1 locking is same as the case of 1 : 1 locking. The last panel (d) shows both the deterministic and stochastic Arnold tongue of 5 : 3 lockings. By the same reason as the case (c), only noise-free and σ = 0.02 cases are computed. In this case, the points where the sixth and the seventh eigenvalues rather than the second and third ones change their values between real and complex are utilized to obtain the Arnold tongues, since the 5 : 3 phase-locking corresponds to the ‘threeperiodic’ behavior of the Markov operator while the 1 : 1 locking corresponds to the ‘one-periodic’ behavior of the operator (see Ref. [5] for details). From these figures we can see that the phaselocking regions (Arnold tongues) shrink by noise although this can be easily expected. We, however, stress that our numerical method which uses the Markov operator and its eigenvalues does not use numerical simulations of the stochastic differential equations (1), (2) and thus we can quantitatively study such a delicate behavior induced by noise in detail.
413
6. Discussion We have studied how the mean period of phaselocked oscillations in a sinusoidally-forced nonlinear oscillator is modulated by an additive noise using a piecewise linearized van der Pol oscillator. The dynamics of the noisy forced oscillator was expressed in terms of the evolution of probability density functions by a Markov operator [4]. This Markov operator is a natural extension of the deterministic (noise-free) Poincaré mapping or a return mapping which is the standard methodology to analyze forced oscillators. It should be note that the operator can be obtained only numerically but heavy numerical (Monte Carlo) simulations of the stochastic differential equations are not necessary. The operator is obtained by solving some integral equations numerically and thus can be computed with a high accuracy. A similar (Markov) approach has been taken by Plesser and Geisel [8] for the analysis of the noisy integrate-and-fire (IF) neuron model. Note that the IF neuron model is essentially same as the piecewise linearized van der Pol oscillator in the relaxation limit and the application of our method to the analysis of IF neuron models is a direct and simple application. Using the Markov operator, it is shown that the mean period of the noisy forced oscillator is modulated by noise not uniformly but critically; the mean period remains near the phase-locked period (the period of the sinusoidal input) against noise and drastically starts changing its value at a critical noise intensity. This behavior of the mean period modulation infers the existence of phase-lockings even in the case of noise presence (although the phase is not locked but varied by the presence of noise in a strict sense). In fact, near the critical change point of the mean period, the eigenvalues of the operator also showed the critical change: alternation of eigenvalue between real and complex values. The breakdown of deterministic phase-lockings is in itself the (saddle-node) bifurcation phenomenon. Although we have already proposed the method to analyze stochastic saddle-node and period-doubling bifurcations of the forced van der Pol oscillator [5] and of a one-dimensional mapping [9] and this method has already been utilized for the analysis of various noisy systems [10], there are still some place for the debate on the validity of our definition of stochastic
414
S. Doi et al. / Physics Letters A 310 (2003) 407–414
bifurcations. In this sense, the result of the present Letter that showed the abrupt change of the mean period at a critical noise level and the apparent coincidence of eigenvalue change, may give some basis on the usefulness and validity of our analysis method. A more rigorous and detailed discussion on the reason of the apparent closeness of both the critical change point of mean periods and the eigenvalues is necessary for the future research.
[3]
Acknowledgement
[7]
This research was partially supported by the Japanese Ministry of Education, Science, Sports and Culture, Grant-in-Aid (#10650429).
[8] [9]
References [1] R. FitzHugh, Biophys. J. 1 (1961) 445; S. Doi, S. Sato, Math. Biosci. 125 (1995) 229. [2] W. Horsthemke, R. Lefever, Noise-Induced Transitions, Theory and Applications in Physics, Chemistry, and Biology, Springer-Verlag, Berlin, 1984;
[4] [5] [6]
[10]
F. Moss, P.V.E. McClintock (Eds.), Noise in Nonlinear Dynamical Systems Vols. 1–3, Cambridge Univ. Press, Cambridge, 1987; T. Kapitaniak, Chaos in Systems with Noise, World Scientific, Singapore, 1988. J. Grasman, Asymptotic Methods for Relaxation Oscillations and Applications, Springer-Verlag, New York, 1987; J. Grasman, J.B.T.M. Roerdink, J. Stat. Phys. 54 (1989) 949. T. Tateno, S. Doi, S. Sato, L.M. Ricciardi, J. Stat. Phys. 78 (1995) 917. S. Doi, J. Inoue, S. Kumagai, J. Stat. Phys. 90 (1998) 1107. A. Lasota, M.C. Mackey, Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics, Springer, Berlin, 1994. A. Buonocore, A.G. Nobile, L.M. Ricciardi, Adv. in Appl. Probab. 19 (1987) 784. H.E. Plesser, T. Geisel, Phys. Rev. E 59 (1999) 7008. J. Inoue, S. Doi, S. Kumagai, in: Proceedings of NOLTA’97, International Symposium on Nonlinear Theory and its Applications, Hilton Hawaiian Village, Hawaii, 1997, pp. 13–16; J. Inoue, S. Doi, S. Kumagai, in: Proceedings of NOLTA’98, Le Regent, Crans-Montana, Switzerland, 1998, pp. 871–874; J. Inoue, S. Doi, S. Kumagai, Phys. Rev. E 64 (2001) 056219. T. Tateno, J. Stat. Phys. 92 (1998) 675; T. Yamanobe, K. Pakdaman, Biol. Cybern. 86 (2002) 155.