Delay-induced instability and oscillations in a multiplex neural system with Fitzhugh-Nagumo networks

1.   Introduction

Over the past decades, there has been an increasing interest and activity in the investigation of multiplex networks since they can be used to describe many real-life systems, such as brain, social groups, and transport networks, et al. ^[1,2,3]. For example, neural networks are rich in certain subgraphs and the interactions of them are crucial for the proper functioning of the brain. In multiplex systems, the processes happening in one group may vitally affect others and a node in one unit is likely part of another unit ^[1]. The interplay of units can lead to plenty of interesting features, which are often different from the behaviors in isolation. For example, the coupling provides a means for driving and modulating of sustainable synchronous oscillations. In fact, the synchronization of neural activity is an important mechanism to transmit and process information in the brain. The dynamical phenomena of multiplex systems play important roles in their functions and have found extensive applications in various fields. For instance, coexisting attractor are important for image processing or can be taken as an additional source of randomness using for information engineering.

Since the speed of the propagation and processing of signals is always finite, time delays are ubiquitous and non-negligible in nature ^{[2,3,4,5,6,7]}. For example, in neural systems, time delays are up to milliseconds for propagation through the cortical networks. Moreover, time delay can be generated by an autapse, which is a synapse from a neuron onto itself. Time delays can not be simply neglected in the dynamical modeling because they have close relationships with the property and function of system. Time delays can give rise to abundant and interesting phenomena, such as instability, synchronization transitions, and chimera states ^[8,9]. The past few decades have witnessed the rapid development of the dynamics of time-delay systems ^{[2,3,6,10,11,12]}. However, this research area is still open and challenging due to the complexity of the theoretical analysis and simulations of infinite-dimensional dynamical systems.

Recently, great efforts have been paid on the dynamics of multiplex neuronal networks with time delays ^{[2,3,6,13,14,15]}. For instance, Nikitin et al. ^[3] investigated the spatio-temporal dynamics of a multiplex network with two different FitzHugh–Nagumo neuronal loops and delayed inter-layer couplings and revealed coexisting partial synchronization patterns of the two-layer network and presented effective control schemes. In previous studies, two-coupled models with specific topological structure are often considered, such as two-coupled rings ^[6,16]. In biological systems, neurons are often lumped into multiple interconnected networks(areas) ^[17]. Actually, some regions act as a relay between other areas of the brain and play important roles in the signal propagation and brain functionality. For instance, parahippocampal regions can be regarded as relay stations, which actively gate impulse traffic between neocortex and hippocampus, having great effects on the propagation of neural activity ^[2].

Motivated by the above discussions, the objective of this work is to study the behaviors of a multiplex neuronal system with three FitzHugh–Nagumo sub-networks and time delays. Each sub-network consists of an arbitrary number of neurons and intra-connections. These sub-networks are connected through the delayed couplings between a single neuron of each sub-network. An autaptic connection is considered in the network, which is a synapse formed by the axon of a neuron on its own soma or dendrites. The structure of such connection has been observed in cerebral cortex, neocortex, cerebellum, hippocampus, substantia nigra, etc. ^[18,19]. For instance, the majority of cortical pyramidal neurons (more than 80$\%$) in the developing neocortex of a human brain have autaptic connections ^[20]. Actually, the autapse provides a time-delayed self-feedback mechanism and can effectively regulate the membrane potential of neurons.

The remaining part of this paper is organized as follows. In Section 2, the model of the multiplex network is presented and the stability and bifurcation of the system are analyzed. Case studies of numerical simulations are given to validate the obtained results and rich dynamical phenomena are revealed in Section 3. Finally, conclusions are made in Section 4.

2.   Network model and stability

The multiplex system is composed of three FitzHugh-Nagumo neural groups with delayed couplings and one autaptic connection, which can be described by a set of delay differential equations as follows

