Asynchronous Chaos and Bifurcations in a Model of Two Coupled Identical Hindmarsh – Rose Neurons Текст научной статьи по специальности «Математика»

Russian Journal of Nonlinear Dynamics, 2021, vol. 17, no. 3, pp. 307-320. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd210305

NONLINEAR PHYSICS AND MECHANICS

MSC 2010: 37G35,37G15,92B20

Asynchronous Chaos and Bifurcations in a Model of Two Coupled Identical Hindmarsh - Rose Neurons

I. R. Garashchuk

We study a minimal network of two coupled neurons described by the Hindmarsh-Rose model with a linear coupling. We suppose that individual neurons are identical and study whether the dynamical regimes of a single neuron would be stable synchronous regimes in the model of two coupled neurons. We find that among synchronous regimes only regular periodic spiking and quiescence are stable in a certain range of parameters, while no bursting synchronous regimes are stable. Moreover, we show that there are no stable synchronous chaotic regimes in the parameter range considered. On the other hand, we find a wide range of parameters in which a stable asynchronous chaotic regime exists. Furthermore, we identify narrow regions of bistability, when synchronous and asynchronous regimes coexist. However, the asynchronous attractor is monostable in a wide range of parameters. We demonstrate that the onset of the asynchronous chaotic attractor occurs according to the Afraimovich-Shilnikov scenario. We have observed various asynchronous firing patterns: irregular quasi-periodic and chaotic spiking, both regular and chaotic bursting regimes, in which the number of spikes per burst varied greatly depending on the control parameter.

Keywords: coupled neurons, synchronization, chaos, Hindmarsh-Rose, bursting

1. Introduction

A reliable model of a neuron should be capable of exhibiting several types of behavior observed for biological neurons, such as resting, tonic spiking and bursting [1-4]. Resting corresponds to simple silence of a neuron without activation. Mathematically such behavior could be well described by a stable equilibrium of a dynamical system. Spiking refers to continuous firing with short intervals of quiescence between the peaks, usually called interspike intervals (ISIs). Bursting refers to a regime of neural activity in which a neuron repetedly fires a discrete group

Received August 4, 2021 Accepted September 5, 2021

This work, except for Section 4, was supported by RFBR grant 20-31-90122. Section 4 was supported by RSF grant 19-71-10048

Ivan R. Garashchuk [email protected]

HSE University

34 Tallinskaya str., Moscow, 123458, Russia

of spikes (a burst), followed by a long period of quiescence before the next burst. Bursting patterns are considered to play an important role in biological processes [1, 2, 5]. Notice that both spiking and bursting can be regular, as well as irregular. A good mathematical portrait of a regular spiking or bursting pattern can be a periodic regime of a dynamical system, represented by a stable limit cycle in the phase space. Irregular firing patterns are represented by chaotic or quasi-periodic attractors.

Here we study a model of coupled neurons where each element is described by the Hind-marsh-Rose model [6]. The Hindmarsh-Rose model is capable of exhibiting all the desirable patterns of neural activity, while, on the other hand, it is one of the relatively simple models derived within the Hodkin-Huxley formalism [7], which makes it convenient to study. Due to this fact, it is one of the most popular neural models, which was thoroughly studied by various researchers, see, e.g., Refs. [1, 4, 8-15] and references therein. Another important property of the Hindmarsh-Rose model is multistability, see Refs. [16, 17] and references therein, which enriches the dynamics but complicates the analysis.

On the other hand, the topic of synchronization in the neural networks is interesting from both biological and mathematical points of view [18-21]. It is supposed that synchronized periodic firing of a group of neurons could lead to pathological behavior, such as epileptic seizures [18], while weakly chaotic almost periodic asynchronous activity is considered as a normal regime of functioning. Therefore, it is interesting to study whether the dynamical regimes of a single neuron would manifest themselves as stable synchronous regimes in a model of coupled neurons, or we would observe new asynchronous types of oscillations for certain ranges of parameters. What is more, since the three-dimensional Hindmarsh-Rose model exhibits multistability in various ranges of control parameters [13, 16, 17], it would be interesting to see if the multistability will be even more common or, conversely, less common, in a model of coupled neurons.

