Stochastic generation of bursting oscillations in the three-dimensional Hindmarsh-Rose model Текст научной статьи по специальности «Математика»

УДК 517.917

Stochastic Generation of Bursting Oscillations in the Three-dimensional Hindmarsh-Rose Model

Lev B. Ryashko* Evdokia S. Slepukhina^

Institute of Mathematics and Computer Science Ural Federal University Lenina, 51, Ekaterinburg, 620083

Russia

Received 15.11.2015, received in revised form 14.12.2015, accepted 12.01.2016 We study the effect of random disturbances on the dynamics of the three-dimensional Hindmarsh-Rose model of neuronal activity. Due to the strong nonlinearity, even the original deterministic system exhibits diverse and complex dynamic regimes (various types of periodic oscillations, oscillations zones with period doubling and adding, coexistence of several attractors, chaos). In this paper, we consider a parametric zone where a stable equilibrium is the only attractor. We show that even in this zone with simple deterministic dynamics, under the random disturbances, such complex effect as the stochastic generation of bursting oscillations can occur. For a small noise, random states concentrate near the equilibrium. With the increase of the noise intensity, random trajectories can go far from the stable equilibrium, and along with small-amplitude oscillations around the equilibrium, bursts are observed. This phenomenon is analysed using the mathematical methods based on the stochastic sensitivity function technique. An algorithm of estimation of critical values for noise intensity is proposed.

Keywords: Hindmarsh-Rose model, neurodynamics, excitability, stochastic sensitivity, stochastic generation of bursting oscillations. DOI: 10.17516/1997-1397-2016-9-1-79-89.

Introduction

The three-dimensional Hindmarsh-Rose (HR) [1] system was developed for the modelling of bursting oscillations. Bursting is a special type of neuronal activity when intervals of periodic spiking (bursts) alternate with intervals of resting [2].

The analysis of the deterministic dynamics of the HR system in various parametric zones was provided in many works [3-8].

The stochastic variant of the HR model was considered by several authors. In particular, effects of noise on the system with periodic external impulse and stochastic resonance [9-12] were studied. Coherence resonance in the stochastic HR model was analysed in [13,14].

Due to the strong nonlinearity, the original deterministic system exhibits diverse complex 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, where the stable equilibrium is a single attractor of the deterministic system. However, even in this seemingly simple case, under the noise, the HR model exhibits rather complex phenomenon of the stochastic generation of bursting. For a weak noise, random states

* Lev.Ryashko@urfu.ru 1 Evdokia.Slepukhina@urfu.r © Siberian Federal University. All rights reserved

concentrate near the equilibrium. If the noise intensity is greater than some critical value, random trajectories can go far from the stable equilibrium, and along with small-amplitude oscillations around the equilibrium, bursts are observed.

The Kolmogorov-Fokker-Planck (KFP) equation is a basic mathematical model for the theoretical analysis of stochastic dynamics. This equation gives the most detailed description of the stochastic attractors in terms of probability density function. However, a direct usage of this equation is very difficult technically even in simple situations, therefore various asymptotics and approximations were developed. For an approximation of KFP solutions, a well-known quasipotential method [15,16] and a stochastic sensitivity function technique [17,18] can be used.

The stochastic sensitivity function (SSF) technique [19] allows us to construct confidence domains that are simple and evident geometrical models for a spatial description of a configu-rational arrangement of random states around the deterministic attractors. In this article, the stochastic generation of bursting oscillations in the HR model is analysed using the stochastic sensitivity function technique and the method of confidence domains.

1. Deterministic model

Consider the three-dimensional Hindmarsh-Rose (HR) [1] model:

x = y — x3 + 3x2 + I — z,

y = 1 — 5x2 — y, (1)

Z = r(s(x — xo) — z),

where x is a membrane potential, variables y, z describe activation and inactivation of ionic currents, I is an external current; 0 < r C 1 is a time scale parameter; s, x0 are other parameters.

Here we fix r = 0.002, s = 4, x0 = —1.6.

Consider the dynamics of the system (1) under variation of the parameter I. Fig. 1 shows a bifurcation diagram of the deterministic system: z-coordinates of equilibrium points and z-coordinates of points of intersection of limit cycles with Poincare section plane x = 0 in dependence of the parameter I value.

