Synchronization in a Multiplex Network of Nonidentical Fractional-Order Neurons

Abstract

Fractional-order neuronal models that include memory effects can describe the rich dynamics of the firing of the neurons. This paper studies synchronization problems in a multiple network of Caputo–Fabrizio type fractional order neurons in which the orders of the derivatives in the layers are different. It is observed that the intralayer synchronization state occurs in weaker intralayer couplings when using nonidentical fractional-order derivatives rather than integer-order or identical fractional orders. Furthermore, the needed interlayer coupling strength for interlayer near synchronization decreases for lower fractional orders. The dynamics of the neurons in nonidentical layers are also considered. It is shown that in lower fractional orders, the neurons’ dynamics change to periodic when the near synchronization state occurs. Moreover, decreasing the derivative order leads to incrementing the frequency of the bursts in the synchronization manifold, which is in contrast to the behavior of the single neuron.

1. Introduction

Fractional calculus has attracted considerable attention due to its wide applications in different fields such as physics and engineering [1,2,3]. Many dynamical systems such as viscoelastic materials can be well expressed by fractional-order differential equations [4]. In the integer-order calculations, the variables are calculated based only on the last value of variables and are therefore memoryless. In contrast, memory is incorporated in fractional calculus by considering a weighted temporal window during integration [5]. This is the most important property of fractional-order equations providing closer models to real systems [6]. This property has led the fractional-order equations to become a proper tool for modeling many dynamical systems, including biological processes [7,8,9,10]. It also can help in understanding the rich dynamics of the neurons’ activities [11]. Many researchers have attempted to propose and analyze fractional-order neuron models [11,12,13,14]. These studies have shown that the neuron model can exhibit various periodic and chaotic firing patterns by varying the fractional order. In addition, the advantages of the Caputo–Fabrizio-type fractional derivatives [15,16] make them very appropriate in modeling real-world problems. The kernel of Caputo–Fabrizio fractional derivative has important features such as nonsingularity, nonlocality and an exponential form. Hence, the newly-defined Caputo–Fabrizio fractional derivative is applied in numerous recent studies [17,18,19,20].

On the other hand, in addition to the pattern of neuronal firing, the study of their collective dynamics is of a great importance [21]. Studies have shown that the interactions among neurons may cause the formation of different collective behaviors such as synchronization, chimera state, spiral waves, etc. [22,23,24]. In fact, the synchronization is an important phenomenon observed in various dynamical systems [25]. For example, the brain’s information processing is accompanied by synchronized activities of the neurons [26]. Fractional-order derivatives can also affect the synchronization behavior of the neurons. Hence, extensive research has been devoted to the investigation of the synchronization problems in neuronal networks [27,28,29], including fractional-order models [30,31,32,33,34]. These studies have revealed many factors affecting the synchronous behavior of neurons. For example, the authors in [35] studied the synchronization of a clustered network model consisting of small-world subnetworks with small heterogeneities. They showed that the global burst synchronization is enhanced by decreasing the heterogeneity. The authors in [36] found that the synchronization between chemically coupled neurons depends on the neurons’ receiving signals rather than the network topology, which is crucial in linearly coupled networks. The paper [37] investigated a neuronal hypernetwork with small-world and random topologies with time-varying links. The investigation for different rewiring frequencies revealed that faster switching causes an enhanced synchronization.

Recent advances in networks theory have led to the utilization of more complex structures [38]. A multilayer framework, which is more suitable for describing many real networks, has achieved a great deal of interest in the last years [39]. This structure allows considering different connections, couplings, or topologies for the nodes in the layers [40]. Hence, it is a proper framework for investigating the behavior of the neurons from different viewpoints. An ability of a multilayer structure is providing a network with nonidentical layers that are closer to real networks. In these networks, two forms of synchronous behavior can be observed. First, the synchronization may occur between the nodes in each layer, while the layers are not necessarily synchronized, which is termed the intralayer synchronization [41]. In contrast, interlayer synchronization occurs when the layers are synchronous [42,43]. However, there exists a gap for corresponding studies on fractional-order models [44].