In this work we study a minimally possible network of two interacting Hindmarsh-Rose neurons. Here we consider the three-dimensional Hindmarsh-Rose systems describing each neuron, and the coupling force caused by the electric link between the neurons [21-23]. We also suppose that the neurons are identical, which is quite a strong condition. Nevertheless, as we show below, even under such circumstances there exists a wide range of values of the control parameter in which no stable synchronous regimes are observed.

The rest of this work is organized as follows. In the next section we introduce a model of two linearly coupled Hindmarsh-Rose neurons and discuss its properties. In Section 3 we describe and discuss bifurcations and dynamical regimes we observe in the model considered. Section 4 is devoted to chaotic and regular bursting solutions in the model. In the last section we briefly discuss and summarize our results.

2. The main system of equation

The model of two identical electrically coupled Hindmarsh -Rose neurons is represented by the following system of equations:

x 1 = y1 — axf + bx\ — z1 + I + D0(x2 — x1),

ii = c — dx\ — yi,

Z1 = r(s(x1 — x0) — z1),

x 2 = i2 — axl + bxl — Z2 + I + Do(xi — x2),

i)2 = c — dx2 — i2,

Z2 = r(s(x2 — x0) — Z2),

2

(2.1)

where xi, i = 1, 2 are the potentials on the membranes of the respective neurons, while yi and zi, i = 1, 2 are the fast and slow ionic currents, respectively. We restrict the system to be composed of two identical elements, therefore we use the same parameter values for both neurons. The system is a fast-slow one, because the time scales of zi variables are defined by the small parameter r. Here we choose the value r = 0.0021, following Ref. [13]. The other parameters, except for external synaptic current I, have the following fixed values: a = 1, b = 3, c = 1, d = 5, s = 4 and x0 = —8/5. The coupling strength between the neurons is chosen to be D0 = 0.1, which corresponds to a relatively strong connection between the pair of neurons, but within a biologically reasonable range [19-21, 23]. The parameter I plays the role of the bifurcation parameter in our work.

For numerical integration of the model (2.1) we use a quasi-constant step size method of the Klopfenstein-Shampine family [24, 25], which is suitable for stiff systems. To calculate the spectrum of Lyapunov exponents we use the standard algorithm by Bennetin et al. [26]. Since the transient processes for the system under study can be quite long (of the order of t ~ we start calculating the Lyapunov exponents only after skipping it.

Because we assume the neurons to be identical, the system (2.1) possesses the symmetry of swapping the neurons, namely, x1 o x2, y1 o y2, z1 o z2. This fact leads to the existence of the synchronization manifold x1 = x2, y1 = y2, z1 = z2 in the phase space. Trajectories embedded in the synchronization manifold correspond to completely identical time series for all the phase variables describing both neurons: Vt, x1(t) = x2(t), y1(t) = y2(t), z1 (t) = z2(t). This manifold is a three-dimensional plane in the six-dimensional phase space. One can notice that inside the synchronization manifold the terms related to coupling of the elements vanish in the respective equations for obi. Thus, all the invariant sets found in the three-dimensional model of a single Hindmarsh-Rose neuron (e.g., in Refs. [8, 9, 13]) exist in the synchronization manifold of the six-dimensional model (2.1). However, they can become unstable in directions transversal to the synchronization manifold and nonattractive (for example, saddle) in the model of two coupled neurons.

3. Dynamical regimes

We study the influence of the bifurcation parameter I on the behavior of the model (2.1), following Ref. [13]. For I > 26 there exists a single stable fixed point, and it is a global attractor. Therefore, we will focus on the interval 0 < I < 26. In this interval there exists only one fixed point (referred to below as FP), which lies inside the synchronization manifold and corresponds to the FP found in Ref. [13]. No new asynchronous FPs were found.

Despite three new dimensions transversal to the synchronization manifold appearing in the six-dimensional system, the stability regions of the FP do not change. It is stable in the following domain: I < 1.2895, 5.3978 < I < 6.1976 and I > 25.261. We start from the right end of the range of I and study changes occurring in the dynamics as I decreases. The FP loses its stability at I = 25.261 via a supercritical Andronov-Hopf bifurcation, giving rise to a stable synchronous limit cycle, as in the three-dimensional case. The limit cycle at I = 25 is presented in Fig. 1. However, unlike in the three-dimensional system, the aforementioned cycle only remains stable up until I = 7.394, while in the three-dimensional case it is stable all the way up to I = 3.220. At I = 7.394 it undergoes a Neimark-Sacker bifurcation, becomes unstable and a quasi-periodic attractor is born.

Fig. 1. Phase portrait of the stable limit cycle (black) and the unstable FP (red) for I = 25.

4.8 4.85 4.9 4.95^.5 5.05 5.1

Fig. 2. Three largest Lyapunov exponents of the continuation of the synchronous limit cycle emerging at large values of I in the range 0 < I < 9.5. We follow the asynchronous attractor after the Neimark - Sacker bifurcation at I = 7.394. Graphs (b) and (c) show enlarged parts around the Neimark-Sacker bifurcation and the transition to asycnhronous chaos, respectively.

The graph of the largest Lyapunov exponents for the continuation of the aforementioned limit cycle in the interval 0 < I < 9.5 is presented in Fig. 2a. The part where the Neimark -Sacker bifurcation occurs is shown in more detail in Fig. 2b. The projections of the phase portraits and the Poincare sections of the stable limit cycle at I = 7.410, before the bifurcation point, and the quasi-periodic attractor at I = 7.388, which emerged via the Neimark-Sacker bifurcation, are presented in Fig. 3 and Fig. 4, respectively. The Poincare section for the quasi-periodic attractor is formed by the intersection with the plane z2 = 5.976 in the direction of increasing z2. It can be clearly seen that, while the periodic attractor lies inside the diagonal synchronization plane, the quasi-periodic one spreads outside of it. The firing pattern of both

0.00

x1

2.00

-10.00

2/i

Fig. 3. Phase portrait of the stable limit cycle for I = 7.410.

H 2.05-

2.10-

(b)

IQP]

