Dynamic analysis of the discrete fractional-order Rulkov neuron map

1.   Introduction

Computational neuroscience studies the function of the brain from different molecular, cellular and behavioral levels ^[1]. The mathematical neuronal models are used to evaluate the qualitative function of neurons ^[2,3]. A wide variety of neural activities, such as different firing patterns, can be represented by neuronal models ^[4,5,6,7,8]. Neuronal models are generally divided into continuous time and discrete time. Continuous time neuronal models are described by the ordinary differential equations (ODE) such as Hodgkin-Huxley ^[9], Fitzhugh-Nagumo ^[10,11], and Hindmarsh-Rose ^[12,13]. For example, Njitack and colleagues examined the effect of a memristive autapse on the dynamics of a two-dimensional Hindmarsh–Rose neuron ^[14]. In another work, the effect of a two-dimensional FitzHugh–Nagumo, along with a 3D Hindmarsh–Rose neuronal model, was evaluated by coupling a multistability memristive synapse ^[15]. Models of continuous time neurons have been extensively studied dynamically, while less attention has been paid to the map-based discrete time models ^[16]. The distinctive features of discrete-time neuron maps are their high potential for complex behavior modeling, memory savings, and simplifying large neural network calculations ^[17,18]. Therefore, many scientists have tried to introduce discrete-time neuron maps ^[19]. For example, with the Euler discretization of the Hodgkin-Huxley neuronal model, the Izhikevich neuronal model was introduced ^[20]. The Nagumo-Sato discrete neuronal model was proposed in 1972 ^[21], and in 1990 a modified version of this model was introduced^[22]. A one-dimensional map-based neuronal model of self-interaction was proposed in 1997 by Pasemann ^[23]. In 2001, Rulkov proposed a simple two-dimensional turbulent neuron model, one of the first phenomenological models to show a variety of basic neural behaviors (bursts of spikes, spikes, silence) ^[24]. Mathematical modeling of bursts and some spikes is impossible in the context of two-dimensional continuous-time systems, but Rulkov's two-dimensional model can demonstrate these behaviors ^[25].

The dynamic behaviors of the coupled Rulkov neuron map have been studied in recent research ^[26,27,28]. Sun et al. investigated the complete synchronization of Rulkov chaotic neural networks using master stability function analysis ^[29]. Franović investigated the effects of synaptic time delay on Rulkov map neurons with chemical coupling ^[30]. Rakshit et al. examined the stabilization, neural synchronization, and dynamics of two coupled Rulkov neurons with electrical and chemical synaptic interactions ^[31]. Wagemakers et al. introduced an electronic implementation of the Rulkov neuron map ^[32]. Li et al. proposed a discrete memristive neuron model based on the Rulkov model, which can better describe the actual firing activity of biological neurons due to biophysical memory ^[33]. Cheng et al. studied the synchronization dynamics of two heterogeneous chaotic Rulkov neurons with electrical synapses ^[34]. It is notable that the target for synchronization stability means energy balance between neurons and the neurons reach energy balance after complete synchronization ^[35,36].

Discrete fractional calculus has been used in the dynamic analysis of map-based systems. The fractional-order derivatives are relatively more accurate than actual order derivatives because they are an effective tool for describing the effect of memory on processing ^[37]. Memory effects in fractional order indicate that all previous models are involved in determining the system modes ^[38]. Over the past several years, researchers have focused on incorporating the fractions in various applications of neural networks, physics, biology, etc. ^{[39,40,41,42,43]}. The discrete fraction order equation is also a powerful tool in modeling chaotic systems ^[44,45,46]. For example, Guo-Cheng et al. reported the delayed logistic equation using the discrete fractional calculus approach and the corresponding discrete perturbation ^[47]. Elsonbaty et al. proposed a susceptible-infected-treatment-recovered-susceptible (SITRS) fractional discrete model to simulate the coronavirus (COVID-19) epidemic ^[48]. Kassim et al. developed a secure image transfer scheme based on synchronizing turbulent time-discrete fractional systems ^[49]. Fractional derivative properties can be expressed in the biological sense as follows ^[50]:

1. The fractional-order derivative represents the inheritance characteristics more accurately compared to integer-order models.

2. FO neuron models correctly represent biological properties in the presence of noise, which full-order models fail to do.

3. Dynamic firing activities of multiple time scales for a single neuron can be treated with FO derivatives.