where $x_i$, $u_j$, and $r_l$ denote the membrane potentials of the $i$-th, $j$-th, $l$-th neuron in the networks $A$, $B$, and $C$, $y_i$, $v_j$, and $s_l$ represent the slow refractory variables of the $i$-th, $j$-th, $l$-th neuron, which describe the time dependence of several physical quantities related to electrical conductances of the relevant ion currents across the membrane, $a_{1i}$, $b_{1i}$, $a_{2j}$, $b_{2j}$, $a_{3l}$, and $b_{3l}$ are positive constants, $d_{im}$, $e_{jp}$, and $h_{lq}$ are the internal connection weights within the sub-networks, $k$ and $\sigma$ are the strength and time delay of the autaptic connection, $c_1$, $c_2$, and $c_3$ represent the strengths of the couplings, $\tau_1$, $\tau_2$, and $\tau_3$ are the coupling time delays between sub-networks, $n_1$, $n_2$, and $n_3$ are the number of the nodes in the sub-networks, $1 \le i, m \le n_1$, $1 \le j, p \le n_2$, $1 \le l, q \le n_3$. In this network, the local kinetics of each node is described by the FitzHugh-Nagumo neuron. The nonlinear activation functions among neurons are assumed to be absolutely smooth and satisfy $f(0) = 0$. Obviously, the origin is the trivial equilibrium of the coupled network. The linearization of Eq (2.1) at the trivial equilibrium of the network can be written in the vector form as follows

where $z = [X^\mathrm{T}, U^\mathrm{T}, R^\mathrm{T}]^\mathrm{T}$, $X = [x_1, y_1, x_2, y_2, \ldots, x_{n_1}, y_{n_1}]^\mathrm{T}$, $U = [u_1, v_1, u_2, v_2, \ldots, u_{n_2}, v_{n_2}]^\mathrm{T}$, $R = [r_1, s_1, r_2, s_2, \ldots, r_{n_3}, s_{n_3}]^\mathrm{T}$, $\beta_{im} = d_{im}f'(0)$, $\mu_{jp} = e_{jp}f'(0)$, $\eta_{lq} = h_{lq}f'(0)$, $\gamma_{g} = c_{g}f'(0)$, $g = 1, 2, 3$,

$M_1 = \begin{bmatrix} A & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & C \end{bmatrix}$, $M_2 = \begin{bmatrix} 0 & 0 & 0 \\ E_1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$, $M_3 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & E_2 & 0 \end{bmatrix}$, $M_4 = \begin{bmatrix} 0 & 0 & E_3 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$, $M_5 = \begin{bmatrix} E_4 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$, $E_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}_{2n_2\times 2n_1}$, $E_2 = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}_{2n_3\times 2n_2}$, $E_3 = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}_{2n_1\times 2n_3}$, $E_4 = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}_{2n_1\times 2n_1}$, $A = \begin{bmatrix} A_1+\beta_{11}G & \beta_{12}G & \cdots & \beta_{1n_1}G \\ \beta_{21}G & A_2+\beta_{22}G & \cdots & \beta_{2n_1}G \\ \vdots & \vdots & \ddots & \vdots \\ \beta_{n_11}G & \beta_{n_12}G & \cdots & A_{n_1}+\beta_{n_1 n_1}G \end{bmatrix}$, $B = \begin{bmatrix} B_1+\mu_{11}G & \mu_{12}G & \cdots & \mu_{1n_2}G \\ \mu_{21}G & B_2+\mu_{22}G & \cdots & \mu_{2n_2}G \\ \vdots & \vdots & \ddots & \vdots \\ \mu_{n_21}G & \mu_{n_22}G & \cdots & B_{n_2}+\mu_{n_2 n_2}G \end{bmatrix}$, $C = \begin{bmatrix} C_1+\eta_{11}G & \eta_{12}G & \cdots & \eta_{1n_3}G \\ \eta_{21}G & C_2+\eta_{22}G & \cdots & \eta_{2n_3}G \\ \vdots & \vdots & \ddots & \vdots \\ \eta_{n_31}G & \eta_{n_32}G & \cdots & C_{n_3}+\eta_{n_3 n_3}G \end{bmatrix}$, $A_i = \begin{bmatrix} a_{1i} & -1 \\ 1 & -b_{1i}\end{bmatrix}$, $B_j = \begin{bmatrix} a_{2j} & -1 \\ 1 & -b_{2j}\end{bmatrix}$, $C_l = \begin{bmatrix} a_{3l} & -1 \\ 1 & -b_{3l}\end{bmatrix}$, $G = \begin{bmatrix} 1 & 0 \\ 0 & 0\end{bmatrix}$.

After some calculations, the characteristic equation of the network reads

