TRANSITORY SHIFT IN THE FITZHUGH--NAGUMO MODEL Текст научной статьи по специальности «Физика»

Russian Journal of Nonlinear Dynamics, 2018, vol. 14, no. 2, pp. 169-177. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd180202

MSC 2010: 34C15

Transitory Shift in the FitzHugh — Nagumo Model

K. E. Morozov

A nonautonomous analogue of the FitzHugh - Nagumo model is considered. It is supposed that the system is transitory, i.e., it is autonomous except on some compact interval of time. We first study the past and future vector fields that determine the system outside the interval of time dependence. Then we build the transition map numerically and discuss the influence of the transitory shift on the solutions behavior.

Keywords: FitzHugh-Nagumo model, transitory system, separatrix, limit cycles, attractors

1. Introduction

There is a great number of mathematical models that describe the single neuron dynamics. The most popular ones are the Hodgkin-Huxley type models that have become widespread due to biological plausibility. The Hodgkin-Huxley system for a single neuron [1] consists of four nonlinear differential equations and, therefore, it is inconvenient for numerical analysis of neuron networks that contain many cells. It is often better to consider a simplification of the Hodgkin -Huxley model (e.g., [1-3]). The most widespread simplification is the FitzHugh - Nagumo model suggested in 1961 by R. FitzHugh [2]. It has the form

X3 t

X = X ^ y lexti ^ ^

Ty = x — by + a.

In system (1.1), x stands for the membrane potential, y is a recovery variable, Iext is the parameter that determines the stimulus current, a, b, T are positive parameters. Values of a, b, T are defined by the inner structure of the system. For example, the approximate values obtained from the study of the giant squid are a = 0.7, b = 0.8, T = 1/0.08 [2, 3].

Received February 09, 2018 Accepted May 11, 2018

This work was supported in part by the Russian Foundation for Basic Research under grants no. 18-01-00306 and no. 16-01-00364.

Kirill E. Morozov kirwamath@gmail.com

Lobachevsky State University of Nizhny Novgorod prosp. Gagarina 23, Nizhny Novgorod, 603950, Russia

Let us rewrite system (1.1) in a more convenient form (with the notations retained). We make a linear shift y — y + a/b:

' x3

x = x---j-y + I, ^^

y = e(x - by),

where e = -¡L, I = Iext - a/b.

System (1.2) does not depend on time explicitly and, therefore, describes only steady-state regimes. Parameters b and e are inner characteristics of the biological object under consideration and we can often neglect their variations over a short period of time. Parameter I determines the magnitude of the external stimulus and it can be significantly changed, which affects the solutions behavior. This motivates us to consider the nonautonomous analogue of the FitzHugh-Nagumo system:

x = x-^-y + I(t), y = e(x - by).

Let us suppose that the function I(t) depends on time only on a finite number of intervals. Without loss of generality, one can assume that the function is time-dependent only on the interval [0;t], outside of which I(t) is constant. In addition, we shall suppose that I(t) is continuous. The concrete form of I(t) for t <E [0; t] depends on the shape of the external stimulus. For instance, it can have the following form:

I(t) = Io(1 - f(t))+ If(t), (1.4)

where f (t) is a cubic spline:

'0, t< 0,

f (t) = {(t/T)2(3 - 2(t/T)), 0 < t < t, (1.5)

1, t>T,

I0, Ii are some constants.

System (1.3) with the function I(t), which satisfies (1.4) and (1.5), is transitory (in the sense of the definition given in [4]). According to [4-6], f (t) is called a transition function. For t < 0 (1.3) describes an autonomous vector field in the plane called the past vector field GP (I(t) = I0, t < 0), for t > t it describes a vector field called the future vector field GF (I(t) = Ii, t > t). We are mainly interested in the effect of nonautonomy on the establishment of dynamics.

We shall first study the behavior of trajectories of the future and past vector fields under variation of parameters. These fields are determined by system (1.2) (the FitzHugh-Nagumo model). A lot of literature is devoted to the system (see [3, 7, 8]), there are educational manuals (e.g., [9]). We shall construct the partition of the (I, b)-plane into regions with topologically different phase portraits for fixed e. It has not been done. Then we shall study the nonautonomous system (Section 3).

2. Investigation of the autonomous system (1.2)

In this section we fix the value of e and study the behavior of system (1.2) under variation of (b, I). Notice that the change I -—I is equivalent to the change x — —x, y — —y. It leads to symmetry of the bifurcation diagram with respect to the axis I = 0. We shall further assume that I ^ 0.

Equilibria of (1.2) lie on the line x — by = 0 and the x-coordinate is defined by the equation

x3 + x - 3/ = 0.

If a = (V)S + |f>> 0, then ,,2) h» a unique eqnHftwiam, if Q < 0, then it has

