The stochastic dynamics of the Hindmarsh{Rose model of neuronal activity 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 noiseinduced bursting, we suggest a constructive approach based on the stochastic sensitivity function technique and the method of con dence domains. It allows us to give a geometric description for 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. We show that the obtained estimations are in a good agreement with direct numerical simulations.
Random perturbations can considerably a ect the behavior of dynamical systems. Neuronal models are very sensitive to noise, and even small stochastic uctuations can lead to signi cant 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 di erent parametrical zones can be found in [3{5]. As for the stochastic case, the e ects 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 e ect 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 con rmed 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 di cult technically 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 con dence domains that are simple and evident geometrical models for a spatial description of a con gurational 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 con dence domains. 2
Consider a nonlinear system of stochastic di erential equations:
dx = f (x) dt + " (x) dw(t): 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 di culties with the direct use of this equation, various asymptotics and approximations 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: t(x; ") = K exp (x (t))>W +(t)(x 2"2 (t)) with the mean value mt = (t) and the covariance matrix D(t; ") = "2W (t). The stochastic sensitivity matrix W (t) is a unique solution of the boundary problem (1) (2) with conditions Here
W_ = F (t)W + W F >(t) + P (t)S(t)P (t)
W (T ) = W (0);
W (t)r(t) = 0:
( (t)); S(t) = G(t)G>(t); G(t) = ( (t)); r(t) = f ( (t)); P (t) = Pr(t); Pr = I rr> 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 con dence ellipse with the center in point (t). The equation of this ellipse in plane t looks like (x (t))>W +(t)(x (t)) = 2k2; where the parameter k determines a ducial probability P = 1 e k. A set of these ellipses for t 2 [0; T ) specify some con dence torus around a deterministic cycle. This torus is a con dence domain in a phase space for the stochastic cycle as a whole [16].
Consider the 3D Hindmarsh{Rose (HR) model [2]: x_ y_ z_ = y = r(s(x 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 x 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 2 (1:2; 4:0). Here, zcoordinates of equilibrium points (green solid for stable equilibrium, red dashed for unstable one), and zcoordinates 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 2 (3:0; 4:0).
= y 5x2
y x0)
z); where w is a standard Wiener process with E(w(t) w(s)) = 0; E(w(t) w(s))2 = jt sj and " is a noise intensity.
Let us study the e ect 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 deterministic 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 su ciently 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 deterministic 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 de ned 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 su ciently small dispersion. With an increase of noise, a signi cant 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 10−3 −1 −1 0 0 1 1 a) 10−2 x x x 1 0 −1 x 1 0 −1 0 200 400 600 800 0 200 400 600 800 t t 0 200 400 600 800 0 200 400 600 800 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 quanti cation of the weight of spiking time in the total time of the observation, consider the numerical characteristic = TTq , 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 di erent 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. 3.7 3.5 3.3 z 3.7 3.5 3.3 3.2 0.4 0.3 0.2 0.1 10−3
I = 3.9 I = 3.7 I = 3.5 I = 3.4 10−2
10−1
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 su ciently far from cycle and approaches to it making several loops that correspond to spikes (Fig. 6). Thus, considering di erent initial points, one can specify a border between these two di erent transient regimes in the phase space. Let us de ne 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 de ne a geometric con guration of the con dence 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 con dence 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 signi cant 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 o 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 3.59 3.58 3.57 x −0.93 −0.92 x
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 con dence ellipses corresponding to di erent values of noise intensity. For small noise, con dence ellipses give a good approximation for the dispersion of random states around the deterministic cycle (see Fig. 8a).
Examining di erent initial points in the plane , one can de ne a pseudo-separatrix, specifying a border between two di erent 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, pseudoseparatrix, and con dence ellipses for two values of the noise intensity (a ducial probability is P = 0:99). For the noise intensity less than critical value (" = 0:01), a con dence 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 con dence 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 numerical simulations.
The e ect of random disturbances on the Hindmarsh{Rose model in the parametric zone of tonic spiking oscillations 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 con dence 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.
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). [1] E. M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge: MIT Press, 2007. [2] J. L. Hindmarsh, R. M. Rose. A model of neuronal bursting using three coupled rst order di erential 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. E ect 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 rst 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 e ect. Chaos, 21:047514, 2011. [16] L. Ryashko, I. Bashkirtseva, A. Gubkin, P. Stikhin. Con dence tori in the analysis of stochastic 3D-cycles.
Mathematics and Computers in Simulation, 80:256269, 2009.