Complexity in time-delay networks of multiple interacting neural groups

1.   Introduction

The last few years have witnessed the rapid development of multiplex networks due to their important and extensive applications in many fields, such as brain, Internet, power grids, food webs, social networks, transport networks, etc[1,2,3,10,13,15]. Multiplex networks consist of interacting individual sub-networks with couplings between them. For example, in the brain, neurons are often lumped into connected groups and the interactions of them are crucial for the proper functioning of the brain[1]. Due to the finite propagation speeds of signals, time delays are inevitable phenomena in mulitplex networks[13,15]. For example, axonal conduction could result in time delays up to milliseconds for the singnal propagation between areas in the brain. As a matter of fact, neglecting time delays usually leads to false or even wrong results[6,9,14].

The analysis of the dynamics of multiplex networks has close relationships with their practical applications. In biological science, for example, the network of coupled cells (each cell is composed of a set of neurons) can exhibit the quadruped gaits of walk, trot, and pace. In particular, by adjusting the coupling strengths and time delays, the primary gaits can be produced from Hopf bifurcation. Experiments show that abnormal neural synchronous oscillations in interacting neural loops are closely related to some diseases, such as schizophrenia, epilepsy, autism, Parkinson's disease and Alzheimer's disease[3,4]. Chaotic motions have been shown in neural systems, which associate with cognition, learning, and retrieval. Multi-stability indicates the possibility to have several attractors for the same parameters when the initial conditions vary. It has also been observed in single neuron and neuronal populations. For instance, in cortical neurons, the transitions between spiking and bursting corresponded to a function of sleep and wakefulness[5].

In interacting networks, the dynamic processes happening in one network may vitally affect another network since a node in one network is likely part of another network[10]. The interplay of networks can lead to plenty of interesting dynamical features. For example, the couplings can provide a means for driving and modulating of sustainable oscillations. Actually, the behaviors of multiplex networks are often different from those observed in isolated network. Time delays, even very short in some cases, have significant influences on the properties of systems and cannot be simply neglected in the dynamical analysis. It is found that time delay can be used as a powerful tool to control various dynamical features, such as synchronization patterns and their transitions, quasi-periodic and chaotic motions, chimera states, et al[6,9,11,14]. For instance, in a two-ring network, the delayed interactions induce complicated waves and the appropriate time delay and transverse coupling can achieve the "backup and recovery" procedure[15]. Thus, the performances of multiplex networks are governed by the assumed dynamics of individual sub-network, their reciprocal interactions and time delays.

Multiplex structures are ubiquitous in biological and physiological systems, such as neural networks. In recent years, the dynamics of multiplex neural networks with time delays have drawn much attention. Nikitin et al.[12] investigated the spatio-temporal dynamics of a multiplex network with two different neuronal loops and delayed inter-layer couplings and revealed coexisting partial synchronization patterns of the two-layer network and presented effective control schemes. Sawicki et al.[13] studied the relay synchronization in a three-layer multiplex FitzHugh-Nagumo neural network with nonlocal coupling topologies within the layers and considered time delay in the interlayer coupling. For more related studies on the dynamics of multiplex neural networks with time delays, the reader is referred to the work in[3,8,15], and the references cited therein. However, few efforts are made on multiple coupled neural networks with time delays. In neural systems, many brain areas, consisting of neural assemblies, are coupled by physical or functional interactions among them. In functional brain networks, for instance, each link indicates a functional interaction between the activity (electrical, magnetic, or hemodynamic/metabolic) of two areas [1]. Particularly, some specific brain areas act as relays between other brain regions, playing important roles in brain functionality and dysfunctions, such as parahippocampal regions[13].

Motivated by the above discussions, the purpose of this paper is to study the dynamical behaviors of a multiplex network that consists of multiple Hopfield neural groups and time delays. Each sub-network can contain different sets of neurons and topological structures. The multiplex network is connected by the couplings between one neuron of each sub-network. The intra-layer time delays are neglected compared with the internal time scale of the neural system, while the inter-layer coupling time delay is larger and cannot be neglected. Then, different time delays are introduced into the couplings between networks.

The rest of this paper is organized as follows. In Section 2, the model of the network is presented and the stability and bifurcations of the network are analyzed. Case studies of illustrative examples are given to support the theoretical results and abundant neural activities are shown in Section 3. Finally, conclusions are drawn in Section 4.

2.   Stability and bifurcation analysis

