Analysis of noise-induced transitions between spiking and bursting regimes in Hindmarsh–Rose neuron model Lev B. Ryashko Evdokia S. Slepukhina lev.ryashko@urfu.ru evdokia.slepukhina@urfu.ru Ural Federal University (Yekaterinburg, Russia) Abstract The stochastic dynamics of the Hindmarsh–Rose model of neuronal ac- tivity is studied. For the parametric zone of tonic spiking oscillations, it is shown that random disturbances transform the spiking dynamic regime into the bursting one. For a quantitative analysis of the noise- induced bursting, we suggest a constructive approach based on the stochastic sensitivity function technique and the method of confidence domains. It allows us to give a geometric description for a distribu- tion of random states around the deterministic attractors and estimate critical values for the noise intensity corresponding to the qualitative changes in stochastic dynamics. We show that the obtained estimations are in a good agreement with direct numerical simulations. 1 Introduction Random perturbations can considerably affect the behavior of dynamical systems. Neuronal models are very sensitive to noise, and even small stochastic fluctuations can lead to significant qualitative changes in properties of such systems. One of the important mechanisms of neuronal behavior is excitability. Under the external impulse, a neuron can generate either individual spikes (i.e. abrupt changes in the electrical potential across a cell’s membrane) or sequences of periodic spikes. A special type of neuronal activity is bursting: a mode when periodic spiking alternates with intervals of quiescence. [1]. The three-dimensional Hindmarsh–Rose (HR) [2] system is one of the simplest models representing the bursting neural activity. The deterministic HR model is fairly well studied. The detailed bifurcation analysis in different parametrical zones can be found in [3–5]. As for the stochastic case, the effects of noise on the system with periodic external impulse, such as stochastic resonance [6, 7] and coherence resonance [8] were studied. The original deterministic HR system describes a wide range of dynamic regimes, such as periodic oscillations of various types, oscillations zones with period doubling and adding, coexistence of several attractors, chaos. This article studies the effect of random disturbances on the dynamics of the HR model in the parametric zone of tonic spiking oscillations. We show that under noise the spiking dynamic regime transforms into the bursting one. This phenomenon is confirmed by changes of distribution of random trajectories in the phase space. The most detailed description of the stochastic attractors in terms of probability density function is given by Kolmogorov-Fokker-Planck (KFP) equation. However, a direct usage of this equation is very difficult technically Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org 306 even in simple cases, therefore various asymptotics and approximations were developed. For an approximation of KFP solutions, a well-known quasipotential method [9,10] and a stochastic sensitivity function technique [14–16] can be used. The stochastic sensitivity function (SSF) technique allows us to construct confidence domains that are simple and evident geometrical models for a spatial description of a configurational arrangement of random states around the deterministic attractors. In this article, for a quantitative analysis of the noise-induced phenomena in the HR model, we suggest a constructive semi-analytical approach based on the stochastic sensitivity function technique and the method of confidence domains. 2 Stochastic sensitivity analysis: theoretical background Consider a nonlinear system of stochastic differential equations: dx = f (x) dt + εσ(x) dw(t). (1) Here, x is n-vector function, w(t) is n-dimensional standard Wiener process, ε is a scalar parameter of noise intensity. Assume that a corresponding deterministic system (ε = 0) has a T -periodic solution x = ξ(t) with an exponentially stable phase curve Γ (limit cycle). Under stochastic disturbances, random trajectories of the system (1) leave the deterministic cycle Γ and form some bundle around it. A detailed probabilistic description of the stochastic trajectories in this bundle is given by Kolmogorov-Fokker-Planck (KFP) equation. In a steady regime, one can consider stationary probability density function ρ(x, ε) governed by the stationary KFP equation. To avoid well-known technical difficulties with the direct use of this equation, various asymptotics and ap- proximations are developed [11–13]. For the approximation of KFP solutions, a well-known quasipotential method [9, 10] and a stochastic sensitivity function (SSF) technique [14, 15] can be applied. Let Πt be a hyperplane that is orthogonal to the limit cycle at the point ξ(t) (0 ≤ t < T ). For this plane, in the neighborhood of the point ξ(t), a Gaussian approximation of the stationary probabilistic distribution can be written [14] as: (x − ξ(t))> W + (t)(x − ξ(t))   ρt (x, ε) = K exp − 2ε2 with the mean value mt = ξ(t) and the covariance matrix D(t, ε) = ε2 W (t). The stochastic sensitivity matrix W (t) is a unique solution of the boundary problem Ẇ = F (t)W + W F > (t) + P (t)S(t)P (t) (2) with conditions W (T ) = W (0), W (t)r(t) = 0. Here ∂f F (t) = (ξ(t)), S(t) = G(t)G> (t), G(t) = σ(ξ(t)), ∂x rr> r(t) = f (ξ(t)), P (t) = Pr(t) , Pr = I − . r> r The eigenvalues λi (t) and the eigenvectors vi (t) of the SSF matrix characterize the distribution of random states in the Poincare section Πt near the point ξ(t) of the cycle. The value M = max λ1 (t) is a useful characteristic of the cycle as a whole. We consider M as a stochastic [0,T ] sensitivity factor of the limit cycle Γ. SSF matrix allows to construct the confidence ellipse with the center in point ξ(t). The equation of this ellipse in plane Πt looks like (x − ξ(t))> W + (t)(x − ξ(t)) = 2k 2 , where the parameter k determines a fiducial probability P = 1 − e−k . A set of these ellipses for t ∈ [0; T ) specify some confidence torus around a deterministic cycle. This torus is a confidence domain in a phase space for the stochastic cycle as a whole [16]. 307 3 Deterministic Hindmarsh–Rose system Consider the 3D Hindmarsh–Rose (HR) model [2]: ẋ = y − x3 + 3x2 + I − z ẏ = 1 − 5x2 − y (3) ż = r(s(x − x0 ) − z), where x is a membrane potential, variables y, z describe ionic currents, I is an external current; 0 < r  1 is a time scale parameter; s, x0 are other parameters. In this paper, we fix r = 0.002, s = 4, x0 = −1.6. The dynamics of the system (3) is considered under variation of the parameter I. Fig. 1 shows a bifurcation diagram of the deterministic system on the interval I ∈ (1.2, 4.0). Here, z- coordinates of equilibrium points (green solid for stable equilibrium, red dashed for unstable one), and z- coordinates of points of intersection of limit cycles and chaotic attractors with Poincare section plane x = 0 (blue solid) are plotted. The system (3) has a single equilibrium which is stable for I < I1 ≈ 1.288. It loses stability due to the subcritical Andronov–Hopf bifurcation when the parameter I passes from left to right through the point I1 . For I < I0 ≈ 1.2677 the stable equilibrium is a unique attractor of the system. As the parameter passes the point I0 , a stable limit cycle appears as a result of a saddle-node bifurcation. Such type of cycle is termed ”burst”. The bursting activity is observed in the system for I0 < I < I2 ≈ 3.292. For I0 < I < I1 the system exhibits a coexistence of the stable limit cycle and the stable equilibrium. Initially, near I0 , a limit cycle emerges with one spike in a burst. Then, with the increase of I, a sequence of period-adding bifurcations occurs: bursts with 2,3,4, ..., 13 spikes appear. A transition from burst with 12 spikes to burst with 13 spikes is accompanied with the transition to chaos. For I > I2 ≈ 3.292, a transition from the bursting chaos to the spiking chaos occurs. With further increase of I, the intermittency of chaotic and regular spiking is observed, then, at I ≈ 3.316, the transition from the chaotic spiking regime to regular one occurs. A sequence of backward period doubling bifurcations takes place: 4-cycle transforms to 2-cycle at I ≈ 3.328, then a transition from 2-cycle to 1-cycle occurs at I ≈ 3.364. Finally, the limit cycle transforms to the stable equilibrium due to the supercritical Andronov–Hopf bifurcation at I ≈ 25.261. Fig. 2 shows the examples of typical regular bursting (Fig. 2a) and tonic spiking (Fig. 2b) attractors. In this paper, we focus on the zone of transition from the bursting regime to the spiking one: I ∈ (3.0, 4.0). Figure 1: Bifurcation diagram: z-coordinates of stable (green solid) and unstable (red dashed) equilibrium points, and z-coordinates of points of intersection of limit cycles and chaotic attractors with x = 0 plane (blue solid). 308 a) x z z 1.5 3.5 3.4 3.4 1 3.3 3.3 3.2 0.5 3.1 3.2 3 0 3.1 2.9 −0.5 3 y −1 −2 x 2.9 −4 1 0.5 −6 0 −8 −0.5 0 200 400 600 800 t 0 200 400 600 800 t −1 b) x z z 1.5 3.5 3.4 1 3.5 3.3 0.5 3.2 3.4 0 3.1 −0.5 3 y −1 −2 x 2.9 −4 1 0.5 −6 0 −8 −0.5 0 200 400 600 800 t 0 200 400 600 800 t −1 Figure 2: Deterministic attractors and time series x(t), z(t) for a) I = 3.24, b) I = 3.4. 4 Stochastic Hindmarsh–Rose system Consider now the stochastically forced HR model: ẋ = y − x3 + 3x2 + I − z + εẇ, ẏ = 1 − 5x2 − y (4) ż = r(s(x − x0 ) − z), where w is a standard Wiener process with E(w(t) − w(s)) = 0, E(w(t) − w(s))2 = |t − s| and ε is a noise intensity. Let us study the effect of random disturbances on the system (4) in the zone of tonic spiking oscillations. Consider the value I = 3.7. Here, the spiking cycle with one spike per period is the attractor of the determin- istic system. Fig. 3 shows random trajectories starting from this deterministic cycle (in the projection on xOz plane) and corresponding x(t) plots for two values of noise intensity. For a sufficiently small noise (ε = 0.01), random trajectories are localized in some small vicinity of the deterministic limit cycle. The type of oscillations remains spiking, with one spike per period (Fig. 3a). This can be considered as noisy tonic spiking oscillations. If the noise intensity is greater than some threshold value (ε = 0.1), random trajectories go far from the determin- istic cycle, with a sharp decrease of z-coordinates, and then an oscillatory transition process of approaching to the deterministic cycle is observed (Fig. 3b). On x(t) plot, one can observe clearly defined intervals of quiescence and spiking phases. This indicates that under noise, the type of oscillations changed from spiking to bursting. Fig. 4 shows how noise changes the dispersion of random trajectories for I = 3.7. Consider Poincare section plane x = x̄, where x̄ is x-coordinate of the equilibrium. In Fig. 4a, z-coordinates of points of intersection of random trajectories with this Poincare section are plotted. For small noise intensities, random states concentrate near the deterministic limit cycle and have a sufficiently small dispersion. With an increase of noise, a significant deviation of z-coordinates in the direction of smaller values is observed. This corresponds to the emergence of bursting oscillations. Fig. 4b demonstrates the details of the distribution of random trajectories. The probability density function P (z) for z-coordinates of points of intersection of random trajectories with x = x̄ is plotted for I = 3.7 for 309 a) z x z 3.8 3.7 1 3.7 3.6 0 3.5 3.5 3.4 −1 3.3 3.3 −1 0 1 x 0 200 400 600 800 t 0 200 400 600 800 t b) z x z 3.8 3.7 1 3.7 3.6 0 3.5 3.5 3.4 −1 3.3 3.3 −1 0 1 x 0 200 400 600 800 t 0 200 400 600 800 t Figure 3: Stochastic generation of bursting oscillations for I = 3.7: a) ε = 0.01, b) ε = 0.1. various noise intensities. For a weak noise, peaks of P (z) are located above the limit cycle. For greater noise intensity values, one can observe a shift of the peak to the zone with smaller z-coordinates. Note that for bursting oscillations, random states with x < −1 are observed (the quiescence zones in Fig. 3b). The value x = −1 can be used as a threshold that separates in the phase space the quiescence from the spiking phase. For a quantification of the weight of spiking time in the total time of the observation, consider T the numerical characteristic ζ = Tq , where Tq is the time spending by the system in the region x < −1 and T is the total time. In Fig. 5, functions ζ(ε) for different I are plotted. For a small noise, random states are concentrated near the deterministic cycle, so ζ = 0. For greater noise intensity values, the bursting oscillations are observed, and ζ increases. Using the function ζ(ε), we can make empirical estimations for critical values of noise intensity corresponding to the transition from the spiking regime to the bursting one. For I = 3.5 we estimate ε∗ ≈ 0.006, for I = 3.7 we have ε∗ ≈ 0.02 and for I = 3.9 we get ε∗ ≈ 0.04. a) b) z P ε = 0.01 3.8 ε = 0.02 3.7 ε = 0.03 ε = 0.05 3.6 ε = 0.1 3.5 3.4 3.3 −3 10 10 −2 ε 0 3.3 3.4 3.5 3.6 3.7 3.8 z Figure 4: Distribution for points of intersection of random trajectories with Poincare section x = x̄ in dependence of noise intensity for I = 3.7: a) z-coordinates of points, b) probability density functions. 310 ζ I = 3.9 I = 3.7 0.4 I = 3.5 I = 3.4 0.3 0.2 0.1 10 −3 10 −2 10 −1 ε Figure 5: Function ζ. The emergence of noise-induced bursting oscillations can be explained by peculiarities of the phase portrait of the deterministic system. Trajectories starting close to the deterministic cycle, tend to it, but a character of a transient process depends on the value of the initial deviation. Indeed, if the deviations are relatively small, the trajectories tend to the cycle monotonically. If we take the initial deviations larger than some threshold, the trajectory goes sufficiently far from cycle and approaches to it making several loops that correspond to spikes (Fig. 6). Thus, considering different initial points, one can specify a border between these two different transient regimes in the phase space. Let us define this border by term pseudo-separatrix. To analyze the mechanism of stochastic generation of bursting oscillations, we apply the stochastic sensitivity function (SSF) technique [14, 15]. SSF allows us to approximate a dispersion of random states around the deterministic attractors. Eigenvalues and eigenvectors of the SSF define a geometric configuration of the confidence domains, i.e., the variance of random states around the stable attractor. For a stable limit cycle, eigenvalues and eigenvectors of the SSF depend on time along the limit cycle. They form a family of confidence ellipses around the stable limit cycle. The eigenvalues of SSF may vary nonuniformly along the limit cycle. This is shown on Fig. 7a, which displays the non-zero eigenvalues of the SSF matrix of the limit cycle for I = 3.7. One can observe a significant overfall of values along the cycle. The largest values of SSF correspond to the spiking phase of the limit cycle. Let us consider a point ξ of the limit cycle from ”transition region”, i.e. a part of the cycle from which the stochastic trajectories go off to the bursting zone in the phase space (see Fig. 7b). Note that the transition region is not a part of the cycle with maximal stochastic sensitivity, but it is a zone in which large sensitivity combines with the proximity to the pseudo-separatrix. z z 3.5 3.59 3.4 3.58 3.3 3.57 3.2 −1 0 1 x −0.93 −0.92 x Figure 6: Deterministic phase trajectories for I = 3.5. 311 a) b) λ1,2 z 1 10 3.78 0 3.77 10 x 3.76 1 10 −1 3.75 0 −2 y 0 10 −2 0 5 10 15 20 25 t −4 −1 Figure 7: Stochastic sensitivity of limit cycle (a) and transition region (b) for I = 3.7. Consider a hyperplane Π that is orthogonal to the cycle at the point ξ. For points of intersection of random trajectories with the plane Π, we can construct confidence ellipses corresponding to different values of noise intensity. For small noise, confidence ellipses give a good approximation for the dispersion of random states around the deterministic cycle (see Fig. 8a). Examining different initial points in the plane Π, one can define a pseudo-separatrix, specifying a border between two different transient regimes in the phase space. To determine the type of the transient regime (spiking of bursting), we use the threshold value x = −1: if a deterministic trajectory has points with x < −1, we assume that the type of the transient regime is bursting. Let us suggest a method allowing to estimate the critical value of the noise intensity ε∗ corresponding to the onset of noise-induced bursting. Fig. 8b shows a point of cycle in the transition region for I = 3.7, pseudo- separatrix, and confidence ellipses for two values of the noise intensity (a fiducial probability is P = 0.99). For the noise intensity less than critical value (ε = 0.01), a confidence ellipse is close to the deterministic cycle and do not intersect the pseudo-separatrix. With an increase of the noise intensity (ε = 0.02), the ellipse expands and intersects the pseudo-separatrix (see Fig. 8b). This means that with the high probability (P = 0.99), stochastic trajectories can go to the bursting zone of the phase space. The noise intensity that corresponds to the intersection of the confidence ellipse with pseudo-separatrix can be used as an estimation of the threshold value ε∗ . For I = 3.7, we get ε∗ ≈ 0.02. The obtained value is in a good agreement with the results of the a) b) ε = 5×10 −5 ε = 0.01 ε = 1×10−4 ε = 0.02 Figure 8: Confidence ellipses method for estimation of critical noise intensity for I = 3.7, P = 0.99: confidence ellipses (solid, dashed) and pseudo-separatrix (dash-dotted). 312 numerical simulations. Conclusion The effect of random disturbances on the Hindmarsh–Rose model in the parametric zone of tonic spiking os- cillations was studied. We showed that under noise, the spiking dynamic regime transforms into the bursting one. For a quantitative analysis of the noise-induced bursting, we suggested a constructive approach based on the stochastic sensitivity function technique and the method of confidence domains. It allows us to describe geometrically a distribution of random states around the deterministic attractors and estimate critical values for the noise intensity corresponding to the qualitative changes in stochastic dynamics. The obtained estimations are in a good agreement with direct numerical simulations. Acknowledgements The work was supported by the Government of the Russian Federation (Act 211, contract No. 02.A03.21.0006) and the Russian Foundation for Basic Research (project No. 16-31-00317 mol a). References [1] E. M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cam- bridge: MIT Press, 2007. [2] J. L. Hindmarsh, R. M. Rose. A model of neuronal bursting using three coupled first order differential equations. Proc R Soc Lond B Biol Sci, 221:87–102, 1984. [3] M. Desroches, T. Kaper, M. Krupa. Mixed-mode bursting oscillations: Dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster. Chaos, 23:046106, 2013. [4] G. Innocenti, A. Morelli, R. Genesio, A. Torcini. Dynamical phases of the Hindmarsh–Rose neuronal model: Studies of the transition from bursting to spiking chaos. Chaos, 17:043128, 2007. [5] A. Shilnikov, M. Kolomiets. Methods of the qualitative theory for the Hindmarsh–Rose model: A case study – a tutorial. Int J Bifurcation Chaos, 18:2141–2168, 2008. [6] A. Longtin. Autonomous stochastic resonance in bursting neurons. Phys Rev E, 55:868–876, 1997. [7] S. Reinker, E. Puil, R. M. Miura. Resonances and noise in a stochastic Hindmarsh–Rose model of thalamic neurons. Bull Math Biol, 65:641–663, 2003. [8] X. Shi, Q.-S. Lu. Coherence resonance and synchronization of Hindmarsh–Rose neurons with noise. Chinese Physics, 14:1088–1094, 2005. [9] M. Dembo, O. Zeitouni. Large deviations techniques and applications. Boston: Jones and Bartlett Publishers, 1995. [10] M. I. Freidlin, A. D. Wentzell. Random Perturbations of Dynamical Systems. New York: Springer, 1984. [11] C. Kurrer, K. Schulten. Effect of noise and perturbations on limit cycle systems. Physica D, 50:311–320, 1991. [12] B. Lindner, L. Schimansky-Geier. Analytical approach to the stochastic FitzHugh–Nagumo system and coherence resonance. Phys Rev E, 60:7270–7276, 1999. [13] G. N. Milshtein, L. B. Ryashko. A first approximation of the quasipotential in problems of the stability of systems with random non-degenerate perturbations. J Appl Math Mech, 59:47–56, 1995. [14] I. A. Bashkirtseva, L. B. Ryashko. Stochastic sensitivity of 3D-cycles. Mathematics and Computers in Simulation, 66:55–67, 2004. [15] I. Bashkirtseva, L. Ryashko. Sensitivity analysis of stochastic attractors and noise-induced transitions for population model with Allee effect. Chaos, 21:047514, 2011. 313 [16] L. Ryashko, I. Bashkirtseva, A. Gubkin, P. Stikhin. Confidence tori in the analysis of stochastic 3D-cycles. Mathematics and Computers in Simulation, 80:256269, 2009. 314