ANALYTICAL PROPERTIES AND SOLUTIONS OF THE FITZHUGH--RINZEL MODEL Текст научной статьи по специальности «Математика»

Russian Journal of Nonlinear Dynamics, 2019, vol. 15, no. 1, pp. 3-12. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd190101

NONLINEAR PHYSICS AND MECHANICS

MSC 2010: 34A05,34A34

Analytical Properties and Solutions of the FitzHugh-Rinzel Model

A. I. Zemlyanukhin, A. V. Bochkarev

The FitzHugh-Rinzel model is considered, which differs from the famous FitzHugh-Nagumo model by the presence of an additional superslow dependent variable. Analytical properties of this model are studied. The original system of equations is transformed into a third-order nonlinear ordinary differential equation. It is shown that, in the general case, the equation does not pass the Painleve test, and the general solution cannot be represented by Laurent series. Using the singular manifold method in terms of the Schwarzian derivative, an exact particular solution in the form of a kink is constructed, and restrictions on the coefficients of the equation necessary for the existence of such a solution are revealed. An asymptotic solution is obtained that shows good agreement with the numerical one. This solution can be used to verify the results in a numerical study of the FitzHugh-Rinzel model.

Keywords: neuron, FitzHugh-Rinzel model, singular manifold, exact solution, asymptotic solution

1. Introduction

The fundamental experimental and theoretical work of Hodgkin-Huxley [1], in which the propagation of nerve impulses along the axon of a squid's nerve cell was studied, excited interest in mathematical modeling of wave modes in biologically active media. This model included four differential equations for variables describing the membrane potential and three types of potential-dependent ion channels, as well as eight auxiliary algebraic equations. The model satisfactorily predicted the response of a separate neuron to external stimulation, but for modeling two- and three-dimensional arrays of connected neurons it turned out to be too complicated.

Received November 28, 2018 Accepted March 05, 2019

Alexander I. Zemlyanukhin [email protected]

Andrey V. Bochkarev [email protected]

Gagarin State Technical University

ul. Politekhnicheskaya 77, Saratov, 410054 Russia

The reduction of this system, the FitzHugh - Nagumo (FN) model [2-4], consists only of two differential equations for fast and slow variables and is widely used for numerical simulation of the active media wave dynamics.

For example [5], in numerical experiments with two-dimensional self-oscillating active media described by the FN model equations, solutions were found and studied that correspond to a phase concave spiral wave, which rotates continuously around a circular obstacle in a finite-sized medium. A concave spiral wave is an anomalous phenomenon for reaction-diffusion processes. Such a wave has a doublet structure, when the main pulse of excitation in the active medium is followed by a second pulse of much smaller amplitude. In contrast to normal spiral waves, the doublet structure waves can be elastically reflected from each other and from the boundaries of the medium, thereby exhibiting soliton properties.

Replacing the constant coefficients of the FN model, which defines the threshold between electrical silence and electrical firing as well as the external current through the membrane, by periodic functions of time leads to a periodic FN model [6]. The authors of [7] proposed a design method based on bifurcation analysis to generate bursting responses in the FN model with a simple periodic external force.

In [8], the FN model is used to study trigger waves, which are a recurring biological phenomenon involved in transmitting information quickly and reliably over large distances. Well-characterized examples of trigger waves include action potentials propagating along the axon of a neuron, calcium waves in various tissues, and mitotic waves in Xenopus eggs. The richness of the dynamic modes of the FN model allowed examining different types of trigger waves — spatial switches, pulses, and oscillations — and to show how they arise. In the recent work [9] it is shown that apoptosis, which is an evolutionarily conserved form of programmed cell death, propagates through the cytoplasm as a trigger wave.

The analytical properties of the initial and perturbed FN models are investigated and classes of their exact and asymptotic solutions are constructed in [10, 11].

In the original version of the FN model, there is no burst generation mode during which periods of rapid action potential spiking are followed by quiescent periods. At the same time, spiking and bursting observed in nerve membranes seem to be important when one investigates information representation models of the brain [7]. To eliminate this drawback, the FN system was supplemented with a third equation for the superslow variable in the FitzHugh - Rinzel (FR) model [12]:

1 3

vt = v--v -w + y + Iext,

wt = 5 (v - pw), (L1)

yt = i (a - v - jy),

where v(t), w(t) and y(t) are, respectively, fast, slow and superslow variables and Iext, P, y, 5, i, a are constants defining the dynamic mode of the model. The subscript hereinafter means differentiation by the corresponding variable.