4. Firing frequency responses are better in fractional designs than in integer models.

5. Dynamic excitatory properties reveal several neurocomputational properties that can be achieved in FO systems, which enrich functional neural mechanisms.

In general, human evolution is done by two genetic systems based on DNA and another based on the transmission of information through the functions of the nervous system. Mechanisms underlying neuron-based inheritance include hippocampal neurogenesis and memory and learning processes that generate new neural complexes and alter brain structure and function ^[51]. From the concept of neuroscience, the discrete fractional order models the neuron's hereditary behaviors as memory in the dynamic model. Given the above, this paper aims to assist researchers in computational neuroscience by introducing a discrete fraction order derivative to the Rulkov neuronal model. The dynamics of the fractional order map is analyzed completely and, in the end, its synchronization is examined. Synchronization of neurons has become an interesting topic in recent years ^[26,52,53] due to its critical applications in the neural processes ^[54,55]. This paper includes the following sections: Section 2 introduces the Rulkov neuron map, demonstrating the system's dynamic and biological behaviors. Section 3 presents a discrete fractional order Rulkov neuron map and dynamically evaluates the effect of fractional order and model parameters. In section 4, the fixed points are calculated, and their stability is examined. Section 5 investigates the synchronization of two fractional-order Rulkov neurons in different fractional orders, and finally, Section 6 gives conclusions.

2.   Mathematical Rulkov neuron map

In this section, Rulkov neuron map is introduced as a discrete system. Rulkov neuron map is a discrete biological neuronal model proposed by Nikolai F. Rulkov in 2001 ^[24,25]. This model, in which $n$ indicates the discrete time, is defined as follows:

where the variable ${x}_{n}$ is the membrane potential of neurons, and the variable ${y}_{n}$ in the model is a system-slowing variable due to the minimal value of $\mu$ ($0 < \mu < < 1$). Unlike the variable ${x}_{n}$, the variable ${y}_{n}$ has no clear biological meaning ^[30]. The parameter $\sigma$ is an external dc current to neurons, and the parameter $\alpha$ is a nonlinear map parameter. The value of $\mu$ is always considered $0.001$. By changing the parameters $\sigma$ and $\alpha$ of the system, different behaviors such as silence ($\sigma = -2 , \alpha = 6.22$), bursting ($\sigma = -1 , \alpha = 4.1$) and chaotic firing ($\sigma = -1 , \alpha = 5.7$) can be seen which are shown in Figure 1 for random initial conditions. Phase diagram for burst mode and chaotic firing mode are shown in parts (d) and (e) of Figure 1.

In continuation, the dynamic analysis of the Rulkov neuron map is performed through bifurcation diagrams and the maximum Lyapunov exponents (MLE) calculations. The positive MLE proves the existence of chaotic behavior. The bifurcation diagram is plotted by changing the parameter $\alpha$ in the range of 4 to 20 with random initial conditions in the range of $[-1, 1]$ for the variable ${x}_{n}$ in part (a) and for the variable ${y}_{n}$ in part (b) of Figure 2. Also, the MLE of the model for parameter $\alpha$ is shown in part (c) of Figure 2 ($\sigma = -0.1 , \mu = 0.001$). It can be seen that with the change of $\alpha$ in the interval $[4, 4.725]$, the neuron's firing is period one (MLE is zero) and then enters chaos in $\alpha = 4.725$ until $\alpha = 8.576$. At $\alpha = 8.576$, the firing changes to periodic and again by period doubling enters another chaotic region in $\alpha = 14.688$, and finally exits from chaos in $\alpha = 19.280$. There is also a periodic window in the interval $[16.213, 16.327]$. The bifurcation diagram for variables ${x}_{n}$ and ${y}_{n}$ with respect to variation of parameter $\sigma$ in the range of $[-2.5, 0]$ are shown in parts (d) and (e) of Figure 2 ($\alpha = 4.1$, and $\mu = 0.001$). The corresponding MLE is drawn in part (f) of Figure 2. By changing σ, the dynamic behavior is periodic for the interval $[-2.5, -1.632]$, $[-1.632, -0.145]$, and then switches to chaotic. The bifurcation diagrams by changing parameter $\mu$ in the interval $[0, 1]$ for $\alpha = 4.1$, and $\sigma = -0.1$ are shown in parts (g) and (h) of Figure 2 for variables ${x}_{n}$ and ${y}_{n}$. By changing the μ parameter, the chaotic behavior is seen in the intervals $[0, 0.087]$ and $[0.377, 0.805]$.

