Coexisting Attractors and Multistate Noise-Induced Intermittency in a Cycle Ring of Rulkov Neurons

Abstract

We study dynamics of a unidirectional ring of three Rulkov neurons coupled by chemical synapses. We consider both deterministic and stochastic models. In the deterministic case, the neural dynamics transforms from a stable equilibrium into complex oscillatory regimes (periodic or chaotic) when the coupling strength is increased. The coexistence of complete synchronization, phase synchronization, and partial synchronization is observed. In the partial synchronization state either two neurons are synchronized and the third is in antiphase, or more complex combinations of synchronous and asynchronous interaction occur. In the stochastic model, we observe noise-induced destruction of complete synchronization leading to multistate intermittency between synchronous and asynchronous modes. We show that even small noise can transform the system from the regime of regular complete synchronization into the regime of asynchronous chaotic oscillations.

1. Introduction

Coexistence of various attractors in phase space (or multistability) is a universal phenomenon in dynamical systems of various nature, including electronics, optics, mechanics, and biology (see [1] and references therein). The choice of attractors in multistable systems is determined by the initial conditions. At present, the phenomenon of multistability is widely studied in many papers devoted to the destruction of synchronous modes using bifurcation analysis, in discrete maps, genetic elements, laser systems, and ensembles of coupled oscillators (see, for example, [2,3,4,5,6]). However, a detailed analysis of synchronous dynamics of interacting units near the boundaries of various types of synchronization in multistable systems was not yet carried out.

In recent years, the interest of many researchers in a study of multistable systems synchronization [7] has shifted from two coupled periodic oscillators to complex systems of interconnected chaotic units (see [8] and references therein). Among extensive research devoted to synchronization of different systems, a study of synchronous dynamics of coupled neurons takes an important place, since the simulation of such systems helps to better understand brain activity mechanisms, as well as to reveal general concepts of key dynamical regimes in coupled systems of different nature. Many applications use identical dissipative neural generators. To simulate cooperative neuron dynamics, numerous models were developed based either on iterative maps [9,10,11,12] or differential equations [13,14,15,16] in various communication configurations.

Interest in the study of collective dynamics of ring-coupled oscillators grew significantly after publication of the Turing’s pioneering paper on morphogenesis [17]. Later, ring geometry was widely explored in a number of physiological and biochemical applications [18,19] because ring-shaped network motifs often occur in biological systems, for example, in peripheral nervous systems and locomotion [20,21,22]. A motif is a link pattern that appears much more often than in randomized network versions. The network motif concept was first used to characterize patterns of relationships in the regulation network of Escherichia coli genes [23].

At the same time, all natural systems, including neurons and the brain, are known to be subject to random fluctuations. Therefore, the study of stochastic phenomena in interacting non-linear systems and development of new methods for their analysis is an attractive problem for modern mathematical modeling [24]. Among them, it is worth mentioning noise-induced transitions [6,25,26], stochastic and coherent resonance [27,28], stochastic generation of patterns [29], noise-induced bifurcations [30] and chaos [31,32,33], and stochastic excitability [34]. In the studies of such phenomena, along with direct numerical simulations [35,36], asymptotics and approximations were actively used [37,38].

In this paper, we are interested in understanding how synchronization occurs in a unidirectional ring of three identical Rulkov neurons coupled by chemical synapses. Neurons in a neural network are known to be connected by either electrical or chemical synapses. While electrical synapses are generally bidirectional, chemical synapses provide one-way communication. Neurons that are initially in a stable equilibrium when they are uncoupled follow a path to chaos as the coupling strength is increased. On this route, a number of attractors are born that coexist in a certain range of the coupling strength and other parameters. When studying the dynamics of a ring of coupled oscillators, some researchers paid special attention to the dynamics of rotational waves circulating along the ring [39,40,41,42,43]. Although in this work we also observe such waves, we mainly focus on the mechanisms for the emergence of multistability and noise-induced intermittency. We will also show how stochastic perturbations induce switching between coexisting states, thus resulting in multistate intermittency.

2. Deterministic System

Consider a ring of three identical neurons coupled unidirectionally as follows

$$\begin{array}{c}{x}_{t+1}=f(\gamma ,x_t)-\sigma \phi (x_t,z_t),\hfill \\ y_{t+1}=f(\gamma ,y_t)-\sigma \phi (y_t,x_t),\hfill \\ z_{t+1}=f(\gamma ,z_t)-\sigma \phi (z_t,y_t),\hfill \end{array}$$