The multiplex network consists of multiple interacting sub-networks with delayed couplings, which can be described by the following delay differential equations

where $u_p$ and $b_{pq}$ are the states of the neurons and connection weights in the $k$-th sub-network for $N_{k-1}+1 \le p$, $q \le N_k$, $f$ represents the nonlinear activation function of the neurons, $c_k$, $\tau_k$, and $g$ are the strengths, time delays and transfer function of the couplings between sub-networks, $n_k$ is the number of the neurons in the $k$-th sub-network, $m$ is the number of sub-networks, $N_m = \sum_{k = 1}^{m}n_k$ is the total number of the neurons of the network, $k = 1, 2, ..., m$. In this network, the local kinetics of each node is described by the Hopfield neuron [7]. Without loss of generality, the nonlinear functions are assumed to be absolutely smooth and satisfy $f_i(0) = g(0) = 0$. The system may have other non-trivial equilibria. The analysis for other equilibrium positions is similar to the case of the trivial equilibrium position after a simple change of variables. The linearization of (1) at the trivial equilibrium of the network can be written in the vector form as follows

where $u(t) = (u_1, u_2, \ldots, u_{N_m})^\mathrm{T}$, $\alpha_{pq} = b_{pq}f'(0)$, $\beta_k = c_kg'(0)$, $A_k = (\alpha_{pq})$, $N_{k-1}+1 \le p$, $q \le N_k$, $k = 1, 2, ..., m$, $p$, $q = 1, 2, ..., N_m$, and

$A = \begin{bmatrix} A_1 & 0 & \cdots & 0 & 0 \\ 0 & A_2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & A_{m-1} & 0 \\ 0 & 0 & \cdots & 0 & A_{m} \end{bmatrix}$, $L_1 = \begin{bmatrix} 0 & 0 & \cdots & 0 & 0 \\ C_1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & 0 \end{bmatrix}$,

$L_2 = \begin{bmatrix} 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 \\ 0 & C_2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & 0 \end{bmatrix}$, $\cdots\cdots$, $L_{m-1} = \begin{bmatrix} 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & C_{m-1} & 0 \end{bmatrix}$,

$L_m = \begin{bmatrix} 0 & 0 & \cdots & 0 & C_m \\ 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & 0 \end{bmatrix}$, $C_k = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}_{n_{k+1} \times n_k}$

The characteristic matrix of the system is

where $I_k$ is a $n_k \times n_k$ identity matrix, $k = 1, 2, \ldots, m$. Then, the characteristic equation of the network reads

where $G_k(\lambda) = |(\lambda+1)I_k-A_k| = \prod_{h = 1}^{n_k}(\lambda-s_{k, h})$, $H_k(\lambda) = |(\lambda+1)\widetilde{I}_k-\widetilde{A}_k| = \prod_{l = 1}^{n_k-1}(\lambda-\widetilde{s}_{k, l})$, $G(\lambda) = \prod_{k = 1}^mG_k(\lambda)$, $H(\lambda) = \prod_{k = 1}^mH_k(\lambda)$, $\widetilde{I}_k$ is a $(n_k-1) \times (n_k-1)$ identity matrix, $\beta = \prod_{k = 1}^m\beta_k$, $\tau = \sum_{k = 1}^m\tau_k$, and $\widetilde{A}_k = (\alpha_{pq})$, $N_{k-1}+1 \le p$, $q \le N_{k}-1$. In fact, $A_k$ and $s_{k, h}$ represent the connection matrix and the $h$-th eigenvalue of the $k$-th sub-network, while $\widetilde{A}_k$ and $\widetilde{s}_{k, r}$ are the connection matrix and the $r$-th eigenvalue of the $k$-th sub-network without the last neuron. The trivial equilibrium of the multiplex network is locally asymptotically stable when all roots of the characteristic equation have negative real parts, and unstable when at least one root of the characteristic equation has a positive real part. Actually, a steady state bifurcation occurs when $\beta = \prod s_{k, h}/\prod\widetilde{s}_{k, l}$ since $\lambda = 0$ is the root of $\Delta(\lambda, \tau) = 0$, where $k = 1, 2, \ldots, m$, $h = 1, 2, \ldots, n_k$, $l = 1, 2, \ldots, n_k-1$. In addition, a Hopf bifurcation occurs when $\Delta(\lambda, \tau) = 0$ has a pair of pure imaginary roots. To determine the boundary of stability, the analysis begins with the case when the network has no time delays, that is,