3.   Fractional-order discrete Rulkov neuron map

In this section, Rulkov neuron map is reformulated as a fractional-order discrete system using the Caputo left difference operator. Discrete fraction order maps have recently been recognized as a new powerful tool in modeling and describing the dynamics of complex nonlinear discrete systems ^[56,57]. Different numerical methods have been developed to solve discrete fraction order maps ^[58,59]. In this paper, we examine the following Rulkov neuron map with discrete fraction order:

where ${\Delta }_{a}^{q}$ is a function of the Caputo delta difference on the left, ${N}_{a} = \left[a, a+1, a+2, \dots \right]$ represents the time scale isolated ($a\in R$ and constant). To obtain the fractional order form of the Rulkov neuron map, we first define the necessary definitions from the discrete fractional calculus and generalize the results to the Rulkov neuron map.

Definition 1. When $u:{N}_{a}\to R$, $q > 0$ and $q\notin N$, the sum of the fraction of order $q$ is defined as follows:

where $a$ is the starting point, and ${t}^{\left(q\right)}$ is the falling function, which is defined as follows:

where $\Gamma \left(t\right)$ is Gamma function, defined as follows ^[60]:

Definition 2. For $q > 0$, $q\in {N}_{a}$ and $u\left(t\right)$ defined on ${N}_{a}$, the difference in the Caputo delta fraction can be expressed as follows ^[58]:

where $m = \left[q\right]+1$ and $m = \left[q\right]+1$.

Theorem 1. For the delta fraction difference equation

the equivalent discrete integral equation can be obtained as follows:

where the initial iteration ${u}_{0}\left(t\right)$ is expressed as follows:

Detailed descriptions of the above results are given in reference ^[61].

Rulkov's neural map of discrete fraction order (Eq. 3.1, and Eq. 3.2) can be written as the following classical map:

According to Theorem 1, we can write the discrete fractional for Rulkov's model in the following form:

In conclusion, the numerical formula of the $n$-dependent discrete Rulkov fractional neuronal model is defined as follows:

The difference between the maps of the integers of Eq. 2.1, and Eq. 2.2 and the fractional order of Eq. 3.14, and Eq. 3.15 is a discrete core function and $x\left(n\right)$ and $y\left(n\right)$ refer to the past information $x\left(0\right), \dots, x\left(n-1\right)$ and $y\left(0\right), \dots, y\left(n-1\right)$. As a result, the current state depends on all past forms, which express the memory effects. We define the proposed model in Eq. 3.14, and Eq. 3.15 as the fractional-order of discrete Rulkov neuron map (FORNM). Also, in Eq. 3.14, and Eq. 3.15, the following equation is used:

By changing the parameters $\alpha$, $\sigma$, and $\mu$ of the FORNM system, different behaviors such as silence ($\alpha = 3.15$, $\sigma = -2$, $\mu = 0.2$, $q = 0.001$), bursts of spikes ($\alpha = 3.5$, $\sigma = -2$, $\mu = 0.2$, $q = 0.001$), and chaotic firing ($\alpha = 6$, $\sigma = -1$, $\mu = 0.3$, $q = 0.001$) can be seen which are shown in Figure 3. When the membrane potential reaches a threshold, the neuron fires and generates a signal that travels to other neurons, which in turn increases or decreases their potential in response to the signal. The neuron model that fires when the threshold is crossed is called a spiking neuron. Neuron spiking occurs when up-down oscillations accompany this spike train ^[62,63]. According to the stated definition, bursts of spikes occurs in part (b) of Figure 3, and this spiking behavior has been previously identical in ^[24,31]. The phase diagram for the bursts of spikes mode is shown in part (d) and for the chaotic firing in part Ⅰ.

In the following, the dynamical analysis of the FORNM is done by plotting the bifurcation diagrams versus the fractional order (q) and the system parameters at different intervals.

3.1.   Effect of the FORNM parameters