(1)

where $\sigma$ is the coupling strength used as a control parameter. In the model (1), the neurons are modeled by the Rulkov function [44]

$$f(\gamma ,x)=\frac{\alpha }{1+x^2}+\gamma ,$$

where $\alpha =4.1$ and $\gamma$ is a coefficient which regulates the neuron dynamics. Since we are interested in a one-way connection, the coupling is realized through a chemical synapse, given as

$$\phi (x,z)=\frac{x-v}{1+exp[-k(z-\theta \left)\right]},$$

where $v=-1.2$, $\theta =-1.55$, and $k=50$.

Figure 1a shows the bifurcation diagram of the isolated Rulkov map $x_{t+1}=f(\gamma ,x_t)$ using $\gamma$ as a control parameter. As seen from the diagram, this model demonstrates a regular behavior for $\gamma >-1.25$ and chaos for $\gamma <-1.25$.

In this paper, we consider the case when isolated neurons are in stable equilibrium modes. In the system $x_{t+1}=f(\gamma ,x_t)$, the stable equilibria are observed for $\gamma >{\gamma }^{*}=0.50795$. As the parameter $\gamma$ passes ${\gamma }^{*}$ from right to left, the equilibrium loses its stability, and the stable 2-cycle appears. Therefore, in what follows we fix $\gamma =0.6$ and study dynamics of the system (1) depending on the coupling parameter $\sigma$.

To analyze the stability of the equilibrium $\overline{x}=\overline{y}=\overline{z}$ in the system (1), we consider the corresponding Jacobi matrix

$$J=\left[\begin{array}{ccc}a-\sigma b& 0& -\sigma c\\ -\sigma c& a-\sigma b& 0\\ 0& -\sigma c& a-\sigma b\end{array}\right],$$

where

$$a=-\frac{2\alpha \overline{x}}{{(1+{\overline{x}}^2)}^2},b=\frac{1}{1+exp[-k(\overline{z}-\theta )]},c=kexp[-k(\overline{z}-\theta )]\frac{x-v}{{\left(1+exp[-k(\overline{z}-\theta )]\right)}^2}.$$

Eigenvalues of the Jacobi matrix can be found explicitly:

$${\lambda }_{1,2,3}=a-\sigma (b+c).$$

The general condition $|{\lambda }_{1,2,3}|<1$ of asymptotic stability of the equilibrium for the considered set of parameters $\alpha =4.1,\gamma =0.6,v=-1.2,\theta =-1.55$, and $k=50$ is equivalent to the inequality $\sigma <{\sigma }^{*}=0.020154$. As the parameter $\sigma$ passes ${\sigma }^{*}$ from left to right, the equilibrium loses its stability, and the stable 2-cycle appears.

Changes in the system (1) behavior with increasing $\sigma$ are illustrated in Figure 1b, where we present the bifurcation diagram of the x-coordinate of the system (1) on the $x=y=z$ line versus $\sigma$. The red line in the figure shows the largest Lyapunov exponent, which is a standard criterion to identify regular ($\mathrm{\Lambda }<0$) and chaotic ($\mathrm{\Lambda }>0$) dynamics. One can see that for weak coupling ($\sigma <0.48$) the system is in a regular mode, whereas for stronger coupling the ring system turns to chaos; a 2-cycle window is observed for $0.74<\sigma <0.84$. As can be seen, even in this restricted subspace of ${\mathbf{R}}^3$, the system exhibits the diversity of regular and chaotic oscillatory regimes.

Next, we illustrate the coexistence of different dynamical regimes in the system (1) for various values of $\sigma$. Note that the system (1) is monostable on the $x=y=z$ line and multistable in ${\mathbf{R}}^3$. In particular, for $\sigma =0.2$ the system (1) exhibits the coexistence of four 2-cycles shown in Figure 2a by color dots. For ease of visualization, the dots are connected by color lines. The cycle in red exhibits complete synchronization ($x=y=z$) (see time series in Figure 2b), while other three cycles exhibit partial synchronization. Namely, the green line in Figure 2a represents a cycle in which x-neuron and y-neuron are identical while z-neuron is in antiphase (see time series in Figure 2c). Here, the antiphase dynamics are manifested in the fact that both x and y values alternate with the z value, i.e., when x and y neurons fire, z-neuron is at rest state and vice versa. The other coexisting states are shown by the blue and pink lines. They indicate, respectively, a cycle in which x-neuron and z-neuron are identical while y-neuron is in antiphase, and a cycle in which y-neuron and z-neuron are identical while x-neuron is in antiphase.