Fig. 4. The phase portrait (a) and the Poincare section (b) of the quasi-periodic attractor for I = 7.388.

neurons corresponding to the synchronous limit cycle is spiking (since it is similar to those in Ref. [13]). For the quasi-periodic regime, the oscillations are not periodic any more, but the firing pattern generally looks similar to the spiking exhibited by the periodic regime, i.e., there is no bursting of several oscillations of fast variables per peak or similar effects.

The quasi-periodic regime remains monostable up to I = 6.1976, when the FP becomes stable again via a subcritical Andronov- Hopf bifurcation, and system (2.1) becomes bistable. The FP is stable for 5.3978 < I < 6.1976. At I = 5.3978 a supercritical Andronov - Hopf bifurcation takes place, the FP becomes unstable and a stable limit cycle arises. A bifurcation diagram for the parameter values around the bistability region is shown in Fig. 5. The quasi-periodic regime remains stable and changes only slightly within this interval. Although the Andronov - Hopf bifurcation takes place at the same I-value as in the three-dimensional model, contrary to the three-dimensional system, the stable synchronous limit cycle exists in a very narrow range of the parameter, coexisting with the asynchronous quasi-periodic attractor. At I = = 5.37596 the synchronous limit cycle undergoes a saddle-node bifurcation (with an unstable limit cycle) and disappears. The phase portrait of this stable synchronous limit cycle is presented in Fig. 6. The cycle is quite small in the phase space, and the firing pattern of the fast variables is spiking. Notice also that the unstable FP belongs to the basin of attraction of the stable limit cycle in the entire range of its existence. It makes the limit cycle self-excited with respect to this FP, according to the classification in Ref. [27]. After the point of the saddle-node bifurcation, I < 5.37596, we could not find any stable limit cycles, and thus, we consider the system to be monostable until the FP becomes stable again. Therefore, the asynchronous attractor becomes monostable and self-excited with respect to the FP, while the FP remains unstable.

The FP becomes stable at I = 1.2895 via a subcritical Andronov - Hopf bifurcation. Let us first discuss the changes occurring to the asynchronous attractor before this point. It remains

(a)2 1.5

Ä1

3

s

0.5

SN

u

-QP UFP

— SFP

— LC

SH

H

V

5.4 5.5 5.6 5.7 5.8 5.9 6 6.1 6.2 I

Fig. 5. (a) Bifurcation diagram for the bistability region. With decreasing I the unstable fixed point (UFP) becomes stable via a subcritical Andronov-Hopf bifurcation (marked as SH). After the region of stability it undergoes a supercritical Andronov - Hopf bifrucation (marked as H), and a stable synchronous limit cycle (LC) emerges. With decreasing I, the synchronous limit cycle disappears via the saddle-node bifurcation (SN) with an unstable limit cycle. The system becomes monostable again after that. The quasi-periodic regime (QP) is stable in the entire range shown here. (b) shows the enlarged part around the saddle-node bifurcation of the synchronous limit cycle.