The dynamic behavior of the FORNM is firstly investigated by changing the parameters of the Rulkov neuron map. It was shown that there are chaotic, periodic, and fixed-point dynamics in the system for different values of parameters $\alpha$, $\sigma$, and $\mu$. Figure 4 shows the bifurcation diagrams wherein the left, middle, and right columns, the fractional order is set to $q = 0.001$, $q = 0.01$, and $q = 0.02$, respectively. The bifurcation parameter in parts (a-c) is $\alpha$ in the interval $[0, 12]$ for $\sigma = -1$ and $\mu = 0.3$, and in parts (d-f) is $\sigma$ in the interval $[-5, 0]$ for $\alpha = 4$ and $\mu = 0.3$, and in parts (g-i) is $\mu$ in the Interval $[0.15, 0.5]$ for $\alpha = 4$ and $\sigma = -1$. The initial conditions are selected randomly. It can be observed that the dynamics are strongly dependent on $q$. Generally, it can be concluded that by increasing $q$, the chaotic behavior is observed for higher bifurcation parameters, and the region of chaotic behavior is shrunken.

3.2.   Effect of the fractional-order

Another important factor in FORNM is the fractional order parameter, which can affect the overall dynamics. Figure 5 presents the bifurcation diagrams of the model according to the fractional order for different parameter values. The bifurcation diagram for $\alpha = 1.9$, $\sigma = -1$, $\mu = 0.3$ in part (a) of Figure 5 shows chaotic behavior for the range $[0, 0.0102]$ and periodic behavior for the range $[0.0102, 0.1]$. In part (b), in which $\alpha = 4.9$, $\sigma = -1$, $\mu = 0.3$ the chaos is seen in the range of $[0, 0.0180]$. The bifurcation diagram in terms of q for $\sigma = -2$, $\alpha = 4,$ and $\mu = 0.3$ is shown in part (c) that has chaotic behavior for the interval $[0, 0.00389]$, while for $\sigma = -0.833$ in part (d), the chaotic behavior happens in $[0, 0.01928]$. The bifurcation diagram for $\mu = 0.33$, $\alpha = 4$, and $\sigma = -1$ is shown in part (e) of Figure 5. In this case, the chaotic behavior is in the range $[0, 0.01809]$. In part (f) where $\mu = 0.45$, two chaotic regions in $[0, 0.2928]$ and $[0.1860, 0.2]$ are observed. As shown in Figure 5, increasing the parameters' values increases the chaotic region.

4.   Fixed points and stability analysis of FORNM

In this section, we calculate the equilibrium points of the proposed FORNM of Eq. 3.14, and Eq. 3.15 and evaluate them in terms of stability. The equilibrium points are calculated by:

i.e.

Putting two equations together:

i.e.

i.e.

According to the method of solving the equation of degree 3 ^[64], by solving Eq. 4.7, we obtain $a = -\sigma$, $b = 1$, $c = -\sigma +\frac{\alpha }{\mu }$, so $p = 1-\frac{{\sigma }^{2}}{3}$, and $q = \frac{\alpha }{\mu }-\frac{2\sigma }{3}-\frac{2{\sigma }^{3}}{27}$. The eigenvalues of Eq. 4.7 are equal to $\Delta = \frac{{\left(\frac{2{\sigma }^{3}}{27}+\frac{2\sigma }{3}-\frac{\alpha }{\mu }\right)}^{2}}{4}-\frac{{\left(\frac{{\sigma }^{2}}{3}-1\right)}^{3}}{27}$, so the equilibrium points of the model in terms of $\Delta$ states are equal to:

● If $\Delta > 0$,

i.e.

● If $\Delta = 0$,

i.e.

● If $\Delta < 0$,

i.e.

Thus, the equilibrium points of the system are calculated, and we use the following theorem to examine the stability conditions of the equilibrium points of the fractional order system.

Theorem 2. The zero equilibrium point of the system:

is asymptotically stable when the eigenvalues ${\lambda }_{i}, i = 1, 2, \dots, n$ of the matrix $M$ have the following two conditions ^[65]:

Therefore, we calculate the Jacobian matrix of Eq. 3.1, and Eq. 3.2, which is as follows:

The eigenvalues of the Jacobin matrix are calculated for equilibrium points at different parameters, and the conditions of Theorem 2 are obtained numerically for the eigenvalues. The parameters that satisfy the Theorem 2 requirements make Eq. 3.1, and Eq. 3.2 always stable.

4.1.   Stability regions of the FORNM