where $P_1(\lambda) = |\lambda I_1-A|$, $P_2(\lambda) = |\lambda I_2-B|$, $P_3(\lambda) = |\lambda I_3-C|$, $\widetilde{P}_1(\lambda) = |\lambda \widetilde{I}_1-\widetilde{A}|$, $\widetilde{P}_2(\lambda) = |\lambda \widetilde{I}_2-\widetilde{B}|$, $\widetilde{P}_3(\lambda) = |\lambda \widetilde{I}_3-\widetilde{C}|$, $P(\lambda) = P_1(\lambda)P_2(\lambda)P_3(\lambda)$, $L(\lambda) = P(\lambda)+k(\lambda+b_{11})\widetilde{P}_1(\lambda)P_2(\lambda)P_3(\lambda)$, $H(\lambda) = (\lambda+b_{11})\widetilde{P}_1(\lambda)P_2(\lambda)P_3(\lambda)$, $W(\lambda) = (\lambda+b_{11})(\lambda+b_{21})(\lambda+b_{31})\widetilde{P}_1(\lambda)\widetilde{P}_2(\lambda)\widetilde{P}_3(\lambda)$, $I_1$, $I_2$, $I_3$, $\widetilde{I}_1$, $\widetilde{I}_2$, and $\widetilde{I}_3$ are identity matrice, $\gamma = \gamma_1\gamma_2\gamma_3$, $\tau = \tau_1+\tau_2+\tau_3$, $\widetilde{A}$, $\widetilde{B}$, and $\widetilde{C}$ represent the connection matrix of three sub-networks without the first neuron respectively, i.e., $\widetilde{A} = \begin{bmatrix} A_2+\beta_{22}G & \beta_{23}G & \cdots & \beta_{2n_1}G \\ \beta_{32}G & A_3+\beta_{33}G & \cdots & \beta_{3n_1}G \\ \vdots & \vdots & \ddots & \vdots \\ \beta_{n_12}G & \beta_{n_13}G & \cdots & A_{n_1}+\beta_{n_1 n_1}G \end{bmatrix}$, $\widetilde{B} = \begin{bmatrix} B_2+\mu_{22}G & \mu_{23}G & \cdots & \mu_{2n_2}G \\ \mu_{32}G & B_3+\mu_{33}G & \cdots & \mu_{3n_2}G \\ \vdots & \vdots & \ddots & \vdots \\ \mu_{n_22}G & \mu_{n_23}G & \cdots & B_{n_2}+\mu_{n_2 n_2}G \end{bmatrix}$, $\widetilde{C} = \begin{bmatrix} C_2+\eta_{22}G & \eta_{23}G & \cdots & \eta_{2n_3}G \\ \eta_{32}G & C_3+\eta_{33}G & \cdots & \eta_{3n_3}G \\ \vdots & \vdots & \ddots & \vdots \\ \eta_{n_32}G & \eta_{n_33}G & \cdots & C_{n_3}+\eta_{n_3 n_3}G \end{bmatrix}$.

According to the characteristic equation of the network, the local stability of the trivial equilibrium of the system is determined by the root distributions of the characteristic equations of each individual sub-network, the coupling strengths and time delays. The stability analysis begins with the case when autaptic time delay is zero, i.e., $\Delta(\lambda, 0, \tau) = P(\lambda)-\gamma W(\lambda)e^{-\lambda \tau} = 0$. The necessary and sufficient conditions for all roots of $\Delta(\lambda, 0, 0) = 0$ having negative real parts can be determined by the Routh–Hurwitz criteria. As the coupling time delay varies, let $\lambda = \pm i\upsilon(\upsilon > 0)$ be a pair of purely imaginary roots. Then, one arrives at

where $P_R(\upsilon) = Re[P(i\upsilon)]$, $P_I(\upsilon) = Im[P(i\upsilon)]$, $W_R(\upsilon) = Re[W(i\upsilon)]$, and $W_I(\upsilon) = Im[W(i\upsilon)]$. Clearly, $P_R(\upsilon)$ and $W_R(\upsilon)$ are even functions, whereas $P_I(\upsilon)$ and $W_I(\upsilon)$ are odd functions. Separating the real and imaginary parts of $\Delta(i\upsilon, 0, \tau) = 0$ and eliminating the harmonic terms give $D(\upsilon) = P_R^2(\upsilon)+P_I^2(\upsilon)-\gamma^2[W_R^2(\upsilon)+W_I^2(\upsilon)] = 0$. Suppose that a positive root $\upsilon_j$ of $D(\upsilon) = 0$ can be found. Then, a set of critical coupling time delay is given by $\tau_{j, n} = (\theta_j+2n\pi)/\upsilon_j$, $n = 0, 1, 2, …$ where $\theta_j\in\left[0, 2\pi \right)$ yields sets of triangle equations

If the polynomial $D(\upsilon) = 0$ has no positive roots, the above critical time delays do not exist. In this case, the network is delay-independent stable or unstable for any given coupling time delay, depending on whether or not the system free of time delay is stable.