The time series of two of four coexisting 2-cycles shown by the red and green lines in Figure 2a are presented in Figure 2b,c, respectively. As seen from Figure 2b, all three neurons are completely synchronized because their trajectories coincide ($x=y=z$). One can see that at the first stage ($t=0$), all three neurons are at the same minimum, which we refer to as “rest state”, whereas at the second stage ($t=1$) they are at the same maximum referred to as “active state”. Since in this 2-cycle all neurons fire on and fire off at the same time, the active state is periodically followed by the rest state. We formally denote this synchronization mode as the [03]-pattern. The number of digits in the square brackets indicates the cycle length (2-cycle), and the digit values indicate the number of neurons in the active (firing) mode at each stage, i.e., at the first stage zero neurons are active and in the second stage all three neurons are active.

Figure 2c shows the time series of one of three coexisting partially synchronized 2-cycles (green line in Figure 2a). One can see that in this state x and y are identical, but z-neuron is in antiphase. Due to the symmetry of attractors, this type of synchronization can be formally identified as [21]-pattern, the same for all three attractors. This means that at the first stage, two neurons fire and one neuron is at rest and at the second stage, two neurons are at rest and one neuron fires, and so on. Thus, depending on the initial conditions, the system can be in different synchronization modes representing either [03]- or [21]-patterns.

As the coupling strength $\sigma$ is increased, the number of coexisting attractors also increases. For example, the system (1) with $\sigma =0.4$ exhibits 16 coexisting 4-cycles. The ($x,y$)-projections of these attractors are illustrated in Figure 3a in different colors.

The time series in Figure 3b–f represent five coexisting patterns. Here, we use a similar pattern classification as for the 2-cycle. Namely, the state with maximum y we refer to as “active state”, and other states we call “rest states”. One can see that the pattern in Figure 3b can be classified as [0003]-pattern because at the first three stages ($t=0,1,2$), all three neurons are in the rest state, and at the fourth stage ($t=3$) they fire, i.e., reach maximum. Moreover, this pattern exhibits complete synchronization. Instead, the patterns in Figure 3c–f display partial synchronization representing, respectively, [0111]-, [1002]-, [0102]-, and [2001]-patterns. Other coexisting 2-cycles also exhibit partial synchronization with combination of the above patterns.

As the coupling strength is further increased, chaotic attractors are born in the coupled system (1). The two coexisting 2-piece chaotic attractors are shown in Figure 4a on the ($x,y,z$) plane for $\sigma =0.6$. Here, the attractor in red corresponds to a completely synchronized chaotic state, while the attractor in blue corresponds to a phase synchronized chaotic state. The time series of these attractors are presented in Figure 4b,c, respectively. We should note that it is not possible to use the same classification scheme for chaotic patterns as for periodic cycles because the neurons fire chaotically.

For $\sigma =0.8$, the coupled neurons are in the periodic window (see Figure 1b), where the system (1) has several coexisting stable solutions illustrated in Figure 5. Specifically, in Figure 5a we show 2-cycle (red line) and two 12-cycles (blue and green dots). The time series of these attractors are presented in Figure 5b–d, respectively. As seen from Figure 5b, the 2-cycle demonstrates complete synchronization with the [30]-pattern. At the same time, the blue 12-cycle exhibits phase synchronization with the [010101010101]-pattern (Figure 5c), while the green 12-cycle represents an asynchronous mode with the [011001100110]-pattern (Figure 5d).

For $\sigma =0.9$, the system exhibits the coexistence of two chaotic attractors shown in Figure 6a. They are a 1-piece chaotic attractor (red line) and a 6-piece chaotic attractor (blue dots). The corresponding time series are presented in Figure 6b,c. In the 1-piece attractor, the neurons are completely synchronized, while in the 6-piece attractor they are phase synchronized.

Thus, the deterministic system (1) of three neurons unidirectionally coupled by chemical synapses in a ring exhibits the coexistence of regular and chaotic attractors, which represent different synchronization states: complete synchronization, partial synchronization, phase synchronization, and asynchronous behavior.

In the following section, we will demonstrate how the system behavior changes under random perturbations.

3. Stochastic System