The stability regions of the equilibrium points of the FORNM model have been investigated by varying the parameters in the model. The fractional order in the FORNM model is examined in the first step. Stability regions of the FORNM model are shown in Figure 6 in the three-dimensional space of $\alpha -\sigma -\mu$ parameters and the two-dimensional space of parameters $\alpha -\sigma$, $\alpha -\mu$, and $\sigma -\mu$. The fractional order in parts (a-c) are $q = 0.001$, $q = 0.01$, $q = 0.025$ respectively. It is shown that the stable regions decrease by increasing the value of fractional order.

In the second step, the effects of $q-\sigma -\mu$ parameters with $\alpha = 15$ constant in part (a), the effects of $\alpha -q-\mu$ parameters with $\sigma = -5$ constant in part (b), and the effects of $\alpha -\sigma -q$ parameters with constant $\mu = 0.3$ in part (c) on the stability of the FORNM model are investigated and shown in Figure 7. In part (a) of Figure 7, the stability of the system is lost as the value of $\sigma$ increases toward zero, and the stability of the system is lost as $\mu$ increases toward one. In part (b) of Figure 7, the stability of the system increases with the increase of $\alpha$ value, and the stability of the system decreases with the increase of $\mu$. In part (c) of Figure 7, by increasing the value of $\alpha$ and $\sigma$, stability is created in the system.

5.   Synchronization of FORNM

Another important aspect of the dynamical analysis is the synchronization in a chaotic state. Two coupled fractional-order Rulkov neuron maps can be formulated as follows:

where subscript 1 corresponds to neuron one and subscript two corresponds to neuron 2, ${d}_{x}$ is the coupling strength of variable $x\left(n\right)$, and ${d}_{y}$ is the coupling strength of variable $y\left(n\right)$. For both neuron models, the initial conditions are randomly chosen between -1 and 1. To determine the synchronization level of two neurons, the average error is calculated in the form of the following equation:

where $N$ represents the number of time series data samples. We investigate the effects of fractional order changes on the synchronization of neurons. Figure 8 shows the synchronization error by changing the coupling strengths ${d}_{x}$ and ${d}_{y}$ in different fractional orders. It has been proven that two integer-order Rulkov maps cannot reach synchronization when connected by electrical synapses ^[66]. Here, it can be seen that the synchronization cannot occur in the fractional-order neurons too. It should be noted that the systems become unstable for coupling coefficients higher than 1. Unsynchronized neuron maps based on previous papers show that increased activity in certain areas is accompanied by decreased activity in other areas ^[67], and this pattern exists during rest and sleep ^[68].

6.   Conclusion

This paper proposed the discrete fractional order Rulkov neuron map in Caputo's concept. The dynamic analysis of the fractional model was performed, and its synchronization was evaluated. Adding discrete fraction calculus to the Rulkov neuron map can consider the effect of memory on neural model dynamics. Memory effects refer to the fact that system states are determined in fractional order by all previous states. At first, the Rulkov neuron map was examined in terms of phase plane, bifurcation diagram, and Lyapunov exponent. It was shown that the biological behaviors of the Rulkov neuron map, such as silence, bursting, and chaotic firing, also exist in the discrete fractional order Rulkov neuron map. Hence, the Rulkov neuron map of discrete fractional order shows the same behaviors as the Rulkov map of integer order. However, changing the fractional order can result in bifurcation. In the next step, the stability regions of the system from a theoretical and numerical point of view were examined and it was shown that increasing the order of the fraction decreases the stability in the FORNM system. In the last stage, the synchronization behavior of two FORNM systems was investigated. The results showed that the fractional-order systems do not achieve synchronization with an increase in the coupling coefficient, which is similar to the integer-order systems. Table 1 shows the reason for the superiority of the result of Rulkov fractional neuron map to the previous works.

For future work, discrete fractional order derivatives can be implemented on DNA-based genetic systems to investigate hereditary factors, so that human evolution can be described from these systems. For the next research topic, different types of discrete fractional order derivatives can be evaluated on the Rulkov neuron map. Due to the fact that the Rulkov discrete neuron map has no synchronization power, different controllers can be used to check synchronization.

Acknowledgments

This work is partially funded by Centre for Nonlinear Systems, Chennai Institute of Technology, India, vide funding number CIT/CNS/2023/RP/005.

Conflict of interest

The authors declare that they have no conflict of interest.