The system (1) has a single equilibrium which is stable for I < I\ ~ 1.288. It loses stability due to a subcritical Andronov-Hopf bifurcation when the parameter I passes from left to right through the point I\.

For I < I0 ~ 1.2677 the stable equilibrium is a unique attractor of the system. When the parameter passes the point I0, a pair of limit cycles (stable and unstable) appears as a result of a saddle-node bifurcation. With the increase of I, the unstable cycle branches off from the stable cycle, decreases in size, and at I = I\ it merges with the equilibrium. The stable cycle of such type 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 [5,6,8].

Fig. 2 shows examples of the burst with one spike for I = 1.268 (Fig. 2a), and the burst with two spikes for I =1.28 (Fig. 2b).

For I2 < I < I6 ~ 25.261, the system exhibits tonic spiking solutions. For I2 < I < I3 « 3.37, a period-doubling cascade of bifurcations with a transition to chaos occurs.

The equilibrium becomes stable again at I4 « 5.398 due to a supercritical Andronov-Hopf bifurcation, and for I4 < I < I5 « 6.198 the stable equilibrium coexist with the stable spiking

Fig. 1. Bifurcation diagram: ¿-coordinates of stable (solid) and unstable (dashed) equilibrium points; ¿-coordinates of points of intersection of limit cycles with x = 0 plane (solid)

cycle. When the parameter passes the point /4, the equilibrium loses its stability due to a subcritical Andronov-Hopf bifurcation.

The limit cycle disappears at /6 where a supercritical Andronov-Hopf bifurcation takes place. For / > /6, the stable equilibrium is a unique attractor of the system.

2. Stochastic model. Generation of bursting oscillations

Consider the stochastic HR model:

X = y — x3 + 3x2 + / — Z + £W,

y = 1 — 5x2 — y, (2)

Z = r(s(x — xo) — 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.

In this paper, we study the effect of random disturbances on the system (2) in the parametric zone / < /0 « 1.2677, where the stable equilibrium is the only attractor of the deterministic system. For parameter values that are close to the bifurcation point /0, the deterministic system is excitable: a small perturbation away from the equilibrium, but larger than some threshold, can result in a large excursion before returning to the resting state. This corresponds to the generation of action potentials. Excitability is one of the essential properties of neuron models.

a)

b)

Fig. 2. Phase portraits (in projection to xOz plane) and time series for a) I = 1.268, b) I = 1.28 (e = 0): stable equilibrium (circle, dashed), stable limit cycle (solid), unstable limit cycle (dash-dotted)

Consider the value I = 1.2. In Fig. 3, random trajectories (in projection on the plane xOz) starting from the stable equilibrium, and corresponding time series for different values of the noise intensity are plotted. Under random disturbances, trajectories of the system (2) leave the stable equilibrium and form a stochastic attractor. For a relatively small noise (e = 0.03), states of this stochastic attractor are concentrated near the equilibrium (see Fig. 3a). When the noise intensity is greater than some threshold value, random trajectories can go sufficiently far from the stable equilibrium, and along with low-amplitude oscillations, bursts are observed (see Fig. 3b for e = 0.1). One can see that a number of spikes in bursts is random.

Consider a process of transition to the stochastic generation of bursting oscillations for different values of the parameter I in the zone I < I0. Fig. 4 shows the details of the distribution of random states in dependence on the noise intensity e for I =1.2 (Fig. 4a) and I = 1.25 (Fig. 4b). Here, z-coordinates of points of intersection of random trajectories with the Poincare section x = x (x is x-coordinate of the equilibrium) are plotted. As one can see, for small noise intensities, random states concentrate near the equilibrium. With an increase of noise, states with large z-coordinates appear, that confirms the emergence of large-amplitude trajectories. Note that for I = 1.25 large-amplitude oscillations are observed for lower noise intensity values than for I = 1.2, because the value I = 1.25 is closer to the bifurcation point I0.

Note that for large amplitude oscillations, random states with x > —1 are observed. The value x = —1 can be used as a threshold that separates in the phase space small-amplitude oscillations near the equilibrium from the spiking phase. For a quantification of the weight of

T

