Dynamical Analysis and Synchronization of a New Memristive Chialvo Neuron Model

Abstract

Chialvo is one of the two-dimensional map-based neural models. In this paper, a memristor is added to this model to consider the electromagnetic induction’s effects. The memristor is defined based on a hyperbolic tangent function. The dynamical variations are analyzed by obtaining the bifurcation diagrams and Lyapunov spectra. It is shown that the most effective parameters on the dynamics are the magnetic strength and the injected current. The memristive Chialvo can exhibit different neural behaviors. It is also proven that, like the primary Chialvo model, the memristive version has coexisting attractors; an oscillating state coexists with a fixed point. In addition, to understand how memristive neurons behave in a network, two memristive Chialvo models are coupled with electrochemical synapses. By connecting two neurons and calculating the synchronization error, we can determine the system’s synchronizability. It is indicated that the electrical coupling is essential for the occurrence of complete synchronization in the network of memristive Chialvo, and the sole chemical coupling does not lead to synchronization.

1. Introduction

Introducing mathematical models for neurons facilitates computational analysis of brain behavior. Computational neuroscience has evolved significantly since Hodgkin and Huxley determined the contribution of sodium and potassium ions in a neuron’s axon and investigated its action potential in 1952 [1]. After the 4D Hodgkin and Huxley model, various mathematical models have also been introduced, and their electrical circuits have been implemented [2,3]. The 3D Hindmarsh–Rose [4], the 2D Morris–Lecar [5], and the 2D FitzHugh–Nagumo [6] models are examples of them in the continuous time domain. There are also other neural models in the discrete-time domain, such as the Rulkov map [7], the Zandi neuron map [8], the Cazelles–Courbage–Rabinovich map [9], and the Chialvo map [10] models. Not only are map-based models fast, but they are also easy to use. In a continuous model with dimensions lower than three, the chaotic oscillation of the neuron cannot be modeled. In contrast, discrete models do not follow a similar rule to illustrate chaotic oscillations. Hence, the networks of continuous neural models are time-consuming, especially when coupled in a high-dimensional network or higher-order one [11,12]. Even though some reduction methods reduce the cost of analyzing a high-dimensional network of oscillations [13], map-based models are generally faster and more efficient for analyzing large networks.

By studying neurons in a network of oscillations, one can analyze different types of synchronization, including partial synchronization [14], phase synchronization [15,16], lag synchronization [17], and complete synchronization [18,19]. A complete synchronization is a state where all the neurons in the network oscillate at the same dynamics simultaneously. Synchronization plays an essential role in the brain’s general functioning, including working memory and motor planning [20,21,22,23]. A variety of types of synchronization have also been observed in patients with epileptic seizures [24], Parkinson’s disease [25], schizophrenia [26], and other brain disorders.

The brain transfers information between neurons via synapses, where neural interactions occur. As the axon fires, the presynaptic neuron releases neurotransmitters in the synaptic cleft (the gap between two synapses). Each neurotransmitter can activate its receptors in the postsynaptic neuron. Then, based on involved ions, the postsynaptic neuron can be either hyperpolarized or depolarized. These neurotransmitters can be electrical or chemical. Also, in some cases, both neurotransmitters can exist simultaneously in the same synapse. As well as chemical and electrical synapses, an axon can receive feedback from itself via autapses [27]. Autapses are self-synapses in the neuron that are common throughout the brain [28], both excitatory and inhibitory [29,30].

Neurons can also be affected by external stimuli, such as temperature [31], noise [32], or electromagnetic changes [33]. The neuron model should be added with a memristor to study the effect of electrical changes on magnetic flux. Prof. Chua first introduced the concept of a memristor or a memory resistor as a fourth fundamental electrical component in addition to resistance, inductance, and capacitance [34]. Meanwhile, the field of studying memristor-based systems became very popular after Strukov et al. implemented it physically in 2008 [35]. By introducing continuous and discrete memristors, many advances were accomplished in chaotic systems [36,37,38] and nonlinear electronic circuits [39,40,41].