Motivated by the above considerations, here, we study the interlayer and intralayer synchronizations in a two-layer network of neurons with fractional-order derivatives of Caputo–Fabrizio type. It is assumed that the layers are nonidentical through different derivative orders. We investigate the effect of different derivative orders on the synchronization of neurons. It is shown that the intralayer synchronization is enhanced by considering fractional-order derivatives. The links between the layers with different derivative orders lead to an interlayer near synchronization which is attained in weaker couplings by decreasing the fractional order. Moreover, the interlayer near synchronization manifold is periodic in lower fractional orders. The fractional order Hindmarsh–Rose model and the considered multiplex network are described in Section 2. The synchronous behavior of the network is investigated in Section 3 in detail. The conclusions are given in Section 4.

2. The Fractional-Order Hindmarsh–Rose Model

To introduce the fractional-order Hindmarsh-Rose model that we will investigate, first, some basic fractional calculus definitions and properties will be stated.

Let $q>0$ and $t\ge 0$. An integral of fractional order q for a continuous function f is defined as [2,3]:

$$I^{q}f\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(q\right)}{\int }_0^t\frac{f\left(\tau \right)}{{(t-\tau )}^{1-q}}d\tau ,$$

where $\mathsf{\Gamma }$ is the Gamma function defined by:

$$\mathsf{\Gamma }\left(q\right)={\int }_0^{\infty }{e}^{-t}{t}^{q-1}dt.$$

For $n-1<q<n$, $n\in \mathbb{N}$, the fractional derivative of order q of Caputo type with a lower limit 0 for an absolutely continuous function $f\in C^{n+1}[[0,\infty ),\mathbb{R}]$ is given by:

$${{}_0^{C}D}_t^{q}f\left(t\right)=\frac{1}{\mathsf{\Gamma }(n-q)}{\int }_0^t\frac{f^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{q-n+1}}d\tau .$$

For $q\in (0,1)$, we have:

$${{}_0^{C}D}_t^{q}f\left(t\right)=\frac{1}{\mathsf{\Gamma }(1-q)}{\int }_0^t\frac{f^{\prime }\left(\tau \right)}{{(t-\tau )}^q}d\tau .$$

Note that, for the above fractional-order operators we have [1,2,3]:

$$I^q\left({{}_0^{C}D}_t^{q}f\left(t\right)\right)=f\left(t\right)-f\left(0\right),t\ge 0,q\in (0,1).$$

The physical meaning of the above property is one of the main reasons for the intensive use of the Caputo-type of fractional-order derivatives in the applied problems.

In 2014, the authors in [12] introduced a fractional-order generalization of an integer-order Hindmarsh–Rose model to qualitatively depict the firing behavior of neurons using a group of simple differential equations. The model has been represented by the following fractional-order equations:

$$\left\{\begin{array}{c}{\displaystyle D_t^{q}x=y-a{x}^3+b{x}^2-z+I,}\hfill \\ {\displaystyle D_t^{q}y=c-d{x}^2-y,}\hfill \\ D_t^{q}z=r[s(x-1.56)-z],\hfill \end{array}\right.$$

(1)

where $q\in (0,1)$, the variables $x,y,z$ represent the membrane action potential, a recovery variable and the slow adaption current, respectively. I represents the external current intensity, and the values of the system parameters are considered as $a=1$, $b=3$, $c=1$, $d=5$, $s=4$ and $r=0.006$.

However, the Caputo fractional derivative is not always suitable to accurately describe the memory effect in a real system [1,2,3] because of the singularity of its kernel at the endpoint of an interval of definition. To overcome this difficulty, a new fractional derivative without any singularity in its kernel has been proposed in [15,16]. The Caputo–Fabrizio fractional operator with a kernel in the form of an exponential function for a function $f\in H^1(0,\infty )$, where $H^1(0,\infty )$ is the space of all twice integrable functions with first derivatives that belong to the same space, can be defined as:

$${{}_0^{CF}D}_t^{q}f\left(t\right)=\frac{M\left(q\right)}{1-q}{\int }_0^t{f}^{\prime }\left(\tau \right)exp\left(-q\frac{t-\tau }{1-q}\right)d\tau ,q\in (0,1),$$

(2)

where $M\left(q\right)$ is a normalization function with $M\left(0\right)=M\left(1\right)=1$.

The derivative (2) has been applied by numerous authors in modeling real-world problems to overcome the limitations of real calculus [17,18,19,20,45]. In fact, this definition allows for the description of mechanical properties related to damage, fatigue and material heterogeneities [45].

In our paper, we will generalize the model (1) to consider a two-layer network of Caputo–Fabrizio fractional-order neurons with different fractional-order derivatives in each layer as follows:

$$\left\{\begin{array}{c}{\displaystyle {{}_0^{CF}D}_t^{q_1}{x}_{1,i}=y_{1,i}-x_{1,i}^3+3x_{1,i}^2-z_{1,i}+I_{ext}+\sigma \sum _{j=1}^N{G}_{ij}^1(x_{1,j}-x_{1,i})+ϵ(x_{2,i}-x_{1,i}),}\hfill \\ {\displaystyle {{}_0^{CF}D}_t^{q_1}{y}_{1,i}=1-5x_{1,i}^2-y_{1,i},}\hfill \\ {\displaystyle {{}_0^{CF}D}_t^{q_1}{z}_{1,i}=r[s(x_{1,i}+1.6)-z_{1,i}],}\hfill \\ {\displaystyle {{}_0^{CF}D}_t^{q_2}{x}_{2,i}=y_{2,i}-x_{2,i}^3+3x_{2,i}^2-z_{2,i}+I_{ext}+\sigma \sum _{j=1}^N{G}_{ij}^2(x_{2,j}-x_{2,i})+ϵ(x_{1,i}-x_{2,i}),}\hfill \\ {\displaystyle {{}_0^{CF}D}_t^{q_2}{y}_{2,i}=1-5x_{2,i}^2-y_{2,i},}\hfill \\ {\displaystyle {{}_0^{CF}D}_t^{q_2}{z}_{2,i}=r[s(x_{2,i}+1.6)-z_{2,i}],}\hfill \end{array}\right.$$

(3)

where the layers consist of $N=100$ neurons with small-world connections where each neuron is coupled to its 20 nearest neighbors, and the probability of random links is $0.1$, $\sigma$ and $ϵ$ denote the strength of the intralayer and interlayer links, respectively, $x_{1,i}$, $y_{1,i}$ and $z_{1,i}$ are the variables in the first layer, while $x_{2,i}$, $y_{2,i}$ and $z_{2,i}$ are the variables in the second one, the system parameters are fixed at $s=4$, $r=0.006$, the external current intensity is $I_{ext}=3.2$. The connection between $i$-th and $j$-th nodes is determined by $G_{ij}^1$ and $G_{ij}^2$ in the first and second layers, respectively, such that $G_{ij}^1=1$ if there is a link between $i$-th and $j$-th nodes, and, $G_{ij}^1=0$, otherwise. The schematic diagram of the proposed two-layer network is shown in Figure 1.

Remark 1.

The two-layer network model (3) is an extension of the integer-order models of the Hindmarsh–Rose type (when $q=1$) proposed in [23,24] and some others. It also generalizes the fractional-order model introduced in [12] considering two layers and different Caputo–Fabrizio fractional-order derivatives in each layer.

The modified Adam—Bashforth method is used for the Caputo-Fabrizio fractional operator to solve fractional-order Equation (3). By this method, a fractional-order system of Caputo–Fabrizio typ:

$${{}_0^{CF}D}_t^{q}X\left(t\right)=F(t,X\left(t\right))$$

$$X\left(0\right)=X_0,$$

can be numerically approached as:

$$X(n+1)-X\left(n\right)=\frac{1-q}{M\left(q\right)\mathsf{\Gamma }\left(q\right)}F(t_n,X_n)-F(t_{n-1},X_{n-1})+\frac{q}{1-q}{\int }_{t_n}^{t_{n+1}}F(\tau ,X\left(\tau \right))d\tau$$

(4)

with a solution given by:

$$X(n+1)-X\left(n\right)=\left(\frac{1-q}{M\left(q\right)\mathsf{\Gamma }\left(q\right)}+\frac{3q\mathsf{\Delta }t}{2M\left(q\right)}\right)F(t_n,X_n)-\left(\frac{1-q}{M\left(q\right)\mathsf{\Gamma }\left(q\right)}+\frac{q\mathsf{\Delta }t}{2M\left(q\right)}\right)F(t_{n-1},X_{n-1}),$$

(5)

where $\mathsf{\Delta }t$ is the step size of time integration.

The Hindmarsh–Rose model (3) represents a variety of dynamics by changing the order of derivatives. The time series of the single Hindmarsh–Rose model for $q=1,0.9,0.8$ and $0.7$ are shown in Figure 2, for fixed parameters $r=0.006,s=4,I_{ext}=3.2$. It can be observed that for $q=1$, the behaviors of the states are chaotic square-wave bursting. By decreasing the derivative order from 1 to $q=0.9$, the waveform is changed from square-wave bursting to triangular-wave bursting. Furthermore, the maximum value of the spikes decreases for the fractional-order derivatives. When the order decreases to $0.8$, the number of spikes in bursts increases, while the frequency of the bursts decreases. For $q=0.7$, the waveform changes to a pseudo-plateau bursting.

3. Behavior of Network

In this section, the behavior of the states of the fractional two-layer network model (3) is investigated for different derivative orders, and the synchronization of the network is studied. To this aim, a numerical solution is established, and the interlayer synchronization and the intralayer synchronization errors are computed as follows:

$$\left\{\begin{array}{c}{\displaystyle E_1=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^T\sum _{j=2}^N\frac{\left|\right|X_{1,j}\left(t\right)-X_{1,1}\left(t\right)\left|\right|}{N-1}dt,}\hfill \\ {\displaystyle E_2=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^T\sum _{j=2}^N\frac{\left|\right|X_{2,j}\left(t\right)-X_{2,1}\left(t\right)\left|\right|}{N-1}dt,}\hfill \\ {\displaystyle E=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^T\sum _{j=1}^N\frac{\left|\right|X_{1,j}\left(t\right)-X_{2,j}\left(t\right)\left|\right|}{N}dt,}\hfill \end{array}\right.$$

(6)

where $X=(x,y,z)$ is the state vector of variables for each layer, $E_1$ and $E_2$ are the intralayer synchronization errors for the first and second layers, respectively, and E is the interlayer synchronization error.

The intralayer and interlayer synchronization errors for different derivatives are represented in Figure 3 and Figure 4, respectively. Since the connections within layers are approximately similar, there is a negligible difference between the intralayer synchronization errors ($E_1\equiv E_2$). Thus, only the error for the first layer ($E_1$) is represented. Figure 3 shows that intralayer synchronization is generally achieved in lower intralayer coupling strengths by increasing the strength of links between layers. Moreover, as the derivative orders decrease, the intralayer synchronization is mainly obtained for weaker intralayer couplings. However, an unexpected change occurs around $ϵ=1$. In almost all cases, the intralayer synchronization error in the area around $ϵ=1$ is less than the error in its surrounding areas.

The intralayer synchronization error for $ϵ=1$ according to intralayer coupling strength is shown in Figure 5a. This figure shows that for identical integer-order derivatives, the synchronization is obtained for $\sigma =0.3$. By changing the order of derivatives to identical fractional orders, the threshold of coupling strength for synchronization is slightly decreased. Interestingly, when the fractional orders of layers are different, this threshold is considerably varied. When $q_1$ decreases to $0.9$, the threshold decreases to $\sigma =0.18$, and for $q_1=0.8$, it is equal to $\sigma =0.9$. With more decreasing of $q_1$, this threshold increases again such that for $q_1=0.7$, the synchronization is obtained for $\sigma =0.11$. For $q_1=0.8$, $q_2=0.9$ and $q_1=0.7$, $q_1=0.8$, the threshold changes to $\sigma =0.18$ again. Therefore, the intralayer synchronization in nonidentical fractional orders is obtained in weaker couplings than the identical orders. The one-dimensional interlayer error is shown in Figure 5b and represents the difference between the results of different derivative orders.

When the order of the derivatives in both layers is the same, the neurons are identical and can reach a complete synchronization state with zero interlayer synchronization error. However, when the order of the derivatives is different, the neurons are nonidentical. Therefore, complete interlayer synchronization cannot be achieved. Instead, the neurons can reach a near synchronization state with a minimum interlayer synchronization error. Figure 4a shows that in the integer-order network, complete synchronization occurs in about $ϵ>0.5$. For fractional-order derivatives, the minimum error for near synchronization is obtained in significantly larger interlayer coupling strengths (Figure 4b–f). A comparison analysis of these figures leads to the conclusion that settings $q_1=0.8$, $q_2=0.9$ and $q_1=0.7$, $q_2=0.8$ imply near synchronization in weaker interlayer couplings than others. Moreover for $q_1=0.7,q_2=0.8$, the neurons become unstable for special coupling strengths, the region of which is shown in black.

As shown in Figure 2, different fractional orders lead to different firing patterns and attractors for the neuron model. When nonidentical neurons are coupled, the coupling strength between them brings their manifolds closer together to eventually reach a near synchronization. We represent the effect of coupling strength on the attractor of the neurons in bifurcation diagrams in Figure 6. In these figures, the maximum values of voltages of the first neuron of each layer are plotted with different colors.

It is observed that only for $q_1=1,q_2=0.9$ and $q_1=0.9,q_2=0.8$, the chaotic firing is preserved by increasing the coupling strength, and the neurons have a chaotic near synchronization state. In other cases, however, strengthening the coupling leads to the bifurcation of the neurons to periodic behavior, and the greater the difference between the two derivative orders, the lower the coupling strength for this bifurcation. In other words, when the two derivative orders have more differences, the chaotic region becomes smaller. On the other hand, the dynamics of the neuron with lower derivative order have more of a tendency to change to periodic. Moreover, it is observed that although the coupling tries to synchronize the neurons, the trajectories do not precisely overlap, and there is a minimal error that does not disappear with increasing coupling coefficient.

The time series and attractors of the first neurons of each layer are represented in Figure 7 in the synchronous manifold. It can be observed that although having a near synchronization state, the attractors of the neurons of each layer are different. Furthermore, as the derivative orders decrease, the frequency of the bursts increases. The highest burst frequency value is in $q_1=0.8,q_2=0.7$ (Figure 7e) case, and the lowest frequency is obtained for $q_1=1,q_2=0.9$ case (Figure 7a). Comparing Figure 2 and Figure 7 we conclude that this result contrasts the dynamics of the single neuron wherein the frequency of burst decreases by decreasing the fractional order.

4. Conclusions

In this paper, we proposed a two-layer multiplex network of neurons considering different Caputo–Fabrizio-type fractional derivative orders in each layer. The synchronization of the neurons within and between the layers was studied by computing the average synchronization error for different derivative orders. It was observed that the nonidentical fractional-order derivative leads to the enhancement of the intralayer synchronization (i.e., the synchronization is achieved in lower intralayer coupling strengths). Moreover, as the derivative order decreases, the needed intralayer coupling strength is decreased. Since the neurons of different layers are nonidentical through different derivative orders, complete interlayer synchronization cannot be obtained. Instead, a near interlayer synchronization can be achieved in which the interlayer synchronization error does not decrease to a certain extent by increasing the coupling coefficient. The results showed that near interlayer synchronization occurred in weaker strengths for lower fractional orders. The dynamics of both layers’ neurons were also considered by varying the interlayer coupling strength. It was observed that when the derivative orders are different, increasing the coupling strength leads to bifurcation in neurons dynamics. Furthermore, the neurons have a periodic near synchronization manifold for lower derivative orders. In addition, the frequency of the bursts increases by decreasing the derivative order, while in the single neuron, the frequency of the bursts in lower derivatives is decreased. An interesting and important future direction of our research is related to the study of the presented model with the help of fractal-fractional derivatives. The established experimental data support as well as molecular level explanations will be very helpful in such direction.