spiking time in the total time of the observation, consider the numerical characteristic n = , where T is the time spending by the system in the region x > —1 and T is the total time. In Fig. 5, functions n(£) for different I are plotted. For a small noise, random states are concentrated near the equilibrium, therefore n = 0. For increasing noise intensity values, the large-amplitude oscillations are observed, and n increases as well.

a)

b)

Fig. 3. Phase portraits (in projection on xOz plane) and time series for I =1.2: a) e = 0.03, b) e = 0.1

Figs. 4 and 5 allow us to estimate a critical value of the noise intensity, corresponding to the onset of bursting generation. For I = 1.2, we get e* « 0.06, and for I = 1.25 we have e* « 0.04.

The emergence of noise-induced large-amplitude oscillations can be explained by peculiarities of the phase portrait of the deterministic system in the zone I < I0. The equilibrium is stable in this parametric region. If initial states are deviated from the equilibrium, the trajectories 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 equilibrium monotonically. If we take the initial deviations larger than some threshold, the transient process has a spike (see Fig. 6). The further increase in the deviations result in the trajectories with two, three and more spikes.

To analyze the mechanism of stochastic generation of bursting oscillations, we apply the stochastic sensitivity function (SSF) technique [17,19]. An essential theoretical background of SSF technique is given in the Appendix.

The stochastic sensitivity matrix reflects the geometry of bundles of stochastic trajectories around the deterministic attractors. If the attractor is an equilibrium, this distribution can be estimated by a confidence ellipsoid. Eigenvectors of the stochastic sensitivity matrix determine

Fig. 5. Function

a spatial arrangement of the ellipsoid while eigenvalues define a size of it.

Fig. 7 shows eigenvalues of the stochastic sensitivity matrix for the equilibrium in the parametric zone I € (1.2,1.288). One can see that the stochastic sensitivity increases unlimitedly with the approaching to the bifurcation point I\ where the equilibrium loses stability.

Random states of the system (2) are distributed non-uniformly around the equilibrium. This non-uniformity can be explained by the significant difference in eigenvalues of the stochastic sensitivity matrix for each fixed I. Thus, the confidence ellipsoids are elongated in the direction of the eigenvector v* corresponding to the largest eigenvalue of the stochastic sensitivity matrix, and they are very narrow in other directions. Thus, the eigenvector v* localizes the main direction for deviations of random trajectories from the equilibrium.

Fig. 8 shows the projection of the line passing through the equilibrium point in the direction v* on the plane xOz, and the intervals MM0, MMi, MM2, starting from the equilibrium M and corresponding to zones with different number of spikes in the transient process. Trajectories starting from any point of MM0 tend to the equilibrium monotonically without spikes. Initial points taken from MMi and MM2 result in large-amplitude trajectories with one and two spikes respectively.

Using these intervals of deviations, we can construct confidence intervals for the distribution of random trajectories, and estimate the critical values of the noise intensity. The value e0 corresponds to MM0, and for e < e0 there are small-amplitude stochastic oscillations near

10"3t-1-1-1-1-1--

1.05 1.1 1.15 1.2 1.25 I

Fig. 7. Eigenvalues of the stochastic sensitivity matrix for the equilibrium

the equilibrium. The values ei and e2 that can be estimated using intervals MM1 and MM2, correspond to the onset of the generation of bursts with one and two spike respectively.

For I = 1.2, we get e0 k 0.0675, e1 k 0.0684, e2 k 0.1084 and for I = 1.25 we have e0 k 0.0391, e1 k 0.0395, e2 k 0.0628. The obtained values are in a good agreement with the results of the direct numerical simulations.

Fig. 8. Intervals corresponding to zones with different number of spikes in the transient process for I = 1.2

Conclusion

The effect of random disturbances on the 3D Hindmarsh-Rose model was studied. We considered the parametric zone where a stable equilibrium is the single attractor in the deterministic system. We found that for sufficiently small noise, random states concentrate near the stable equilibrium, but if the noise intensity is greater than some threshold value, along with small-amplitude oscillations around the equilibrium, large-amplitude oscillations are observed. Thus, the stochastic generation of bursting oscillations occurs. In this paper, we provide the analysis of this phenomenon using the stochastic sensitivity technique, and suggest the method of estimation of critical values for noise intensity.