where $G(\lambda) = \lambda^n+g_1\lambda^{n-1}+g_2\lambda^{n-2}+\cdots+g_{n-1}\lambda+g_n$, $H(\lambda) = \lambda^{n-m}+h_1\lambda^{n-m-1}+h_2\lambda^{n-m-2}+\cdots+h_{n-m-1}\lambda+h_{n-m}$, $n = N_m$. For the case when all characteristic roots $s_{k, h}$ of each sub-network and all characteristic roots $\widetilde{s}_{k, l}$ of each sub-network without the last neuron have negative real parts and $\beta < 0$, all coefficients $e_i > 0$ hold. It follows that $\Delta(\lambda, 0) > 0$ holds for $\lambda > 0$. Thus, the multiplex network free of time delays is stable. In this case, the trivial equilibrium of the network free of time delays is stable when the isolated sub-networks and the isolated sub-networks without the last neuron are stable and the number of the inhibitory couplings between sub-networks is odd. In general case, a set of necessary and sufficient conditions for all roots of (4) having negative real parts is given by Routh-Hurwitz criteria. As the time delay varies, let $\lambda = \pm i\omega(\omega > 0)$ be the roots of the characteristic equation of the multiplex network, one arrives at

where $G_R(\omega) = Re[G(i\omega)] = Re[\prod_{k = 1}^mG_k(i\omega)]$, $G_I(\omega) = Im[G(i\omega)] = Im[\prod_{k = 1}^m G_k(i\omega)]$, $H_R(\omega) = Re[H(i\omega)] = Re[\prod_{k = 1}^mH_k(i\omega)]$ and $H_I(\omega) = Im[H(i\omega)] = Im[\prod_{k = 1}^mH_k(i\omega)]$. By separating the real and imaginary parts and eliminating the harmonic terms, one arrives at $F(\omega) = G_R^2(\omega)+G_I^2(\omega)-\beta^2[H_R^2(\omega)+H_I^2(\omega)] = 0$. If $F(\omega) = 0$ has any positive roots $\omega_j$, a set of critical time delays can be determined by $\tau_{j, r} = (\theta_j+2r\pi)/\omega_j$, r = 0, 1, 2, ... where $\theta_j\in \left[ 0, 2\pi\right)$ yields sets of triangle equations $cos\theta_j\! = \![G_R(\omega_j)H_R(\omega_j)+G_I(\omega_j)H_I(\omega_j)]/[\beta|H(i\omega_j)|^2]$ and $sin\theta_j\! = \![G_R(\omega_j)H_I(\omega_j)-G_I(\omega_j)H_R(\omega_j)]/[\beta|H(i\omega_j)|^2]$. Then, the trivial equilibrium of the multiplex network will undergo finite stability switches and becomes unstable at last when time delay increases from the zero to the positive infinity [9]. On the other hand, if $F(\omega) = 0$ has no positive root, then the network is delay-independent stable or unstable for any given time delay, depending on whether or not the multiplex network free of time delay is stable.

Based on the above results, the stability of the trivial equilibrium of the network can be obtained by the characteristic values $s_{k, h}$ and $\widetilde{s}_{k, l}$ of the isolated sub-networks and the product of the coupling strength $\beta$ and the sum of the time delays $\tau$. Thus, the stability of the network can be determined by the characteristic values of the sub-networks and the strengths and time delays of the couplings.

Let $m = 3$, $n_1 = 2$, $n_2 = 3$, $n_3 = 4$, $b_{11} = 0.8$, $b_{12} = -1.8$, $b_{21} = 1.3$, $b_{22} = -1.5$, $b_{33} = 0.8$, $b_{34} = -3.8$, $b_{35} = -1.8$, $b_{43} = -0.2$, $b_{44} = -1.8$, $b_{45} = 2.3$, $b_{53} = 1.3$, $b_{54} = 2.1$, $b_{55} = -1.5$, $b_{66} = 0.8$, $b_{67} = -3.8$, $b_{68} = 0.5$, $b_{69} = -1.8$, $b_{76} = -0.2$, $b_{77} = -1.8$, $b_{78} = -1$, $b_{79} = 2.3$, $b_{86} = -0.3$, $b_{87} = 1.2$, $b_{88} = -1.8$, $b_{89} = 0.1$, $b_{96} = 1.3$, $b_{97} = 2.1$, $b_{98} = 0.1$, and $b_{99} = -1.5$. In this case, the network consists of three sub-networks with different nodes. The numbers of the nodes in the three sub-networks are two, three, and four. Figure 1 gives the region indicating stability of the origin of the delayed triplex network. As shown in Figure 1, by calculating the values of the critical time delays for fixed $\beta$ and then changing the values of $\beta$, the black and red bifurcation curves can be obtained. Along the solid black and red bifurcation curves, the characteristic equation of the system has a pair of imaginary eigenvalues. More specifically, the characteristic equation of the system adds a pair of conjugate roots with positive real parts for each crossing the black curves, but reduces a pair for each crossing the red curves. For example, when $\beta = -1.4$, as $\tau$ increases from zero to infinity, the trivial equilibrium of the network loses its stability for crossing the black curves and remains stability for crossing the red curves, but it becomes unstable at last.