three equilibria. The transition from one equilibrium to three occurs through the saddle-node bifurcation and Q is a control parameter in the normal form on the center manifold.

Let (xo, y0) be the coordinates of the equilibrium state. The linearization matrix in the neighborhood of (x0,y0) has the form

A = (1—x2 —1Y

\ e —be J

Then the characteristic equation for (x0,y0) takes the form

A2 — ¿A + A = 0,

where 5 = tr(A) = 1 — x2 — be, A = det(A) = e(1 — b(1 — x0)).

Note that A = e/b((1 — b)/b + x0)) > 0 in the case of a single equilibrium. Indeed, it is not hard to see that A vanishes if only Q = 0. In the domain determined by the inequality Q > 0 A preserves the sign and, therefore, it is positive. So, if system (1.2) has a unique equilibrium, then this equilibrium is a focus (node) that changes its stability when 5 changes the sign (via the Andronov-Hopf bifurcation). The Andronov-Hopf bifurcation occurs in the case of three equilibria as well (in the neighborhood of one of the foci).

Necessary conditions of the Andronov-Hopf bifurcation are determined by the relations 5 = 0, A > 0.

Let us compute the first Lyapunov coefficient L1. To de this, we first make a linear substitution in system (1.2) transferring the equilibria (x0,y0) to the origin:

(xnew — x x0,

ynew = y — V0.

Then we obtain

ix = (1 — x0)x — y — x0x2 — x3/3, • b M

y = ex — eby.

We can now apply the Bautin formula [10] and compute Li

L\ = —^—(26 — eb2 — 1).

4A3/2

Let us now determine the values of the parameters (b,I) that nullify Li. We have the system

1 - x0 = eb, (2.2)

eb2 < 1, (2.3)

2b - 1= eb2; (2.4)

r x0 , 1 - b

I = yH--(2-5)

Conditions (2.2) and (2.3) are conditions of a complex focus, (2.4) implies that Li = 0, (2.5) is the equation that defines the x-coordinate of the equilibrium.

It follows from (2.3) and (2.4) that b < 1. Further, we obtain from (2.4) that e < 1. Expressing b from (2.4) under these assumptions, we have

o* - £

Let us now use (2.2) and (2.5). Excluding x0, we finally obtain the value of I* > 0:

I* =4/3(1 - e)3/4.

To compute the second Lyapunov coefficient L2, we first make the following change in system (2.1)

y - (1 - x0)x

X —> x, y —> -—-,

y y/A

and normalize the time variable t' = \fAt. We have

xo 2 1 3

x = — y--— x--— X ,

y V^ 3A/A

y A 3A

Then L2 takes the form

L2 = ^¿s~2 i-16^ " + 3A(X " -To)) = ^~2 (nbef - 16(be)2 + 3A(be)).

Let us substitute the expressions for A and b* in this system

L2(e) = (16(1 - «3 - 16(1 - VT^e)2 + 6^1^7(1 - Vl^e)2).

Denote 1 - Vl -1 = t, 0 < t < 1,

Hence, L2 < 0 for 0 < t < 1. Thus, there are no more than two limit cycles born in the case of Li = 0. Moreover, the inner cycle is unstable, and the outer one is stable.

Thus, the following proposition holds.

Proposition. 1) for e ^ 1 there is a unique limit cycle born from a complex focus as a result of the Andronov - Hopf bifurcation;

2) for 0 < e < 1 there are two points of the parameter plane corresponding to the Bautin bi-

e

£

furcation (the generalized Andronov - Hopf bifurcation): (b,

The Bautin point is the point where the double cycle curve begins. This curve can be found numerically. For e = 1 the curve degenerates, and disappears for e > 1.

In addition to the local bifurcations, there are bifurcations of separatrix loops in the region defined by Q < 0, which can lead to the birth or disappearance of limit cycles. For a fixed value of e, curves corresponding to the separatrix loops can be found numerically with high accuracy.

Figure 1 shows the complete bifurcation diagram for e = 0.2 as an example. The red curves correspond to the Andronov - Hopf bifurcation, the black curves show the double cycle bifurcation,the green curves show the saddle-node bifurcation and, finally, the blue curves correspond to the global bifurcations of the separatrix loops. The blue lines are located at a distance of the order 10_4, but the calculation accuracy allows us to separate them. Figure 2 shows the enlarged fragments of the marked areas of the diagram. In the domain Q > 0 (one equilibrium) the form of the diagram is preserved for e £ (0,1). The qualitative phase portraits are shown in Fig. 3 for all values of parameters. Unstable limit cycles are colored red, stable ones are colored blue.

System (1.2) can be written in the form

ex = x - ^ -y + 1, y = (x - by).

If e is supposed to be small, then the phase space is divided into regions of fast and slow motions. This limiting case is well studied, the techniques of [12] are used in the study (see also [2, 3]).

1.2

X .53,1.13)