Currently, the FR model (1.1) is fairly well known. In [13] various networks of diffusively coupled identical neurons modeled by a system of FR coupled differential equations were considered. Synchronization conditions in a network in which the central element modeling the pacemaking neuron is linked to the group of uncoupled neurons were presented. Using the Poincare return mappings, the authors of [14] were able to examine in detail the bifurcations that underlie the complex activity transitions between: tonic spiking and bursting, bursting and mixed-mode oscillations, and finally, mixed-mode oscillations and quiescence in the FR model.

This paper is devoted to clarifying the analytical properties of the FR model. In the 2nd part of the paper, the original system of equations reduces to a single 3rd-order equation. In the 3rd part, its particular exact solution is constructed. In the 4th part, the exact solution of the original system is found. Finally, in the last part, an asymptotic solution is obtained and good agreement with the numerical solution is established.

2. Reducing the system to a single equation

We reduce system (1.1) to one equation for the fast variable v(t). To do this, we differentiate term by term the first equation of (1.1) with respect to t and replace the first-order derivatives vt, wt and yt with the right-hand sides of the corresponding equations of system (1.1). Repeat this operation with the resulting equation and replace the 2nd and 3rd lines of the original system with two new equations obtained in this way:

1 3

vt = v - -v° -W+y + /ext,

vtt = fi (v,w,y), (2.1)

Vttt = f2 (v,W,y).

The functions fi and f2 do not contain derivatives and are not given here because of their bulkiness. The first two equations of system (2.1) are linear in the functions (w,y) and can be rewritten in the form

W " y = f3 (^vt), (2.2) (/5 + v2 - 1) w - (yP + v2 - 1) y = f4 (v, vtt).

Assuming that the determinant of system (2.2)

A = -/5 + YP (2.3)

is nonzero, we find (w,y) and substitute the result in the 3rd equation of (2.1) to get:

vttt + {v2 + a - 1) vtt + 2v(vtf + (av2 + c) vt + ^bv3 + ev - (/ext7 + <?) P»S = 0, (2.4) where

a = /35 + YP, b = (35yp, c = b - a + 5 + p, e = (ft + y) - b. (2.5) Hereinafter we will assume that the parameters of the model satisfy the condition

IextY + ^ = 0. (2.6)

The substitution v (t) = aot-p into the leading terms of Eq. (2.4):

vttt + v2 vtt + 2v(vt )2 (2.7)

1 Q

gives a fractional value p = ^ for the pole order of its exact solution and —1, ^ for corresponding Fuchs indices. Therefore, Eq. (2.4) fails the Painleve test and its general solution cannot be decomposed into a Laurent series. However, it is easy to show that the corresponding Puiseux series contains three arbitrary constants, that is, Eq. (2.4) passes the Painleve test in a weak form [15]. This fact gives us the possibility of finding single-valued partial solutions in a closed form.

3. Finding an exact limited solution of the single equation

After replacing the variable v (t) = \Ju (t), equation (2.4) is reduced, taking into account (2.6), to the homogeneous equation

u2uttt + u (u (u + a - 1) - ) utt + 7 (ut)3 ~

V 2 ' 4 (3.1)

12

— -u (a — u — 1) (Ut)2 + u2 {an, + c) ut H—u3 (bu + 3e) = 0, 23

the solution of which has a simple pole.

In accordance with the singular manifold method [16-18], we will seek a solution to Eq. (3.1) in the form

« = + «i (t) • (3.2)

Substituting (3.2) into (3.1) and collecting terms with the largest (in absolute value) negative degrees of F(t), we express the functions u0 and u\ in terms of F(t). Substitution (3.2) takes the form

3 3 Ft 3 Ftt

u = - H------, (3.3)

2 2 F 4: Ft

therefore, the function u(t) satisfies the Riccati equation

ut = A + Bu + Cu2, (3.4)

where

33

5 = 2, C = -2/3, A = ----(5F), (3.5)

and (SF) is the Schwarzian derivative for the function F(t):

(SF) = *j(3-6)

Note that the expression (3.6) is in the general case a function of time. Let us show that in our case it is a constant.

Let us express from (3.4) ut, utt and uttt in terms of the function u and the Schwarzian derivative:

2 3 3

ut = --u2 + 2u----(SF),

8 3 / 3 \

Utt = gti3 - 4u2 + 6m - 3 - -(SF)t + i u - -J (SF),

uttt = -yt/4 + yt/3 - 24ii2 + 24m - 9 - ^{SF)tt +

(3.7)

+ ( « - I) (SF)t - -4(SF)2 - f^u2 - 8u + 6 ) (SF).

Substituting (3.7) into (3.1) and demanding that the coefficients with the same powers of u be zero, we have

7(SF) + 2 - 8 (a + c - b) = 0, 1

(3.8)

(SF)t + -a (SF) - 3a - 4 (c + e) = 0,

12

(SF)tt + (a - 2) (SF)t + -(SF)2 + (c + 2) (SF) + 2c + 3 = 0,

[(SF) + 2] [3(SF)t + ((SF) + 2) (a - 4)] = 0,

[(SF) + 2]3 = 0.

From the last equation of system (3.8) it follows that

(SF) = -2, (3.9)

therefore, the Schwarzian derivative (3.6) does not depend on the variable t. The remaining equations of system (3.8) are satisfied if we take

a = - (c + e),

b = - e.

(3.10)

Equalities (3.10) are conditions on the coefficients of Eq. (3.1), under which it is equivalent to the homogeneous Riccati equation

o

(3-11)

2

Ut = 2 u — -it, 3

the general solution of which contains the integration constant Ci and gives an exact kink-shaped solution to Eq. (2.4):

v = v u =

3

1 + Cl e_2i

1/2

3 / , ( ln Cl

- il + tanh (t -

1/2

(3.12)

The solution (3.12) of Eq. (2.4) can be obtained in a simpler way. Noting that (2.4) does not explicitly contain an independent variable, by replacing

x = v, g (x) = vt (3.13)

the order of Eq. (2.4) can be reduced to the second one:

929xx + 9(9x)2 + (•x2 + a - l) ggx + 2xg2 + {ax2 + c) g + x (j;bx2 + e^

H =0.

(3.14)

Compensation of the leading terms of Eq. (3.14) with the substitution g = a0x-p is achieved at p = -3. This means that the solution of Eq. (3.14) can be found in the form of a 3rd-order polynomial in x:

g = b0 + b1x + b2 x2 + b3x3. (3.15)

The equations obtained after substituting (3.15) into (3.14) and grouping in powers of x are satisfied when the equalities

bo = 0, h = 1, b2 = 0, 63 =

(3.16)

are fulfilled along with conditions (3.10). Subsequent substitution of (3.16) into (3.15), then (3.15) into (3.13) gives a 1st-order equation for the function v(t):

13

Vt = V - -V"

(3.17)

the general solution of which coincides with (3.12).

4. Finding an exact solution of the initial system

Let us go back from Eq. (2.4) to the original system (1.1). Conditions (3.10), under which there is an exact solution (3.12) of Eq. (2.4), lead to equalities for the coefficients of system (1.1):

7 = —P, » = S. (4.1)

When conditions (4.1) are fulfilled, the determinant (2.3) becomes zero and Eq. (2.4) loses its connection with system (1.1) from which it was derived. Substituting conditions (4.1) directly into system (1.1) leads to its degeneration:

1 3 a

vt = v--v — w + y + — , 3 p

wt = 5v — 5Pw, (4.2)

yt = 5v — 5Py — 5a.

In fact, after substituting

y = w~jj (4-3)

in (4.2), the 3rd equation can be eliminated from the system, since it identically coincides with the 2nd one:

1 3

Vf = v--V ,

3 ' (4.4)

wt = 5v — 5pw.

The 1st of the two remaining equations of the system (4.4) does not differ from (3.17) and has the same general solution (3.12). The 2nd equation of (4.4) has a general solution

where C2 is the integration constant. The integral in (4.5) after the change of the variable 2 = exp (P5t) is expressed in terms of the generalized hypergeometric function depending on z.

Equalities (3.12), (4.5) and (4.3) give an exact solution to system (1.1) under conditions (4.1), and this fact is easily verified by direct substitution. This solution contains only two arbitrary constants Ci and C2 because between the variables w and y there is an algebraic dependence (4.3) instead of a differential one and the true order of the degenerate system (4.2) is 2. Thus, in the process of solving Eq. (2.4) using the singular manifold method, conditions (4.1) were found under which not only Eq. (2.4), but also the original system (1.1) has an exact solution. Note that the exact particular solution to Eq. (2.4) is included in the exact solution to system (1.1) as an integral part, despite the fact that the equation cannot be derived from the system when conditions (4.1) are fulfilled.

The graphs of the functions v(t) and w(t) of the exact solution (3.12), (4.5) for different values of P, 5, a and initial conditions are shown in Figs. 1a-1d.

As we see, in all cases the system asymptotically tends to a stationary state:

v = y/3, w = —, y = —(4.6)

which is a solution to system (1.1), if the time derivatives are set to zero in it. Linear stability analysis shows that state (4.6) is stable to small deviations of functions v(t), w(t) and y(t) from stationary values at P5 > 0.

(a) <7 = 2, ¿=1, /3 = 2, u(0) = 1, w(0) =0 (b) a = 2, Ö = 1, ß = 2, u(0) = 1, iü(0) = 1.5

3

2

0

0 1 2 3 4 5 6 7 t

(c) cr = 1, Ô = 1, /3 = 0.5, u(0) = 1, w(0) =0 (d) ct = 1, Ô = 1, /3 = 0.5, i;(0) = 1, w(0) = 1.5 Fig. 1. Plots of the exact solution (3.12), (4.5); v(t) — solid line, w(t) — dashed line.

5. Asymptotic solution

Equation (2.4) under condition (2.6) has a solution in the form of a simple exponent v (t) =

= Ae mt when the equalities

e = m3 + (1 — a) m2 + cm, b = 3m (a — 3m)

(5-1)

are true. Such a solution is not limited, but it has a physical meaning when considering the evolution of a system from some initial instant of time if m> 0. As t the system tends

to a stationary state

2 / Ö

v = 0, w = 0, y = —a I m + 1 - —-

3 V — m

We will try to find a solution to Eq. (2.4) in the form

v (t) = A (z) e

—mt

(5-2)

(5-3)

where A(z) and z = f (t) are unknown functions. Substituting (5.3) into (2.4), we group the result in powers of the function A(z) and its derivatives.

Equating constant coefficients at A and A3 to zero, we obtain conditions (5.1). The coefficients at Az, Azz and A2Az give an overdetermined system of linear homogeneous equations for f (t):

fttt + (a - 3m - 1) ftt + (c + 2(1 - a) m + 3m2) ft = 0,

3ftt + (a - 3m - 1) ft = 0, (5.4) ftt + (a - 6m) ft = 0. System (5.4) becomes consistent and has a solution

f (t) = Ci + C2e(6m-a)t, (5.5) where C1 and C2 are arbitrary constants, under the conditions

a = i(15m-l), c=i(33m2-12m + l). (5.6)

The remaining nonzero terms containing Azzz, A2Azz and A(Az)2 form a variable coefficient 3rd-order nonlinear equation:

o 1, . ,,, r -,

2(Az )2 + AAzz = 0. (5.7)

- 7-e"5(m+1)tA

~~~ (3m - 1) C2

In the case of m > —1, we can consider the asymptotic reduction of equation (5.7) in the limit t —

Azzz = 0. (5.8)

The general solution of Eq. (5.8) allows us, taking into account (5.3), (5.5) and (5.6), to write down the asymptotic solution of Eq. (2.4):

v (t) = 5ie"^(5m_1)i + B2e-{-im~l)t + B3e~mt, (5.9)

containing three arbitrary constants B\, B2, B3.

Expression (5.9) is bounded when t > 0, m ^ 1/4 and can be used as an analytical approximation of a numerical solution. The graphs of the asymptotic solution (5.9) for two sets of constants Bi, B2, B3 and the corresponding numerical solution of Eq. (2.4) with the same initial conditions are shown in Figs. 2a-2b. The numerical solution is obtained on the basis of the Maple computing environment using the Rosenbrock method. Note that, in the case of m ^ —1, the solution in the form (5.3) is unboundedly increasing in absolute value as t — and, therefore, has no physical meaning.

6. Conclusion

We have considered the analytical properties of the FitzHugh-Rinzel model describing the bursting activity of a neuron. We have found that the model's system of equations is reduced to a 3rd-order nonlinear differential equation, which does not pass the Painleve test. With the help of direct substitution, an exact exponential solution of the equation was found. The exact kink-shaped solution of the equation was obtained in two ways: by the singular manifold method and by the order reduction method. The conditions were found under which there are exact solutions for both the equation and the original system. If these conditions are met, then the differential order of the system is reduced by one. It is shown that the asymptotic solution of the equation in the form of a sum of three exponential functions agrees well with the numerical solution. Graphs of exact and asymptotic solutions for different sets of parameters are given.

0.3-

0.4

0.3

0.2

0.2

0.1-

0.1-

0

20

-0.1-

0

-0.2-

0 5 10 15 20

t

(a)B1 = -2,B2 = -l,B3 = 3

(b)Bi = l,B2 = 2,B3 = -2.6

Fig. 2. Plots of numeric (solid line) and asymptotic (dashed line) solutions.

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] Nagumo, J., Arimoto, S., and Yoshizawa, S., An Active Pulse Transmission Line Simulating Nerve Axon, Proc. of the IRE, 1962, vol. 50, no. 10, pp. 2061-2070.

[4] FitzHugh, R., A Kinetic Model for the Conductance Changes in Nerve Membranes, J. Cell. Comp. Physiol., 1965, vol.66, suppl. 2, pp. 111-117.

[5] Mornev, O. A., Tsyganov, I.M., Aslanidi, O.V., and Tsyganov, M. A., Beyond the Kuramoto-Zel'dovich Theory: Steadily Rotating Concave Spiral Waves and Their Relation to the Echo Phenomenon, JETP Lett., 2003, vol.77, no. 6, pp. 270-275; see also: Pis'ma v Zh. Eksper. Teoret. Fiz., 2003, vol. 77, no. 6, pp. 319-325.

[6] Llibre, J. and Vidal, C., Periodic Solutions of a Periodic FitzHugh - Nagumo System, Internat. J. Bifur. Chaos Appl. Sci. Engrg, 2015, vol.25, no. 13, 1550180, 6pp.

[7] Tsuji, Sh., Ueta, T., Kawakami, H., and Aihara, K., A Design Method of Bursting Using Two-Parameter Bifurcation Diagrams in FitzHugh - Nagumo Model, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 2004, vol. 14, no. 7, pp. 2241-2252.

[8] Gelens, L., Anderson, G. A., and Ferrell, J.E., Jr., Spatial Trigger Waves: Positive Feedback Gets You a Long Way, Mol. Biol. Cell, 2014, vol. 25, no. 22, pp. 3486-3493.

[9] Cheng, X., and Ferrell, J.E., Jr., Apoptosis Propagates through the Cytoplasm As Trigger Waves, Science, 2018, vol.361, no. 6402, pp. 607-612.

[10] Kudryashov, N. A., Asymptotic and Exact Solutions of the FitzHugh - Nagumo Model, Reggul. Chaotic Dyn, 2018, vol. 23, no. 2, pp. 152-160.

[11] Kudryashov, N.A., Rybka, R. B., and Sboev, A. G., Analytical Properties of the Perturbed FitzHugh-Nagumo Model, Appl. Math. Lett., 2018, vol.76, pp. 142-147.

[12] Rinzel, J., A Formal Classification of Bursting Mechanisms in Excitable Systems, in Mathematical Topics in Population Biology, Morphogenesis and Neurosciences, E. Teramoto, M. Yamaguti (Eds.), Lecture Notes in Biomath., vol.71, Berlin: Springer, 1987, pp. 267-281.

[13] Belykh, V. N. and Pankratova, E. V., Chaotic Synchronization in Ensembles of Coupled Neurons Modeled by the FitzHugh - Rinzel System, Radiophys. Quantum El, 2006, vol.49, no. 11, pp. 910-921; see also: Izv. Vyssh. Uchebn. Zaved. Radiofizika, 2006, vol.49, no. 11, pp. 1002-1014.

[14] Wojcik, J. and Shilnikov, A., Voltage Interval Mappings for an Elliptic Bursting Model, arXiv:1310.5232 (2013).

[15] Conte, R., Singularities of Differential Equations and Integrability, in Introduction to Methods of Complex Analysis and Geometry for Classical Mechanics and Nonlinear Waves (Chamonix, 1993), D.Benest, C.Froeschle (Eds.), Gif-sur-Yvette: Frontieres, 1994, pp. 49-143.

[16] Weiss, J., Tabor, M., and Carnevale, G., The Painleve Property for Partial Differential Equations, J. Math. Phys., 1983, vol. 24, no. 3, pp. 522-526.

[17] Weiss, J., The Painleve Property for Partial Differential Equations: 2. Backlund Transformation, Lax Pairs, and the Schwarzian Derivative, J. Math. Phys., 1983, vol. 24, no. 6, pp. 1405-1413.

[18] Kudryashov, N. A., From Singular Manifold Equations to Integrable Evolution Equations, J. Phys. A, 1994, vol. 27, no. 7, pp. 2457-2470.