Appendix

Stochastic differential equation is a standard mathematical model for dynamical systems with random disturbances. Consider a nonlinear system of stochastic differential equations:

dx = f (x) dt + ea(x) dw(t). (3)

where x is n-vector function, w(t) is m-dimensional standard Wiener process, f (x) and a(x) are sufficiently smooth functions of appropriate dimensions, e is a scalar parameter of noise intensity. The stochastic HR system (2) is a special case of the system (3) for n = 3, m = 1,

' y — x3 + 3x2 + I — z ' 1

f = 1 — 5x2 — y , a = 0

r(s(x — xo) — z) 0

Under stochastic disturbances, random trajectories of the system (3) leave the deterministic attractor and form some probabilistic distribution around it. The detailed description of this distribution is given by Kolmogorov-Fokker-Planck (KFP) equation [21]. If a character of the transition process is not significant, and the main interest is a steady regime, one can consider stationary probability density function p(x, e) governed by the stationary KFP equation. In general case, for systems with n > 2, to solve this equation analytically can be very difficult. Therefore, various asymptotics and approximations are developed [22-24]. For the approximation of KFP solutions, a well-known quasipotential method [15,16] and a stochastic sensitivity function (SSF) technique [17-19,25] can be applied.

Consider the case when the system (3) have an exponentially stable equilibrium x. The corresponding approximation of the quasipotential near x gives an asymptotic of the stationary

distribution of random states of the system (3) in the Gaussian form:

(x — X, W-1(x — x))

p(x, e) = K exp — -

2e2

For the exponentially stable equilibrium x, the matrix W is the unique solution of the equation

FW + WFT = -S, F = f (x), S = GGT, G = a(x).

dx

The matrix W connects the stochastic input (e2) and the stochastic output (a covariance matrix D): D(e) = e2W. Thus, the matrix W can be considered as a coefficient of stochastic sensitivity of the system to random disturbances. We will term it the stochastic sensitivity matrix. This matrix characterizes a spatial arrangement of the stationary distributed random states of the stochastic system (3) around the deterministic equilibrium x. Eigenvectors vi of the matrix W determine a spatial arrangement of a confidence ellipsoid, whilst eigenvalues Xi define its size. In the special case when the coordinates of random states are uncorrelated and the variances in all directions are the same, confidence surfaces are spheres.

Consider confidence intervals in the direction of vi. The ends of the confidence intervals are points of intersection of the confidence ellipsoid with a line passing in the direction of vi. Their

coordinates can be written as x12 = x±aivi, where ai = ^/2Xiek, the coefficient k is determined

3

by the three-sigma rule k = . The confidence intervals allow to estimate critical values of

*

* aJ noise intensity: e* = . .

i V2Xk

The work was supported by Government of the Russian Federation (Act 211, contract 02.A03.21.0006).

References

[1] J.L.Hindmarsh, R.M.Rose, A model of neuronal bursting using three coupled first order differential equations, Proc. Royal Soc. Lond, B, Biol. Sci., 221 (1984), no. 1222, 87-102.

[2] E.M.Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, Cambridge, MIT Press, 2007.

[3] X.-J.Wang, Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle, Phys. D, 63 (1993), no. 1-4, 263-274.

[4] J.M.Gonzalez-Miranda, Complex bifurcation Structures in the Hindmarsh-Rose Neuron Model, Int. J. Bifurcation Chaos, 17 (2007), no. 9, 3071-3083.

[5] G.Innocenti et al., Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos, Chaos, 17(2007), 043128.

[6] A.Shilnikov, M.Kolomiets, Methods of the qualitative theory for the Hindmarsh-Rose Model: A case study - A Tutorial, Int. J. Bifurcation Chaos, 18(2008), no. 8, 2141-2168.

[7] 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(2013), no. 4, 046106.

[8] R.Barrio et al., Macro- and micro-chaotic structures in the Hindmarsh-Rose model of bursting neurons, Chaos, 24(2014), no. 2, 023128.

[9] A.Longtin, Autonomous stochastic resonance in bursting neurons, Phys. Rev. E, 55(1997), no. 1, 868-876.

10] S.Reinker, E.Puil, R.M.Miura, Resonances and Noise in a Stochastic Hindmarsh-Rose Model of Thalamic Neurons, Bull. Math. Biol., 65(2003), no. 4, 641-663.