Fig. 6. Phase portrait of the stable limit cycle (black) and the unstable FP (red) for I = 5.390 after the Andronov - Hopf bifurcation.

quasi-periodic up to I & 4.9, where the asynchronous chaotic attractor emerges. The onset of the asynchronous chaos occurs according to the Afraimovich-Shilnikov scenario. The torus loses smoothness and, through the destruction of the torus, the chaotic attractor emerges (for details, see Ref. [28]). The emergence of chaotic dynamics can be seen in the graph of the largest Lyapunov exponents in Figs. 2a, 2c. The Poincare section of the quasi-periodic attractor before the loss of smoothness at I = 4.981 is presented in Fig. 7. The shape of the closed curve in the Poincare map has become more complicated compared to that shown in Fig. 3. The Poincare section of the asynchronous chaotic attractor with a small positive largest Lyapunov exponent (A1 = 0.900 • 10"3) right after the torus loses smoothness at I = 4.906 is presented in Fig. 8. Let us note that, because the position of the attractor is changing in the phase space with variation of the parameter, we have to shift the plane used to construct the Poincare section and to preserve visual consistency of the graphs. The Poincare section of a "more developed" chaotic attractor further after the bifurcation point at I = 4.786 with X1 = 0.94 • 10_2 is presented in Fig. 9.

0.0

2.0 -6.0

Vi

Fig. 7. Projection of the Poincare section of the quasi-periodic regime for I = 4.981 on the (x±, y1, x2) space. The Poincare section is formed by the intersection with the z2 = 4.87 plane in the direction of increasing z2.

Fig. 8. Projection of the Poincare section of the chaotic regime right after the loss of smoothness at I = = 4.906 on the (x1, y1, x2) space. The Poincare section is formed by the intersection with the z2 = 4.834 plane in the direction of increasing z2.

Fig. 9. Projection of the Poincare section of the developed chaos for I = 4.786 on the (x1, y1, x2) space. The Poincare section is formed by the intersection with the z2 = 4.765 plane in the direction of increasing z2.

The dynamics remains chaotic with windows of periodicity (since the attractor is not pseudo-hyperbolic) all the way to I = 1.2895. Notice also that, while for I > 1.2895 the system is monostable, for parameter values close to the bifurcation point, if the initial conditions are chosen in the vicinity of the unstable FP, the transient process before convergence to the chaotic attractor can be very long (up to t ~ 106), and the system can exhibit behavior resembling almost periodic self-oscillations for quite a significant time length. After the stable point goes through a subcritical Andronov - Hopf bifurcation, we observe one more narrow bistability region: while the chaotic attractor still exists, the FP is already stable. The graph of the largest Lyapunov exponents of the continuation of the asynchronous attractor and the real parts of the largest

eigenvalues of the Jacobian matrix at the FP around the bistability region are shown in Fig. 10. One can see that the bistability takes place in 1.2760 < I < 1.2895. At I — 1.2760 the chaotic asynchronous attractor faces a crisis, and the stable FP becomes a global attractor. The system is monostable for I < 1.2760.

xlO"3

1 1 1 Attractor i - crisis ' | \ '1 \ '1 \ 1 -0 -X1 -A2 — ■e1 2 — -634

\

i i " i i i i i i ~ i i i if/ Hop! 1 1 F J-

1.27 1.275 1.28 1.285 1.29 1.295 1.3

I

Fig. 10. Chart of Lyapunov exponents showing the bistability region. Here A1 and A2 correspond to the largest Lyapunov exponents of the asynchronous attractor, while e12, e34 represent the real parts of the complex-conjugated pairs of eigenvalues of the Jacobian matrix at the FP. If I decreases, the FP becomes stable via a subcritical Andronov-Hopf bifurcation (marked as a Andronov-Hopf point), where the pair ei 2 crosses the imaginary axis. The asynchronous attractor undergoes a crisis at I — 1.2760, after which the largest Lyapunov exponents coincide with the eigenvalues of the FP. The bistability takes place between the Andronov - Hopf point and the point of crisis of the asynchronous attractor.

4. Chaotic and regular bursting solutions

In this section we discuss chaotic and regular bursting solutions in system (2.1). If we start from higher values of the parameter I, then the stable synchronous limit cycle represents a spiking pattern equivalent to the one in the single neuron case. After the Neimark-Sacker bifurcation, the activation pattern of the quasi-periodic regimes is still spiking, however, interspikes intervals become aperiodic and the spikes heights can vary (see an example for I — 4.981 in Fig. 11). In addition, the oscillations are not synchronous any more, although much of the time they are close to synchronous.

The chaotic attractor emerging from the quasi-periodic one through the Afraimoivich -Shilnikov scenario at first exhibits chaotic spiking behavior — the ISIs and the spikes heights become irregular, and anomalously long quiescence intervals become possible, which are not observed in the quasi-periodic regime. Moreover, desynchronization becomes much more apparent than in the quasi-periodic case (see an example for I — 4.786 in Fig. 12). It is also interesting that the anomalously long resting intervals for individual neurons might not align given that the oscillations are asynchronous.

H

2 1 0 -1.

0 50 100 150 200 250 t

Fig. 11. Quasi-periodic spiking at I = 4.981 (corresponds to the phase portrait in Fig. 7).

-1

80

159

239

4150

4200

t

4250

4300

Fig. 12. Chaotic spiking and anomalous resting periods at I = 4.786 (corresponds to the phase portrait in Fig. 9) at different time intervals. The resting period of the second neuron is significantly longer in (a), while that of the first one, in (b).

As I decreases further, the frequency of appearance of such anomalous intervals of quiescence increases. From I = 4.786 to I = 4.3080 the frequency of occurrence of such events increases about tenfold. Also, the length of the largest quiescence periods increases approximately twofold. An example of two closely lying, but relatively short anomalous quiescence periods at I = 4.3080 are presented in Fig. 13a, while a deep and long resting interval is shown in Fig. 13b. The long periods of quiescence tend to align much better than their narrower and less deep counterparts, which tend to occur quite randomly for each individual neuron and may not align at all, as shown in Fig. 13a. It is also interesting that in Fig. 13b the spiking patterns to the left and to the right of the long resting period differ significantly (ISIs and heights of the spikes).

Fig. 13. Chaotic spiking regime at I = 4.3080. (a) Shorter anomalously long ISIs and (b) a long anomalous quiescence period

250 500 " 750 0 150 300 . 450

Fig. 14. (a) Intermediate regime between spiking and bursting at I = 3.7880 and (b) bursting (with long bursts) at I = 3.4680

If we decrease I further, the frequency of occurrence of long resting periods interrupting the spiking keeps increasing. Eventually, at some point the firing pattern starts resembling something in between spiking, often interrupted by the long quiescence periods, and extremely long bursts (see, e.g., the typical "bursts" at I = 3.7880 in Fig. 14a). The first group of spikes in Fig. 14a contains 24 peaks, and the second one is even larger. At lower values of I the firing pattern finally starts looking closer to what is usually considered as long bursts (see, e.g., Fig. 14b showing the activation pattern of the voltage variables at I = 3.4680).

A burst is called anomalous if the ISIs or heights of the spikes within the burst are changing nonmonotonously, and regular otherwise. Examples of regular bursts are given in Figs. 15a, 15b in which regimes from windows of periodicity at I = 3.188 and I = 2.4280 are presented. A clear difference between a regular and an anomalous burst can be seen if one compares Figs. 14a, 14b with Fig. 15a. A full period of the regime in Fig. 15a consists of two bursts, one of which has 16 spikes, while the other has 17. Although we do not strictly define the phase of oscillations in this paper, one can at least visually observe that, except for the quiescence periods, the voltages of both neurons oscillate almost in antiphases: the burst of 16 peaks of the first neuron occurs simultaneously with the 17-peaks burst of the second neuron and vice versa. Moreover, within each of the bursts, the spikes of the first neuron coincide with the interburst ISIs of the second neuron and vice versa. The same effect can be observed in Fig. 15b, but now the bursts consist of 11 and 10 spikes, respectively. The projection of the phase portrait of the regime at I = 2.4280 on (xl, y1, x2) (which corresponds to the bursting pattern in Fig. 15b) is shown in Fig. 16. The long trail emanating far from the main area, in which the largest part of the attractor (limit cycle) is lying, corresponds to the long quiescence period between the bursts. It passes close to the synchronization manifold, and it explains why the behavior of the variables x1 and x2 resembles synchronous behavior as they pass the resting period in the graphs in Figs. 14 and 15. The same effect takes place for the chaotic attractors, but the line density in the phase space is too high for chaotic regimes to visually analyze the phase portraits; therefore, we present only an example for a periodic regime.

The general tendency is the following: as I decreases from the onset of chaotic bursting towards the point of the subcritical Andronov-Hopf bifurcation and the attractor crisis, the bursting pattern of the regular regimes in the windows of periodicity and the chaotic bursting in the regions of chaos gradually simplify. For the regular regimes, the number of spikes per burst and the full period of the limit cycles decreases. For chaotic regimes, of course, the number of spikes per burst is irregular, but the average number of spikes per burst declines with decreasing I. Examples of chaotic bursting regimes with less spikes per burst at I = 2.3680 and I = 1.7080

0 150 300 t 450 0 150 300 f 450

Fig. 15. Regular bursting in windows of periodicity at (a) I = 3.1880 and (b) I = 2.4280

Vi xi

Fig. 16. Phase portrait of the regular bursting pattern in the window of periodicity at I = 2.4280. The diagonal plane represents the projection of the synchronization manifold onto the (x1, y1, x2) space. The black line and the red dot correspond to the limit cycle and the unstable fixed point, respectively. The seeming intersection of the limit cycle and the synchronization manifold does not actually happen — it is the effect of projection of the six-dimensional phase space onto the three-dimensional subspace, in which we plot the graph.

Fig. 17. Chaotic bursting at (a) I = 2.3680 and (b) I = 1.7080 RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2021, 17(3), 307-320

Fig. 18. (a) Regular bursting at I = 1.2800 and (b) chaotic bursting at I = 1.2780

are presented in Figs. 17a, 17b. One can clearly notice that the difference between such regimes and the one in Fig. 14b is quite drastic, and the average number of spikes per burst differs significantly.

Finally, at the region of bistability close to the crisis of the asynchronous attractor, we can observe bursting patterns with very low numbers of spikes per burst, such as the limit cycle at I = 1.2800, the full period of which consists of two bursts with 2 spikes each (see Fig. 18a). In the regions of chaos within the same parameter area, the average number of spikes per burst is comparable. For example, in Fig. 18b chaotic bursting at I = 1.2780, exhibiting mostly 2-4 spikes per burst is demonstrated.

5. Conclusion

In this work we have studied a pair of coupled neurons described by the Hindmarsh-Rose model with a linear coupling. We have shown that the bifurcation structure of this model is significantly different from that of the three-dimensional model of a single neuron with the same parameters. While in some regions of the bifurcation parameter the dynamics of the system is asymptotically the same as in the three-dimensional case, there exists a very wide range of values of the control parameter in which the dynamics is qualitatively different from that in the three-dimensional case.

We have found a substantially new asynchronous attractor that does not have an analogue in the three-dimensional model. At first it emerges from a synchronous stable limit cycle via the Neimark-Sacker bifurcation, giving rise to a quasi-periodic attractor. This dynamical regime corresponds to the irregular quasi-periodic spiking pattern of neural activity. Then it becomes chaotic according to the Afraimovich-Shilnikov scenario, which is a common universal scenario of the onset of asynchronous chaos in models of coupled nonlinear oscillators [29, 30]. On the other hand, the invariant sets found in the model of a single neuron become transversally unstable and nonattractive in a wide range of the bifurcation parameter, while new asynchronous attractors are monostable. In the three-dimensional system the onset of chaos occurs via a period-doubling cascade of a limit cycle, which does not happen in the six-dimensional case. Moreover, in the studied range of parameters we did not observe any stable synchronous chaotic regimes at all. The asynchronous chaotic regime persists in a significantly wider range of the parameter than the chaotic one in the single neuron model. Discarding narrow windows of periodic oscillations, in the coupled neurons model it exists approximately in the range of 1.29 < I < 4.9, while the single neuron exhibited chaotic behavior only in 3.12 < I < 3.36 [13].

Moreover, we have not found any stable synchronous bursting solutions for the studied values of the parameter, i. e., all bursting solutions existing for a single neuron and lying in the synchronization manifold become unstable. However, we have observed a variety of bursting solutions corresponding to asynchronous regimes in the chaotic region and regions of periodic windows. Some of these regimes are qualitatively similar to those found in Ref. [13], while others are not. The lengths and numbers of spikes per burst vary over a great range from minimal two spikes per burst (see Fig. 18) to extremely long "bursts", which could be borderline called spiking with sudden long resting intervals between them (see Figs. 13 and 14). Furthermore, a significantly wider range of chaotic dynamics, compared to the single neuron model, implies that in the system of two neurons the bursting solutions exist in a much wider range of the control parameter. Among the asynchronous regimes we have also found irregular quasi-periodic and chaotic spiking, while no regular spiking patterns have been observed.

In addition, the six-dimensional system is bistable in a narrower range of the parameter than the three-dimensional one, which may be beneficial from the point of view of biological applications. On the other hand, many firing patterns found for a single neuron cannot be observed in the model of coupled oscillators, because the corresponding synchronous regimes become unstable, while qualitatively different asynchronous activity patterns are observed instead. However, it might be beneficial for the biological applications if we consider that regular synchronous firing can lead to pathological behavior, such as epileptic seizures, while weakly chaotic asynchronous regimes are desirable for normal functioning of the neural networks [18].

The author thanks Dmitry Sinelshchikov and the reviewer for their valuable comments and suggestions that helped to improve the manuscript.

References

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

[2] Izhikevich, E. M., Desai, N. S., Walcott, E. C., and Hoppensteadt, F. C., Bursts As a Unit of Neural Information: Selective Communication via Resonance, Trends Neurosci., 2003, vol. 26, no. 3, pp. 161— 167.

[3] Izhikevich, E. M., Neural Excitability, Spiking and Bursting, Internat. J. Bifur. Chaos Appl. Sci. Engrg, 2000, vol. 10, no. 6, pp. 1171-1266.

[4] Bursting: The Genesis of Rhythm in the Nervous System, S.Coombes, P. C.Bressloff (Eds.), Singapore: World Sci., 2005.

[5] Lisman, J. E., Bursts As a Unit of Neural Information: Making Unreliable Synapses Reliable, Trends Neurosci, 1997, vol.20, no. 1, pp. 38-43.

[6] Hindmarsh, J. L. and Rose, R. M., A Model of Neuronal Bursting Using Three Coupled First Order Differential Equations, Proc. R. Soc. Lond. Ser. B Biol. Sci., 1984, vol. 221, no. 1222, pp. 87-102.

[7] 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.

[8] Shilnikov, A. and Kolomiets, M., Methods of the Qualitative Theory for the Hindmarsh - Rose Model: A Case Study. A Tutorial, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 2008, vol. 18, no. 8, pp.2141-2168.

[9] Barrio, R., Angeles Martinez, M., Serrano, S., and Shilnikov, A., Macro- and Micro-Chaotic Structures in the Hindmarsh - Rose Model of Bursting Neurons, Chaos, 2014, vol. 24, no. 2, 023128, 11 p.

[10] Rosenblum, M.G. and Pikovsky, A. S., Controlling Synchronization in an Ensemble of Globally Coupled Oscillators, Phys. Rev. Lett., 2004, vol.92, no. 11, 114102, 4p.

[11] Pikovsky, A., Rosenblum, M., and Kurths, J., Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge Nonlinear Sci. Ser., vol. 12, New York: Cambridge Univ. Press, 2001.

[12] Izhikevich, E. M., Which Model to Use for Cortical Spiking Neurons?, IEEE Trans. Neural Netw., 2004, vol. 15, no. 5, pp. 1063-1070.

[13] Innocenti, G., Morelli, A., Genesio, R., and Torcini, A., Dynamical Phases of the Hindmarsh - Rose Neuronal Model: Studies of the Transition from Bursting to Spiking Chaos, Chaos, 2007, vol. 17, no. 4, 043128, 11 p.

[14] Huerta, R., Rabinovich, M. I., Abarbanel, H. D. I., and Bazhenov, M., Spike-Train Bifurcation Scaling in Two Coupled Chaotic Neurons, Phys. Rev. E, 1997, vol.55, no.5, R2108-R2110.

[15] Holden, A.V. and Fan, Y.-Sh., From Simple to Simple Bursting Oscillatory Behaviour via Chaos in the Rose - Hindmarsh Model for Neuronal Activity, Chaos Solitons Fractals, 1992, vol.2, no.3, pp.221-236.

[16] Erichsen, R. and Brunnet, L.G., Multistability in Networks of Hindmarsh - Rose Neurons, Phys. Rev. E, 2008, vol.78, no. 6, 061917, 6 p.

[17] Malashchenko, T., Shilnikov, A., and Cymbalyuk, G., Six Types of Multistability in a Neuronal Model Based on Slow Calcium Current, PLoS One, 2011, vol. 6, no. 7, e21782, 10 p.

[18] Sacramento, J.F., Chew, D. J., Melo, B.F., Donega, M., Dopson, W., Guarino, M.P., Robinson, A., Prieto-Lloret, J., Patel, S., Holinski, B. J., Ramnarain, N., Pikov, V., Famm, K., and Conde, S. V., Bioelectronic Modulation of Carotid Sinus Nerve Activity in the Rat: A Potential Therapeutic Approach for Type 2 Diabetes, Diabetologia, 2018, vol.61, pp. 700-710.

[19] Castanedo-Guerra, I. T., Steur, E., and Nijmeijer, H., Synchronization of Coupled Hindmarsh - Rose Neurons: Effects of an Exogenous Parameter, IFAC-PapersOnLine, 2016, vol.49, no. 14, pp. 84-89.

[20] Yu, H. and Peng, J., Chaotic Synchronization and Control in Nonlinear-Coupled Hindmarsh - Rose Neural Systems, Chaos Solitons Fractals, 2006, vol. 29, no. 2, pp. 342-348.

[21] Eteme, A. S., Tabi, C.B., and Mohamadou, A., Synchronized Nonlinear Patterns in Electrically Coupled Hindmarsh - Rose Neural Networks with Long-Range Diffusive Interactions, Chaos Solitons Fractals, 2017, vol. 104, pp. 813-826.

[22] Zoli, M., Torri, C., Ferrari, R., Jansson, A., Zini, I., Fuxe, K., and Agnati, L.F., The Emergence of the Volume Transmission Concept, Brain Res. Rev., 1998, vol. 26, nos. 2-3, pp. 136-147.

[23] Eteme, A. S., Tabi, C. B., and Mohamadou, A., Long-Range Patterns in Hindmarsh - Rose Networks, Commun. Nonlinear Sci. Numer. Simul., 2017, vol.43, pp. 211-219.

[24] Klopfenstein, R. W., Numerical Differentiation Formulas for Stiff Systems of Ordinary Differential Equations, RCA Rev., 1971, vol.32, pp.447-462.

[25] Shampine, L.F. and Reichelt, M.W., The MATLAB ODE Suite: Dedicated to C.William Gear on the Occasion of His 60th Birthday, SIAM J. Sci. Comput., 1997, vol. 18, no. 1, pp. 1-22.

[26] Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.-M., Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems: A Method for Computing All of Them: P. 1: Theory, Meccanica, 1980, vol. 15, no. 1, pp. 9-20.

[27] Dudkowski, D., Jafari, S., Kapitaniak, T., Kuznetsov, N. V., Leonov, G. A., and Prasad, A., Hidden Attractors in Dynamical Systems, Phys. Rep., 2016, vol. 637, pp. 1-50.

[28] Afraimovich, V. S. and Shilnikov, L.P., Invariant Two-Dimensional Tori, Their Breakdown and Stochasticity, Amer. Math. Soc. Transl. (2), 1991, vol.149, pp. 201-212; see also: Methods of Qualitative Theory of Differential Equations, E. A. Leontovich-Andronova (Ed.), Gorky: Gorky Gos. Univ., 1983, pp. 3-26.

[29] Garashchuk, I. R., Sinelshchikov, D. I., Kazakov, A. O., and Kudryashov, N. A., Hyperchaos and Multistability in the Model of Two Interacting Microbubble Contrast Agents, Chaos, 2019, vol. 29, no. 6, 063131, 16p.

[30] Stankevich, N. V., Dvorak, A., Astakhov, V., Jaros, P., Kapitaniak, M., Perlikowski, P., and Kapitaniak, T., Chaos and Hyperchaos in Coupled Antiphase Driven Toda Oscillators, Regul. Chaotic Dyn., 2018, vol.23, no. 1, pp. 120-126.