Let $m = 4$, $n_k = 3$, $b_{ij} = a(i\neq j)$, $b_{ii} = 0$, $c_k = c$, and $\tau_k = \tau_s$. This network consists of four neural loops each with three nodes and bidirectional interconnections. The characteristic equation of the network yields $\Delta(\lambda, \tau) = G_1^4(\lambda)-\beta e^{-\lambda\tau}H_1^4(\lambda) = 0$, where $G_1(\lambda) = (\lambda +1-2\alpha)(\lambda +1+\alpha)^2$ and $H_1(\lambda) = (\lambda+1)^2-\alpha^2$. Thus, the characteristic equation can be rewritten as

where $\Delta_1(\lambda, \tau) = (\lambda +1+\alpha)^4$, $\Delta_2(\lambda, \tau) = [(\lambda +1-2\alpha)(\lambda +1+\alpha)-\beta_1e^{-\lambda\tau_s}(\lambda +1-\alpha)]$, $\Delta_3(\lambda, \tau) = [(\lambda +1-2\alpha)(\lambda +1+\alpha)+\beta_1e^{-\lambda\tau_s}(\lambda +1-\alpha)]$, and $\Delta_4(\lambda, \tau) = [(\lambda +1-2\alpha)^2(\lambda +1+\alpha)^2+\beta_1^2e^{-2\lambda\tau_s}(\lambda +1-\alpha)^2]$. It is obvious that the network undergoes $a$ steady state bifurcation when $a = -1$ and $1-\alpha-2\alpha^2\mp\beta_1\pm\alpha\beta_1 = 0$. Figure 2 gives the stable region of the trivial equilibrium of the network in the parameter plane of the intra-connection weight $\alpha$ of the sub-network and coupling strength $\beta_1$ between sub-networks.

As shown in Figure 2, there exists only delay-independent stability region, where the trivial equilibrium is always locally stabile for all time delays. In particular, one can check that the distribution of the characteristic roots of the four-coupled network is the same as the two-coupled model consisting of a pair of bidirectional rings. Then, the grey region shown in Figure 2 is also the stable area of the two-coupled ring neural networks.

In neural systems, due to the input-output relation of neurons, the activation functions $f$ and $g$ are often assumed to satisfy $f'(0) = g'(0) = 1$; $f'(\upsilon) > 0$ and $g'(\upsilon) > 0$ for all $\upsilon\in\mathbb{R}$; $-\infty< \lim_{\upsilon\to\pm\infty}f(\upsilon), g(\upsilon)< +\infty$; $\upsilon f''(\upsilon)< 0$ and $\upsilon g''(\upsilon)< 0$ for all $\upsilon\neq0$. Typically, the sigmoid functions, such as hyperbolic tangent function, are widely used in neural networks. By constructing a Lyapunov functional $V(t) = \sum_{p = 1}^{N_m}u_p^2(t)+\sum_{j = 1}^m|c_j|\int_{t-\tau_j}^tg^2[u_{N_j}(s)]ds$, one has the following derivatives

where $b = max\left\lbrace|b_{pq}| \right\rbrace$, $c = max\left\lbrace|c_{j}| \right\rbrace$, and $n_r = max\left\lbrace|n_k| \right\rbrace$. It is easy to check that $f'(\upsilon)\le 1$ and $g'(\upsilon)\le 1$. Moreover, $f[\upsilon_i(t)] = \varphi_i(t)\upsilon_i(t)$, $g[\upsilon_i(t)] = \phi_i(t)\upsilon_i(t)$, where $\varphi_i(t) = \int_0^1f'[\rho\upsilon_i(t)]d\rho$ and $\phi_i(t) = \int_0^1g'[\rho\upsilon_i(t)]d\rho$. Then, there exist $\varphi^*\in\left( 0, 1\right]$ and $\phi^*\in\left( 0, 1\right]$ such that $\varphi_i(t)\le \varphi^*\le 1$ and $\phi_i(t)\le \phi^*\le 1$, respectively. Hence, $\overset{.}V(t)< 0$ holds true for $n_rb+c< 1$ when $\upsilon\neq0$. Therefore, the trivial equilibrium of the network is globally asymptotically stable when $n_rb+c< 1$ for all time delays.