When the autaptic time delay is taken into account i.e., $\sigma > 0$, it can be regarded as a parameter for fixed value of the coupling time delay. In order to find the boundary of stability, one needs to consider the case when the characteristic equation has a pair of purely imaginary roots $\lambda = \pm i\omega(\omega > 0)$. Then, one obtains

where $L_R(\omega) = Re[L(i\omega)]$, $L_I(\omega) = Im[L(i\omega)]$, $H_R(\omega) = Re[H(i\omega)]$, $H_I(\omega) = Im[H(i\omega)]$, $W_R(\omega) = Re[W(i\omega)]$, and $W_I(\omega) = Im[W(i\omega)]$. Clearly, $L_R(\omega)$, $H_R(\omega)$, and $W_R(\omega)$ are even functions, whereas $L_I(\omega)$, $H_I(\omega)$, and $W_I(\omega)$ are odd functions. Separating the real and imaginary parts of the characteristic equation with two different delays yields

Eliminating the autaptic time delay in the above equation gives

If $F(\omega) = 0$ has any positive roots $\omega_j$, then the characteristic equation has purely imaginary roots $\lambda = \pm i\omega$. In this case, the characteristic equation of the network has sets of critical autaptic time delays $\sigma_{j, n} = (\varphi_j+2n\pi)/\omega_j$, $j = 1, 2, \ldots$, $n = 0, 1, 2, \ldots$, where $\varphi_j\in\left[0, 2\varphi\right)$ and $\varphi_j$ satisfies

After some calculations, one arrives at $Re[\lambda(\sigma)|_{\lambda = i\omega}]' = 0.5\omega F'(\omega)/|S(i\omega)|^2$, where, $S(i\omega) = L'(i\omega)-kH'(i\omega)e^{-i\omega\sigma}+k\sigma H(i\omega)e^{-i\omega \sigma}-\gamma W'(i\omega)e^{-i\omega\tau}+\gamma \tau W(i\omega)e^{-i\omega\tau}$. Thus, the variation direction of its real part with respect to the autaptic time delay can be determined by the sign of $F'(\omega)$.

Based on the above results, one can check that the crossing real parts of the roots of the characteristic equation at $\sigma_{2j-1, n}$ corresponding to $\pm i\omega_{2j-1}$ must be from the left to the right, and the crossing at $\sigma_{2j, n}$ corresponding to $\pm i\omega_{2j}$ must be from the right to the left, where $\omega_1 > \omega_2 > \cdots > \omega_j > \omega_{j+1} > \cdots > 0$, $j = 1, 2, \ldots$. Thus, as the internal time delay $\sigma$ varies from zero to the infinity, the characteristic equation always adds a new pair of conjugate roots with positive real parts for each crossing at $\sigma_{2j-1, n}$, but reduces such a pair for each crossing at $\sigma_{2j, n}$. In addition, more roots of the characteristic equation change their sign of real parts from the negative to the positive at $\sigma_{2j-1, n}$ than those changing the sign of real parts from the positive to the negative at $\sigma_{2j, n}$ with an increase of autaptic time delay in a given long interval. Then, the system can undergo finite stability switches and become unstable at last when the autaptic time delay increases from zero to the infinity ^[5]. On the other hand, when Equation (2.8) has no positive roots, the stability of the system is independent of the autaptic time delay because the signs of the real parts of the characteristic roots do not change when the autaptic time delay varies.

3.   Illustrative examples

In this section, the nonlinear activation function between neurons is chosen as hyperbolic tangent function, which is a typical sigmoid function and has been widely used in neuronal networks.