11] V.V.Osipov, E.V.Ponizovskaya, Multivalued stochastic resonance in a model of an excitable neuron, Phys. Lett. A, 271(2000), no. 3, 191-197.

12] J.Baltanas, J.Casado, Noise-induced resonances in the Hindmarsh-Rose neuronal model, Phys. Rev. E, 65(2002), 041915.

13] S.Xia, L.Qi-Shao, Coherence resonance and synchronization of Hindmarsh-Rose neurons with noise, Chinese Physics, 14(2005), no. 6, 1088-1094.

14] H.Gu et al., Experimental observation of the stochastic bursting caused by coherence resonance in a neural pacemaker, Neuroreport, 13(2002), no. 13, 1657-1660.

15] M.I.Freidlin, A.D.Wentzell, Random Perturbations of Dynamical Systems, New York, Springer, 1984.

16] M.Dembo, O.Zeitouni, Large deviations techniques and applications, Boston, Jones and Bartlett Publishers, 1995.

17] I.A.Bashkirtseva, L.B.Ryashko, Stochastic sensitivity of 3D-cycles, Mathematics and Computers in Simulation, 66(2004), no. 1, 55-67.

18] I.Bashkirtseva, L.Ryashko, Analysis of excitability for the FitzHugh-Nagumo model via a stochastic sensitivity function technique, Phys. Rev. E, 83(2011), no. 6, 061109.

19] I.Bashkirtseva, L.Ryashko, Sensitivity analysis of stochastic attractors and noise-induced transitions for population model with Allee effect, Chaos, 21(2011), no. 4, 047514.

20] I.I.Gikhman, A.V Skorokhod, Stochastic differential equations and its applications, Kiev, Naukova dumka, 1982 (in Russian).

21] C.W.Gardiner, Handbook Of Stochastic Methods For Physics, Chemistry, And The Natural Sciences, New York, Springer, 1983.

22] B.Lindner, L.Schimansky-Geier, Analytical approach to the stochastic FitzHugh-Nagumo system and coherence resonance, Phys. Rev. E, 60(1999), no. 6, 7270-7276.

23] C.Kurrer, K.Schulten, Effect of noise and perturbations on limit cycle systems, Phys. D, 50(1991), 311-320.

24] G.N.Milshtein, L.B.Ryashko, A first approximation of the quasi-potential in problems of the stability of systems with random nondegenerate perturbations, J. Appl. Math. Mech, 59(1995), no. 1, 47-56.

25] I.A.Bashkirtseva, L.B.Ryashko, Quasipotential method in the study of local stability of limit cycles to random perturbations, Izv. vuzov. Prikl. nelineinaya dinamika, 9(2001), no. 6, 104113 in (Russian).

Стохастическая генерация пачечных колебаний в трехмерной модели Хиндмарш—Роуз

Лев Б. Ряшко Евдокия С. Слепухина

В работе изучается воздействие случайных возмущений на динамику трехмерной модели Хи-ндмарш-Роуз нейронной активности. Благодаря сильной нелинейности даже исходная детерминированная система демонстрирует весьма разнообразные и сложные динамические режимы (периодические колебания разных типов, удвоение и добавление периода колебаний, сосуществование нескольких аттракторов, хаос). В данной статье 'рассматривается параметрическая зона, в которой единственным аттрактором является устойчивое равновесие. Показывается, что даже в этой зоне с простой детерминированной динамикой под влиянием случайных возмущений в системе может наблюдаться такое сложное явление, как стохастическая генерация пачечных колебаний. При малых шумах случайные состояния концентрируются вблизи устойчивого равновесия. При увеличении интенсивности шума траектории могут проходить далеко от равновесия и наряду с малоамплитудными осцилляциями вокруг равновесия наблюдаются пачечные колебания. Проводится анализ этого явления с помощью математических методов, основанных на технике функций стохастической чувствительности, и предлагается алгоритм оценки критических значений интенсивности шума, вызывающего пачечные колебания.

Ключевые слова: модель Хиндмарш-Роуз, нейродинамика, возбудимость, стохастическая чувствительность, стохастическая генерация пачечных колебаний.