3.   Case studies

(1) Let $c_1 = -1.4$, $c_2 = 0.4$, $c_3 = 2.5$ and other parameters of the network are the same as the first case in the above section. Based on the obtained analysis, one can conclude that the trivial equilibrium of the network is locally asymptotically stable when $\tau\in\left[ 0, \tau_{1, 0}\right) \cup\left(\tau_{2, 0}, \tau_{1, 1} \right)\cup\cdots\cup\left(\tau_{2, r}, \tau_{1, r+1} \right)\cup\cdots\cup\left(\tau_{2, 7}, \tau_{1, 8} \right)$ and unstable when $\tau\in\left( \tau_{1, 0}, \tau_{2, 0}\right) \cup\left(\tau_{1, 1}, \tau_{2, 1} \right)\cup\cdots\cup\left(\tau_{1, r}, \tau_{2, r} \right)\cup\cdots\cup\left(\tau_{1, 8}, +\infty \right)$, where $\tau_{1, r} = 1.319, 4.656, 7.992, 11.328, \cdots$ and $\tau_{2, r} = 2.456, 6.155, 9.854, 13.553, \cdots$. On the other hand, let $c_1 = 1.5$, $c_2 = 0.8$, and $c_3 = 1.4$. In this case, the trivial equilibrium of the network free of time delays is unstable. After some calculations, the network is locally asymptotically stable when $\tau\in\left( \tau_{2, 0}, \tau_{1, 0}\right) \cup\left(\tau_{2, 1}, \tau_{1, 1} \right)\cup\left(\tau_{2, 2}, \tau_{1, 2} \right)\cup\left(\tau_{2, 3}, \tau_{1, 3} \right)$ and unstable when $\tau\in\left[ 0, \tau_{2, 0}\right) \cup\left(\tau_{1, 0}, \tau_{2, 1} \right)\cup\left(\tau_{1, 1}, \tau_{2, 2} \right)\cup\left(\tau_{1, 2}, \tau_{2, 3} \right)\cup\left(\tau_{1, 3}, +\infty \right)$ where $\tau_{1, r} = 2.702, 5.943, 9.183, 12.424, \cdots$ and $\tau_{2, r} = 0.838, 4.631, 8.242, 12.218, \cdots$. Figure 3 shows the amplitudes of the bifurcated periodic oscillations of the coupled neurons between sub-networks. The results given in Figure 3 reach a very good agreement with the conclusions shown in Figure 1.

(2) $m = 3$, $n_1 = n_2 = n_3 = 3$, $b_{12} = b_{31} = b_{23} = a$, $b_{45} = b_{54} = b_{56} = b_{65} = b_{46} = b_{64} = a$, $b_{78} = b_{87} = b_{79} = b_{97} = b_{89} = b_{98} = b_{77} = b_{88} = b_{99} = a$, $c_k = c$, and $\tau_k = \tau_s$. In this case, the network consists of three-node sub-networks with different topological structures: a uni-directional loop, a bidirectional loop, and a bidirectional loop with self-connections. Figure 4(a) gives the stability region of the trivial solution of the system in the parameter plane of coupling strength and time delay when $a = 0.4$. Let $c = -1.3$. Figure 4(b) shows the region indicating the stability of the trivial solution of the system in the parameter plane of internal connection weight and coupling time delay. Along the vertical line, the characteristic equation of the system has a zero root. Along the curves, the characteristic equation of the system has a pair of imaginary eigenvalues. As shown in Figure 4(b), the vertical line intersects with the curves, which corresponds to the interaction of a steady state bifurcation and Hopf bifurcation.