(1) $n_1 = n_2 = n_3 = 3$, $a_{11} = 0.8$, $a_{21} = a_{31} = 0.3$, $a_{12} = a_{22} = a_{32} = 0.4$, $a_{13} = a_{23} = a_{33} = 0.5$, $b_{11} = b_{21} = b_{31} = 0.8$, $b_{12} = b_{22} = b_{32} = 0.5$, $b_{13} = b_{23} = b_{33} = 1.2$, $k = 0.5$, $d_{21} = d_{32} = d_{13} = 0.15$, $e_{21} = e_{32} = e_{13} = 0.15$, $h_{21} = h_{32} = h_{13} = 0.15$, $c_1 = c_2 = c_3 = -0.2$, $\tau_1 = \tau_2 = \tau_3 = \tau_s$, and other parameters are zero. After some calculations, one can check that the trivial equilibrium of the multiplex network free of time delays is locally asymptotically stable. Solving the polynomial $D(\upsilon) = 0$ gives $\upsilon_1 = 0.615$ and $\upsilon_2 = 0.563$. There follows, two sets of the critical coupling time delays can be obtained $\tau_{1, n} = 0.65, 10.87, 21.09, 31.32, 41.54, \ldots$ and $\tau_{2, n} = 7.25, 18.42, 29.58, 40.75, 51.92, \ldots$. Then, a pair of roots of the characteristic equation is crossing the imaginary axis from the left to the right when $\tau = \tau_{1, n}$ and from the right to the left when $\tau = \tau_{2, n}$. Thus, the trivial equilibrium of the multiplex network free of autaptic time delay is locally asymptotically stable for $\tau\in\left[0, \tau_{1, 0}\right) \cup\cdots\cup(\tau_{2, n}, \tau_{1, n+1})\cup\cdots\cup(\tau_{2, 3}, \tau_{1, 4})$ and becomes unstable for $\tau\in(\tau_{1, 0}, \tau_{2, 0})\cup\cdots\cup(\tau_{1, n}, \tau_{2, n})\cup\cdots\cup(\tau_{1, 4}, +\infty)$.

Figure 1(a) gives the stable trivial equilibrium of the network when $\tau_s = 0.1$. Figure 1(b) shows that periodic oscillations arising from Hopf bifurcation come into being when $\tau_s = 1$. Figure 1(c) illustrates that the trivial equilibrium remain stable since all roots of the characteristic equation have negative real parts when $\tau_s = 3$. Figure 1(d) gives that the trivial equilibrium loses its stability again and periodic oscillations occur because the network adds a pair of characteristic roots with positive real parts when $\tau_s = 4$. The solid, dashed, and dotted curves are the responses of the first neurons in each sub-network. As shown in Figure 1, the system exhibits stable rest state and bifurcated periodic oscillations, which coincide with the theoretical results.

Figure 1.  The responses of the multiplex network free of autaptic time delays.

DownLoad: Full-Size Img PowerPoint

Let $\tau_s = 10$. After some calculations, the quasi-polynomial $F(\omega) = 0$ has two positive real roots $\omega_1 = 1.09$ and $\omega_2 = 0.62$. There follows, two sets of critical autaptic time delays can be obtained $\sigma_{1, n} = 4.62, 10.36, 16.11, \ldots$ and $\sigma_{2, n} = 10.07, 20.20, 30.34, \ldots$. Based on the conclusions in the above section, a pair of roots of the characteristic equation is crossing the imaginary axis from the left to the right when $\sigma = \sigma_{1, n}$, and from the right to the left when $\sigma = \sigma_{2, n}$. The critical autaptic time delays can be ranked as $0 < \sigma_{1, 0} < \sigma_{2, 0} < \sigma_{1, 1} < \sigma_{1, 2} < \ldots$. Hence, the trivial equilibrium of the multiplex neuronal network is locally asymptotically stable for $\sigma\in\left[0, \sigma_{1, 0}\right) \cup(\sigma_{2, 0}, \sigma_{1, 1})$ and becomes unstable for $(\sigma_{1, 0}, \sigma_{2, 0})\cup(\sigma_{1, 1}, +\infty)$. It is obvious that the trivial equilibrium of the multiplex network undergoes three stability switches and will become unstable at last. Figure 2 gives the responses of the multiplex network when the autaptic time delay varies. As shown in Figure 2, the network exhibits the switches between the stable trivial state and periodic oscillations, which reach an agreement with the obtained analytical conclusions.

Figure 2.  The responses of the multiplex network when $\tau_s = 10$.

DownLoad: Full-Size Img PowerPoint

(2) $n_1 = n_2 = n_3 = 3$, $a_{11} = a_{21} = a_{31} = 1$, $a_{12} = a_{22} = a_{32} = 1.1$, $a_{13} = a_{23} = a_{33} = 1.3$, $b_{11} = b_{21} = b_{31} = 0.65$, $b_{12} = b_{22} = b_{32} = 1.12$, $b_{13} = b_{23} = b_{33} = 0.8$, $k = 0.2$, $d_{21} = d_{32} = d_{13} = 0.15$, $e_{21} = e_{32} = e_{13} = 0.15$, $h_{21} = h_{32} = h_{13} = 0.15$, $c_1 = c_2 = c_3 = 0.15$, $\sigma = 0.05$, $\tau_1 = \tau_2 = \tau_3 = \tau_s$, and other parameters are zero.