Now, suppose that the three ring-coupled Rulkov neurons undergo independent random perturbations:

$$\begin{array}{c}{x}_{t+1}=f(\gamma ,x_t)-\sigma \phi (x_t,z_t)+\epsilon {\xi }_{1,t}\hfill \\ y_{t+1}=f(\gamma ,y_t)-\sigma \phi (y_t,x_t)+\epsilon {\xi }_{2,t},\hfill \\ z_{t+1}=f(\gamma ,z_t)-\sigma \phi (z_t,y_t)+\epsilon {\xi }_{3,t},\hfill \end{array}$$

(2)

where ${\xi }_{i,t}(i=1,2,3)$ are uncorrelated white Gaussian noises with parameters $\mathrm{E}{\xi }_{i,t}=0$ and $\mathrm{E}{\xi }_{i,t}^2=1$ and $\epsilon$ is the noise intensity.

Let us consider first stochastic effects for $\sigma =0.2$ at which the deterministic system exhibits two patterns, [03] and [21]. Figure 7 shows the time series of the difference $(x-y)$ of the system (2) for different values of the noise amplitude. Under weak noise, the difference $(x-y)$ undergoes small-amplitude oscillations near zero, i.e., the system remains in complete synchronization (Figure 7a). However, stronger noise induces intermittent switching between the completely synchronization regime with [03]-pattern and partial synchronization with [21]-pattern (Figure 7b). One can see that the switching frequency increases as the noise amplitude $\epsilon$ is increased (Figure 7c).

Consider now the stochastic system (2) for $\sigma =0.8$. We remind that the deterministic system is in the periodic window exhibiting regular dynamics with the coexistence of several cycles shown in Figure 5. Let the initial state of the stochastic system (2) be at the deterministic 2-cycle (red dots in Figure 5a). For weak noise ($\epsilon =0.02$), random states leave this cycle and form some distribution around the stable solution (small light blue zones around the red dots in Figure 8a and red lines in Figure 8b). As the noise amplitude is increased ($\epsilon =0.03$), random solutions begin to hop between basins of coexisting attractors and form complex multimodal stochastic oscillations (green zones in Figure 8).

In Figure 8, we see areas where the random states are highly concentrated and areas with a blurred distribution. The concentrated areas correspond to the location of points of the deterministic cycles while the blurred areas appear in zones of chaotic transients. The example of such chaotic transient in the deterministic system (1) is illustrated in Figure 9 by x-time series. As can be seen, the system behaves chaotically before approaching the deterministic 2-cycle.

Note that such stochastic transformations in the dynamics are accompanied by transitions from order to chaos. These transitions occur for a certain noise level, as seen from Figure 10, where we plot the largest Lyapunov exponents versus the noise amplitude for random solutions starting at the points of the completely synchronized 2-cycle. Negative values of $\mathrm{\Lambda }$ localize $\epsilon$-zones of order. As $\epsilon$ exceeds some threshold, the Lyapunov exponents become positive, so that the system dynamics becomes chaotic.

4. Conclusions

The topic of the presented research relates to a new actively developed area of the mathematical modeling and analysis of complex processes in neuron networks. Strong non-linearity of mathematical models and difficulties of obtaining analytical solutions require new methods of computational mathematics, which include bifurcation analysis, description of coexisting attractors and their basins of attraction. In this paper, we have explored complex dynamics of a cyclic ring of Rulkov neurons. Such a ring of three coupled neurons is the elementary cell of a complex neural network, so-called network motif.

Despite its apparent simplicity, the system of three neurons connected by chemical synapses exhibits extremely rich dynamics, including the coexistence of a multitude of attractors. In the deterministic neural model we have found the coexistence of various synchronous regimes: complete synchronization, partial synchronization, and phase synchronization, as well as asynchronous oscillations. As the coupling strength is increased, the system dynamics transforms from stable equilibria to periodic oscillations and chaos.

We have demonstrated that adding stochastic perturbations to each neuron drastically changes the neural dynamics. Specifically, random noise induces intermittent switching between coexisting synchronous and asynchronous regimes. Even relatively weak noise is able to transform the regular synchronous mode to chaotic asynchronous oscillations.

We believe that complex dynamics observed in the small motif of three coupled neurons will help to better understand the functionality of larger neural networks with different topologies. Our study sheds light on the complex mechanisms of the mutual effect of the coupling and random perturbations on the oscillatory behavior of neural systems.