Figure 5(a) shows that two non-trivial equilibria (0, 0, 0, 0.553, -0.553, 0, 0, 0, 0) and (0, 0, 0, -0.553, 0.553, 0, 0, 0, 0) come into being when $a = -1.1$ and $\tau_s = 1.5$. In this case, the system undergoes a pitchfork bifurcation for $a = -1$. Figure 5(b) gives two coexisting non-trivial fixed points and a branch of periodic oscillations arising from the trivial equilibrium when $a = -1.1$ and $\tau_s = 1.7$. This phenomenon holds true due to the interaction of the pitchfork bifurcation and Hopf bifurcation since there exists a Hopf bifurcation arising from the trivial equilibrium for $\tau_s = 1.66$. Figure 5(c) illustrates the coexistence of two branches of periodic oscillations bifurcating form non-trivial fixed points and a branch of periodic oscillation arising from the trivial equilibrium when $a = -1.1$ and $\tau_s = 2.1$. This fact indicates the system undergoes a Hopf bifurcation arising from the non-trivial equilibrium points for $\tau_s = 2$. As shown in Figure 5, the dashed, dotted and solid curves correspond to bifurcated solutions under different initial conditions.

(3) $m = 3$, $n_1 = n_2 = n_3 = 3$, $b_{11} = b_{44} = b_{77} = 1.18$, $b_{22} = b_{55} = b_{88} = 0.47$, $b_{33} = b_{66} = b_{99} = 1.5$, $b_{31} = b_{64} = b_{97} = 2.9$, $b_{32} = b_{65} = b_{98} = 0.7$, $b_{13} = b_{46} = b_{79} = -2.5$, $b_{12} = b_{45} = b_{78} = 0$, $b_{23} = b_{56} = b_{89} = 2.98$, $b_{21} = b_{54} = b_{87} = -22$, $c_1 = c_2 = c_3 = 0.2$, and $\tau_1 = \tau_2 = \tau_3 = \tau_s$. Figure 6 illustrates a pair of coexisting period-6 oscillations and a pair of chaotic motions when $\tau_s = 0.17$. As shown in Figure 6, the blue, red, green, and purple curves represent four phase trajectories under different initial conditions, respectively. By increasing the time delay $\tau_s = 0.25$, the period-6 orbits disappear and two period-3 solutions come into being as shown in Figure 7. When $\tau_s = 0.43$, Figure 8 gives the coexistence of two chaotic motions and two period-4 responses. With a slightly increase of $\tau_s$, as shown in Figure 9, a pair of period-2 oscillations occur when $\tau_s = 0.47$. Figure 10 gives the coexistence of two pairs of period-4 responses. As shown in Figures 6-10, by changing the values of time delays, the system exhibits rich multi-stability phenomena, such as the coexistence of multi-period oscillations and chaotic attractors and different coexisting multi-period motions.

Figure 11 depicts the detailed bifurcation diagrams on the Poincaré section as the function of the coupling time delays under different initial conditions. 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 in 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é sections are defined by $\sum = \left\lbrace(\tau_s, u_1)\colon(u_3 = 0, \overset{.}u_3> 0) \right\rbrace$. It is shown that the time delays play important roles in the dynamics of the system and can lead to different types of multi-period oscillations, chaotic attractors and multi-stability phenomena. Moreover, there exist period-doubling bifurcations, as depicted in Figure 11.

4.   Conclusions

In biological and physiological systems, multiplex networks often represent neurons in different areas of the brain. In this paper, the stability, bifurcations and complicated oscillations of a multiplex neural network with multiple Hopfield-type neural sub-networks and delayed couplings are studied. Each sub-network contains different sets of nodes and topologies. Time delays are considered in the couplings between sub-networks. The delay-independent and delay-dependent stability of the trivial equilibrium of the multiplex network are discussed by decomposing and analyzing the associated characteristic equation. It is shown that the stability of the network is determined by the root distributions of the characteristic equations of the isolated sub-networks and the product of coupling strengths and the sum of time delays. The regions indicating the stability and instability of the network in the parameter planes are shown. Different types of bifurcations are analyzed, such as steady state bifurcation, Hopf bifurcation and the interaction of them. The coexistence of two non-trivial equilibrium points and periodic oscillations and three coexisting periodic oscillations are obtained. By taking the coupling time delay as parameter, the bifurcation diagrams on Poincaré section are given and interesting dynamical phenomena are shown, such as the coexistence of different patterns of multi-period orbits and strange attractors. The revealed results in this paper can provide promising information for understanding the mechanisms of rhythms, complicated firing patterns, and multi-stability in large-scare complex neural systems.

Acknowledgments

The authors thank the anonymous reviewers for their helpful comments and suggestions that have helped to improve the presentation.