A Poincaré section is defined as the projection of solutions of the system. The points on Poincaré section depend on the behaviors of the system. If the final motion of the system is periodic, there is only one point on the Poincaré section. For a period-n motion, n points will appear in the section, but the numbers of points become infinite for non-periodic motions such as chaotic responses. The Poincaré section is defined by $S = \left\lbrace (x_1(t), x_1(t-1)):(y_1(t) = 0, \overset{.}y_1(t) > 0)\right\rbrace$. Figure 3 gives the responses of the multiplex network when $\tau_s = 0.1$. The phase plots depict that a pair of period-3 orbits and a pair of chaotic motions come into being under different initial conditions, as shown in Figures 3(a), 3(b), 3(d), and 3(e). The round and asterisk points appear in Figures 3(c) and 3(f) and verify these coexisting period-3 orbits and chaotic attractors.

Figure 3.  The responses of the multiplex network when $\tau_s = 0.1$. (a) and (b) Phase trajectories of a pair of period-3 oscillations under initial conditions IC1(0.1, 0.4, 0.8, 0.8, 0.8, 0.7, -0.3, -0.75, -0.25, -0.9, -1.25, -0.05, -0.01, -0.6, -0.45, -0.5, -1, -0.2) and IC2(-0.1, -0.4, -0.8, -0.8, -0.8, -0.7, 0.3, 0.75, 0.25, 0.9, 1.25, 0.05, 0.01, 0.6, 0.45, 0.5, 1, 0.2); (c) Poincaré section plot; (d) and (e) Phase trajectories of a pair of chaotic attractors under initial conditions IC3(0.1, 0.4, 0.8, 0.8, 0.8, 0.7, 0.3, 0.75, 0.25, 0.9, 1.25, 0.05, -0.01, -0.6, -0.45, -0.5, -1, -0.2) and IC4(-0.1, -0.4, -0.8, -0.8, -0.8, -0.7, -0.3, -0.75, -0.25, -0.9, -1.25, -0.05, 0.01, 0.6, 0.45, 0.5, 1, 0.2); (f) Poincaré section plot.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the coexistence of a pair of period-3 orbits and a pair of period-2 oscillations when $\tau_s = 0.3$. The points on the Poincaré plots, as shown in Figures 4(c), and 4(f), coincide with the results in the phase plane. Figure 5 illustrates the multiple coexisting chaotic motions when $\tau_s = 0.4$. As shown in Figures 3–5, the network exhibits different types of multiple coexisting attractors when the coupling time delays vary. That is, the time delays have great effects on the performance of the multiplex network and can be used to regulate the complicated dynamical behaviors of the system, such as multi-stability.

Figure 4.  The responses of the multiplex network when $\tau_s = 0.3$. (a) and (b) Phase trajectories of a pair of period-3 orbits under initial conditions IC1 and IC2; (c) Poincaré section plot; (d) and (e) Phase trajectories of a pair of period-2 motions under initial conditions IC3 and IC4, respectively; (f) Poincaré section plot.

DownLoad: Full-Size Img PowerPoint

Figure 5.  The responses of the multiplex network when $\tau_s = 0.4$. (a) and (b) Phase trajectories of two chaotic responses under initial conditions IC1 and IC2; (c) Poincaré section plot; (d) and (e) Phase trajectories of two chaotic orbits under initial conditions IC3 and IC4; (f) Poincaré section plot.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

This paper has studied the stability, bifurcation, and multi-stability coexistence of a FitzHugh–Nagumo neuronal network consisting of three populations with delayed couplings between one neuron of each sub-network and autapstic connection. By regarding the sum of coupling time delays and autaptic time delay as the parameters, the stability and bifurcation of the network equilibrium have been studied. It is shown that the stability of the network can be determined by the root distributions of characteristic equations of each individual sub-network, the product of coupling strengths, autaptic time delay, and the sum of coupling time delays. Illustrative examples are given to validate the theoretical analysis and a variety of complex phenomena are observed, such as the coexistence of different multi-period orbits and chaotic motions. The revealed results may provide promising and useful information for understanding the mechanisms of rhythms and complexity of real interacting neural systems.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant Nos. 11872169 and 12172119, Natural Science Foundation of Jiangsu Province under Grant No. BK20191295, and Fundamental Research Funds for the Central Universities under Grant No. B200202114. The authors thank the anonymous reviewers for their helpful comments and suggestions that have helped to improve the presentation.

Conflict of interest

The authors declare there is no conflict of interest.