1 3

(2.24,0.i

13 -5 6

0.8

0.4

0.4

0.8

1.2

1.6

2.0

Fig. 1. Bifurcation diagram of (1.2) for e = 0.2.

©

(a)

0.0604

0.0602

0.0600

4

1 ">

lZ

13

(b) 2.5 1.5 0.5 -0.5 -1.5 -2.5

1.2634 xlO"3

1.2638

1.2642

3

Ç

10 N 6

7

\

1.365 1.367 1.369 1.371

1.373

Fig. 2. Enlarged fragments of the parameter plane. ©

©

Fig. 3. Qualitative phase portraits of (1.2) for e = 0.2.

3. Investigation of the nonautonomous system (1.3)

The nonautonomous system (1.3) is defined in the extended phase space R3. The dynamic system for the past vector field Gp is defined in the semispace R2 x R_ (in R2 x R+r (t > t) in the case of the future vector field GF).

If A is any invariant set of the past vector field GP, then we call the set {(x; y; t): (x, y) G A, t ^ 0} a backward-invariant set in the extended phase space. For instance, equilibria and periodic orbits of Gp correspond to backward-invariant sets for the full vector field. If A is any invariant set of the future vector field GF, then we call the set (x; y; t): (x,y) G A,t ^ t a forward-invariant set in the extended phase space. Limit cycles of the past and future vector fields correspond to cylindric integral surfaces in the extended phase space which are backward-invariant and forward-invariant, respectively. In a similar way one can determine backward-hyperbolic and forward-hyperbolic sets [4].

Since the nonautonomous portion of the dynamics is confined to a compact interval of time, then its influence on the solutions behavior can be described in terms of the transition map T : R2 ^ R2, defined as follows.

Definition. Let x(x0,t) be a solution of (1.3) that satisfies the following Cauchy condition x(x0,0) = xo. Then

T(xo) = x(xo, t). (3.1)

Let l(t),t ^ 0 be a semitrajectory of the autonomous field GP located in some basin of attraction. If T(l(0)) belongs to the attraction basin of the same type of attractor of the future vector field Gf, then we say that there is a preservation of the limit regime. Otherwise, there is a change of regime.

Since the transitory shift changes the value of I, the pairs of parameters (b, I) that correspond to the past and future vector fields have the form (b, I0) and (b, Ii).

Let e = 0.2 be fixed. Suppose that I0 = 0. We first consider b G (0,1). From the analysis of the bifurcation diagram (1) and phase portraits (3) it follows that GP has a phase space structure that corresponds to region 1 of (1), GF can have three topologically different structures of the phase space depending on the value of the parameter b,b < 1: 1) a stable limit cycle surrounding an unstable equilibrium (region 1); 2) a stable equilibrium and a stable limit cycle, separated by an unstable cycle (region 2); 3) a unique globally stable equilibrium (region 3).

In the first case, the self-oscillating regime is preserved, the time dependence does not affect the qualitative behavior of solutions. In the second case, there is a change of the regime for some solutions (the self-oscillatory regime changes to the stationary one). The domain of the phase space that corresponds to the change is determined by the unstable limit cycle of the future vector field and, consequently, does not depend on the transition map (the form of time dependence). The qualitative behavior of solutions changes for all solutions in the third case.

In the case of three equilibria (b > 1), more combinations of topological structures of the past and future vector fields are possible. However, some of them lead to qualitatively similar changes in the solutions behavior. It follows from the analysis of (1) that the field Gp can have a phase space structure that corresponds to the region 14 (trajectories of GP tend to a globally stable limit cycle), to the regions 7 or 10 (trajectories of GP tend to a stable limit cycle or to one of two stable equilibria), to the region 6 (trajectories of GP tend to one of two stable equilibria). In the first of these cases, one stable equilibrium state may appear in the future vector field, i.e., it is possible to change the self-oscillatory regime to the stationary one (e.g., if (b, Ii) belongs to the region 3). Domains of the phase space that correspond to the change are determined only