Because of being nonlinear and having hysteresis loops, memristors are applicable in various scientific fields. Memristor-based neural networks and analog memristor arrays have lately been employed in machine learning or pattern recognition approaches [42,43,44] and also in image encryption [45]. Numerous studies have also been conducted on the memristive model of both continuous and discrete mathematical neural models. Memristive Hindmarsh–Rose [46,47,48], memristive Hodgkin–Huxley [49], memristive Rulkov [50,51,52], and memristive Zandi [53] are some examples. The neuronal models can be used to simulate neuronal populations and reproduce their behaviors, which are relevant to many brain disorders and cognitive tasks. However, most proposed models are flow-based, resulting in a high computational cost. The map-based models can decrease this cost considerably and are a good substitute for flow-based models, especially when investigating neuronal networks. Therefore, improving the map-based models by adding the memristor can be helpful for this purpose. Hence, this paper investigates the effect of an external electromagnetic field on the Chialvo map-based model.

Here, we consider the hyperbolic tangent function as the memductance function. The electromagnetic effects are mainly examined through the bifurcation diagram and Lyapunov exponents for different parameters. It is shown that magnetic strength significantly influences the neuron’s dynamics. The effects of the initial condition are also examined, and coexisting attractors are shown. Moreover, by coupling two memristive Chialvo models, the synchronous and asynchronous behaviors of the model are analyzed by computing the synchronization error. Both electrical and chemical synapses are investigated to see how each coupling affects its synchronizability. The dependence of the coupled neurons’ dynamics on their initial conditions is also revealed.

The dynamics of the memristive Chialvo are studied in Section 2, consisting of the system’s fixed points and stability. The effect of different system parameters (especially the magnetic strength) on the system dynamics is also investigated through numerically simulating bifurcation diagrams and Lyapunov exponents, the neuron time series, and phase space. The effects of initial conditions in the system are also considered through the basin of attraction diagram. In Section 3, the synchronizability of the system is examined when two memristive neurons are chemically and electrically coupled. Computing synchronization errors shows that chemical coupling strength is essential for complete synchronization to occur. Finally, the conclusion of the paper is provided in Section 4.

2. Neuron Dynamics

Take n as the time step, the 2D map model of Chialvo is described as follows [10]:

$$\{\begin{array}{c}x\left(n+1\right)=x^2\left(n\right)\exp \left(y\left(n\right)-x\left(n\right)\right)+I\\ y\left(n+1\right)=ay\left(n\right)-bx\left(n\right)+c \end{array}$$

(1)

where $x$ represents the activation variable and $y$ represents the recovery variable. The four remaining parameters, i.e., $I$, $a$, $b$, and $c$, are, respectively, the injected ion current, the time of recovery ($a<1$), the activation dependence of the recovery process ($b<1$), and the offset. The model can illustrate both chaotic and periodic behaviors of the neuron. Here, to analyze the electromagnetic effects on this map-based neuron, the $x$ parameter in Equation (1) is added with a discrete memristor. A memristor with hyperbolic tangent flux controls, which is bounded, smooth, and nonlinear, is described in Equation (2) [50].

$$\begin{array}{c}F\left(n\right)=\tan \mathrm{h}\left(\varphi \left(n\right)\right)x\left(n\right)\\ \varphi \left(n+1\right)=\varphi \left(n\right)+ϵx\left(n\right)\end{array}$$

(2)

In Equation (2), $\varphi$ is the flux variable, $ϵ$ is the time scale factor for generating the induced electromotive force in finite transients, and $x$ is the voltage of the memristor. Accordingly, the modified 3D memristive Chialvo model (3D-MCM) can be formulated as follows:

$$\{\begin{array}{c}x\left(n+1\right)=x^2\left(n\right)\exp \left(y\left(n\right)-x\left(n\right)\right)+I+k\tan \mathrm{h}\left(\varphi \left(n\right)\right)x\left(n\right)\\ y\left(n+1\right)=ay\left(n\right)-bx\left(n\right)+c \\ \varphi \left(n+1\right)=r\varphi \left(n\right)+ϵx\left(n\right)\end{array}$$

(3)

in which $k$ is the magnetic strength parameter and $r$ is the scale factor of flux.

The fixed point of the 3D-MCM is the point that all state variables ($\mathit{X}=\left\{x,y,\varphi \right\}$) satisfy $\mathit{X}\left(n+1\right)=\mathit{X}\left(n\right)={\mathit{X}}^{\ast }$, which leads to the following:

$$\{\begin{array}{c}{x}^{\ast }=x^{\ast }{}^2\exp \left(y^{\ast }-x^{\ast }\right)+I+k\tan \mathrm{h}\left({\varphi }^{\ast }\right)x^{\ast }\\ y^{\ast }=a{y}^{\ast }-b{x}^{\ast }+c \\ {\varphi }^{\ast }=r{\varphi }^{\ast }+ϵx^{\ast }\end{array}$$

(4)

Considering the parameters in

Table 1. The value of parameters for 3D-MCM.

$\mathit{a}$	$\mathit{b}$	$\mathit{c}$	$\mathit{k}$	$\mathit{ϵ}$	$\mathit{r}$	$\mathit{I}$
$0.89$	$0.18$	$0.28$	$0.145$	$1$	$0.95$	$0.005$

, the fixed point obtained by solving Equation (4) is $S = \left(x^{\ast }, y^{\ast }, {\varphi }^{\ast }\right) = \left(0.005, 2.536, 0.109\right)$. To analyze the stability of the model equilibria, the Jacobian of the 3D-MCM is calculated as follows:

$$J_s=\left[\begin{array}{ccc}k\tan \mathrm{h}\left(\varphi \right)-x^2\exp \left(y-x\right)+2x\exp \left(y-x\right)& x^2\exp \left(y-x\right)& -kx(\tan \mathrm{h}{\left(\varphi \right)}^2-1)\\ -b& a& 0\\ 1& 0& r\end{array}\right]$$

(5)

The fixed point of the 3D-MCM is stable when all Jacobean eigenvalues at it are less than unity. The eigenvalues of $J_s$ at the fixed point ($S$) are $L1=0.1403$, $L2=0.9509$, $L3=0.8899$ showing that $S$ is stable. The stability analysis of other fixed points can be conducted similarly.

To demonstrate the dynamics of the 3D-MCM as a result of the variation of four system parameters ($k$, $ϵ$, $r$, and $I$), the bifurcation diagrams (Figure 1 and Figure 2) and Lyapunov exponents (Figure 3) are obtained. As shown in Figure 1, two initial conditions $\left(1.0, 0.8, 0.2\right)$ and $\left(0.5, 0.2, 0.3\right)$ represent the bursting (purple) and quiescent (blue) state of the system, respectively, indicating that the 3D-MCM has coexisting attractors when the same parameter and different initial conditions are chosen. The straight blue line in the bifurcation diagrams shows the stable fixed-point attractor $s=\left(0.00, 2.54, 0.109\right)$, which does not bifurcate by changing the parameters. Chaotic behavior and several periodic windows can be observed at $0.1417<k<0.1682$ in Figure 1a, where the largest periodic window is approximately $0.1458<k<0.1480$. For $k>0.1682$, the system has two different fixed points, while for $k < 0.1417$, both initial conditions lead to the single fixed point $s=\left(0.00, 2.54, 0.109\right)$. In Figure 1b for $ϵ<0.096$, in Figure 1c for $0.7103<r$, and in Figure 1d for $0.004069<I<0.01438$, chaotic behavior with periodic windows can be seen. Both initial conditions lead to one fixed point in Figure 1c for $0.7103>r$ and Figure 1b for $ϵ<0.096$. Additionally, in Figure 1d for $I< 0.004069$ and $I>0.01781$, a single fixed point is obtained, whereas coexisting fixed points are evident for $0.01438<I<0.01781$, which both satisfy Equation (3), based on the parameters. To have better observation, the bifurcation diagrams are plotted with respect to two bifurcation parameters, and the system’s dynamics are plotted with particular colors in Figure 2 with (1.0, 0.8, 0.2) as the initial condition. The purple color shows the chaotic oscillation, the red one is the periodic oscillation, and the blue regions show the quiescent mode of the 3D-MCM. Figure 2a depicts the system bifurcation in the I-k plane in which some periodic oscillation (red region) can be observed in the chaotic zone (purple regions), while there is no quiescent mode (blue region) in the chaotic zone. In Figure 2b, which illustrates the system bifurcation in the r-ϵ plane, chaotic and quiescent zones are somehow separated by periodic regions.

The Lyapunov exponents corresponding to each panel of Figure 1 are presented in Figure 3. Using $\left(1.0, 0.8, 0.2\right)$ as the initial condition, chaotic ranges with at least one positive LE indicate aperiodic behavior in the mentioned ranges.

Figure 4 illustrates the periodic and chaotic dynamics of the 3D-MCM for the parameters in Table 1 and two different values for $k$ in the time domain and phase space. Figure 4a shows the periodic orbits of the 3D-MCM at $k=0.148$, and Figure 4b shows the chaotic dynamics at $k=0.152$ with the initial condition $\left(1.0, 0.8, 0.2\right)$. Throughout the figures, the transient parts of the 3D-MCM are represented by the bright gray points, while dark blue ones indicate the steady states.

By analyzing 3D-MCM dynamics, it has been shown that the model demonstrates a wide range of neuronal behaviors, including quiescence, periodic oscillation, and chaotic oscillation. The map-based behaviors are more efficient than those of other neuron models, especially in large-scale networks. As the 3D-MCM is a memristor-based model with memory, it can be used in learning approaches or other fields needing memory. The proposed 3D-MCM is compared to some other memristor-based and non-memristor-based neuron models in

Table 2. Comparison of the proposed model with some other memristive neuron models.

Memristive Model	Memductance Function	Discrete	Summary	Year
Morris–Lecar with memristive autapse [54]	Polynomial	No	The dynamics of the proposed Morris–Lecar neuron model with memristive autapse are investigated theoretically with the help of time series, phase space, and the bifurcation diagram. It demonstrates that memristive autapse can cause multistability. Also, this memristor can enhance the synchronizability of the coupled memristive autaptic neurons.	2023
Memristive Rulkov neuron [55]	Sinusoid	Yes	The HP memristor adds the Rulkov neuron to analyze the memristor-based neuron. Phase portraits, bifurcation structures, and spectral entropy complexity are then used to study its dynamics. It shows that the resistance of the memristor can postpone the bifurcation point. Additionally, the memristive model has a wider range of complexity.	2022
Memristive Hindmarsh–Rose [56]	Time-delayed	No	The effect of a time-delay memristor is examined by using the 2D Hindmarsh–Rose model. When the external source is a DC source, it has been demonstrated that the memristor may produce chaotic or hyperchaotic currents by adjusting the time-delay parameter. The bifurcation diagram, Lyapunov exponents, and phase portraits are used to study the dynamical behaviors of the suggested memristor model.	2021
Memristive improved Izhikevich [57]	Polynomial	No	The memristive version of the improved Izhikevich model is analyzed to see how the external magnetic field affects the neuron’s firing pattern. Model dynamics are only examined using time series and bifurcation diagrams.	2020
3D Memristive Chialvo Model	Hyperbolic tangent	Yes	In this paper, the 3D memristive Chialvo model is proposed to study the impact of hyperbolic tangent fluctuation on neuron dynamics. The stability of the system’s fixed point is first analyzed. The variation of the model based on system parameters, particularly the magnetic strength, is then examined using single and double-parameter bifurcation diagrams and corresponding Lyapunov exponents, in addition to time series analysis.	-

. Most studies have demonstrated that memristor resistance has a significant impact on system dynamics.

3. Synchronizability of the 3D-MCM

An analysis of synchronizability can be done by coupling at least two neurons. Hence, the effect of flux $\varphi$ on the collective behavior of 3D-MCM synchronization is investigated by connecting a pair of neurons with electrochemical synapses, as follows:

$$\{\begin{array}{c}{x}_i\left(n+1\right)=x_i^2\left(n\right)\exp (y_i\left(n\right)-x_i\left(n\right)+I+k\tanh \left({\varphi }_i\right)x_i\left(n\right)\\ +g_{el} \left(x_{i^{\prime }}\left(n\right)-x_i\left(n\right)\right)\\ +g_{ch}\left(V_s^{\ast }-x_i\left(n\right)\right)\left(\frac{1}{1+\exp \left(-\beta \left(x_{i^{\prime }}\left(n\right)-{\theta }_s\right)\right)}\right)\\ y_i\left(n+1\right)=a{y}_i\left(n\right)-b{x}_i\left(n\right)+c\\ {\varphi }_i\left(n+1\right)=r{\varphi }_i\left(n\right)+ϵx_i\left(n\right)\end{array}$$

(6)

In Equation (6), which shows the two coupled 3D-MCM, $x_i$ is the states of $i$-th neuron and $x_{i^{\prime }}$ denotes the other neuron states, where $i=1, 2$. Also, $g_{el}$ and $g_{ch}$ are, respectively, electrical, and chemical coupling. Here, the coupled model parameters are mentioned in

Table 3. The value of parameters for two coupled 3D-MCM.

$\mathit{a}$	$\mathit{b}$	$\mathit{c}$	$\mathit{k}$	$\mathit{ϵ}$	$\mathit{r}$	$\mathit{I}$	${\mathit{V}}_{\mathit{s}}^{\ast }$	$\mathit{\beta }$	${\mathit{\theta }}_{\mathit{s}}$
$0.89$	$0.18$	$0.28$	$0.145$	$1$	$0.95$	$0.005$	$-1.4$	$50$	$-1.4$

.

The synchronizability of 3D-MCM is determined by the synchronization error, which is defined as:

$$\mathrm{error}={\langle \sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left({\varphi }_1-{\varphi }_2\right)}^2}\rangle }_n$$

(7)

in which ${\langle .\rangle }_n$ indicates the mean over samples. The effects of chemical and electrical coupling strengths on synchronizability are examined through error. First, the error is calculated concerning only one coupling to investigate how each coupling strength affects the coupled neurons’ synchronizability. Then, the effect of both chemical and electrical coupling strengths on its synchronizability is analyzed.

Figure 6a shows the error in the absence of chemical coupling ($g_{ch}=0$) and varying $g_{el}$, while Figure 6b shows the error for the opposite condition ($g_{el}=0$ and varying $g_{ch}$) when the initial conditions are chosen as $\left(x_1\left(0\right),y_1\left(0\right),{\varphi }_1\left(0\right)\right)=\left(1.0,0.8,0.2\right)$ and $\left(x_2\left(0\right),y_2\left(0\right),{\varphi }_2\left(0\right)\right)=\left(0.5,0.2,0.3\right)$. In Figure 6a, error $=0$ indicates that the neurons are synchronous when $g_{el}>0.0436$, while in Figure 6b, the error is always positive, which shows that the neurons cannot be synchronized in the absence of $g_{el}$.

The synchronization error of the network is also calculated when both chemical and electrical couplings vary. Black regions in Figure 7a indicate zero error, hence the complete synchronization of two neurons. As $g_{el}$ increases, the synchronization error gets lower. It confirms that the network cannot be fully synchronized without the $g_{el}$ parameter. Figure 7b indicates the error in the $k$-$g_{el}$ plane when $g_{ch}=0.00012$. According to Figure 1a, in this diagram, $k$ is chosen in the chaotic range of a single 3D-MCM. The black spots in the error plane, where error $=0$, illustrate that by adjusting chemical coupling strength, two neurons can be synchronized whether the single neuron is in its quiescent mode or oscillates chaotically or periodically. The error in Figure 7a indicates that if electrochemical coupling in the synapses is considered, for larger $g_{ch}$, larger $g_{el}$ is required for full synchronization. In addition, in Figure 7b, it is observed that as $k$ increases, the dynamics of the system become more ordered and smaller $g_{el}$ is needed to achieve zero synchronization error.

Figure 8 shows four time series of two MCMs with different coupling strengths. In Figure 8a, the chemical coupling is absent, but complete synchronization has occurred; however, in Figure 8b, when the electrical coupling is absent, synchrony in neurons does not occur; one neuron oscillates chaotically, and the other is placed in the quiescent state. In addition, two different patterns of asynchronous oscillation are illustrated in Figure 8c,d.

The network basin of attraction is also plotted for coupled 3D-MCMs. Figure 9 shows the oscillators’ dynamics with initial conditions $\left(x_1\left(0\right),0.8,0.2\right)$ and $\left(x_2\left(0\right),0.2,0.3\right)$. The initial conditions in the purple region are $x_1\left(0\right)$ and $x_2\left(0\right),$ which leads the networks to oscillate chaotically but not necessarily synchronously. The ones in the blue region lead to the quiescent state. The system parameters are set as mentioned in Table 3.

4. Discussion and Conclusions

The Chialvo neuron map model under electromagnetic fluctuation was studied, as described in this paper. A hyperbolic tangent version of the memristor was added to the activation variable of the 2D Chialvo map, and its dynamics were analyzed. The fixed point of the memristive Chialvo was calculated, and its stability was demonstrated. Moreover, the effects of model parameters were examined, revealing that the magnetic strength can mainly cause the model to bifurcate. Also, the occurrence of chaotic behavior by adjusting parameters was evident through plotting the bifurcation diagram and calculating Lyapunov exponents. Additionally, the memristive model’s time series and phase space were plotted to demonstrate the periodic and chaotic patterns of the neuron. Based on the initial condition, the system has coexisting attractors (bursting or quiescent mode). The synchronizability of the memristive Chialvo model was also examined by coupling two identical neurons with electrochemical synapses. The computation of synchronization error showed that in the absence of electrical coupling, complete synchronization cannot arise; however, the electrical coupling can individually synchronize the two coupled neurons without chemical coupling. In the end, the effect of different coupling strengths on the neurons’ dynamics was illustrated in the time domain.