0.8

0.6

0.4

0.2

0

0.2

-0.4

-0.6

-0.8

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

Fig. 4. Images of limit cycles of the past vector field under T (I0 = 0, I\ = 0.5, b = 1.363, e = 0.2, t =1).

by invariant sets of the future vector field. In the second case, there may also be a change of the self-oscillatory regime to the stationary one and vice versa, but the domains of the phase space that determine the change depend on the form of the time dependence. Besides, it is possible here that solutions change one stationary state to another ("left" focus to "right" and vice versa, see Fig. 4). Finally, in the third case, the future vector field cannot have any stable limit cycles and the regime is preserved here.

In Fig. 4 the images of limit cycles of the past vector field under the transition map T are colored blue. Green domains correspond to the change of the self-oscillatory regime to the stationary one, the red ones correspond to the preservation of the (stationary) regime. Solutions mapped into the red domain change their stationary state ("left" focus to "right").

4. Conclusion

In the present paper we consider a generalization of the FitzHugh - Nagumo system, which is widely used in modeling the neuron activity. The generalization implies that, instead of an autonomous system, we have a nonautonomous one such that time dependence is confined on a compact interval of time. We call such systems transitory dynamical systems. For t < 0 the right-hand sides of the system under consideration define an autonomous vector field, called the past vector field, for t > t they define another autonomous vector field, called the future vector field.

The investigation of transitory systems started fairly recently [4-6] and the first examples were Hamiltonian transitory systems [4]. For these systems, the measure of transport between the cells of the past and future vector fields is introduced. Then, quasi-Hamiltonian systems were considered (the Duffing-Van der Pole equation [5], the Josephson equation [6]). In [4, 5], the question of quantitative characterization of a possible regime change was also considered.

In the present paper, we study a nonconservative system that is not close to an integrable one. Possible structures of the past and future vector fields were found for a fixed value of e. A bifurcation diagram was constructed. It was established that there are two limit cycles at most if the system has a unique equilibrium and three cycles at most if the system has three equilibria. We build the transition map numerically and then discuss the influence of time dependence on the regime establishment. In some cases, we can understand this influence

without knowing the form of the transition function. It is stated that the transitory shift can

lead to the appearance of spike oscillations.

Useful conversations with Morozov A.D. are gratefully acknowledged.

References

[1] Hodgkin, A. L. and Huxley, A. F., A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve, J. Physiol., 1952, vol. 117, no. 4, pp. 500-544.

[2] FitzHugh, R., Impulses and Physiological States in Theoretical Models of Nerve Membrane, Bio-phys. J., 1961, vol.1, no. 6, pp. 445-466.

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

[4] Mosovsky, B. A. and Meiss, J. D., Transport in Transitory Dynamical Systems, SIAM J. Appl. Dyn. Syst., 2011, vol. 10, no. 1, pp. 35-65.

[5] Morozov, A.D. and Morozov, K.E., Transitory Shift in the Flutter Problem, Nelin. Dinam., 2015, vol. 11, no. 3, pp. 447-457 (Russian).

[6] Morozov, A.D., Morozov, K.E. Transitory Shift in the Pendular Type Equations, Nelin. Dinam., 2016, vol. 12, no. 4, pp. 577-589.

[7] Rocsoreanu, C., Georgescu, A., and Giurgiteanu, N., The FitzHugh-Nagumo Model: Bifurcation and Dynamics, Dordrecht: Springer, 2000.

[8] Ringkvist, M. and Zhou, Y., On the Dynamical Behaviour of FitzHugh - Nagumo Systems: Revisited, Nonlinear Anal, 2009, vol.71, nos.7-8, pp. 2667-2687.

[9] Prokin, I.S., Simonov, A.U., and Kazancev, V. B., Mathematical Modelling of Neurodynamical Systems, Nighny Novgorod: Nighegorodsk. Univ., 2012 (Russian).

[10] Bautin, N.N., Behavior of Dynamical Systems near to the Boundaries of Stability, Moscow: Nauka, 1984 (Russian).

[11] Andronov, A.A., Leontovich, E. A., Gordon, 1.1., and Maier, A. G., Theory of Bifurcations of Dynamic Systems on a Plane, New York: Wiley, 1973.

[12] Mishenko, E. F. and Rozov, N.Kh., Differential Equations with Small Parameters and Relaxation Ocsillations, New York: Plenum, 1980.