A fractional-order discrete memristor neuron model: Nodal and network dynamics

1.   Introduction

Fractional order systems can be found in a variety of disciplines, including physics, electronics, biology, and engineering ^[1,2,3,4,5]. As compared to integer order systems, such non-integer-order systems are ubiquitous and used to describe a wide range of real-time dynamics. Surprisingly, non-integer-order systems impersonate the dynamics of the real-time systems more precisely than integer order systems. Due to this fact, a wide range of studies has been dedicated to understanding the dynamics of non-integer-order systems in recent years. The researchers discovered that non-integer order models are the best approach for describing processes with memory and genetic characteristics in neural networks, and that they also facilitate information processing and stimulus anticipation ^[6,7]. Furthermore, fractional order calculus with neural networks improves computation accuracy and flexibility and it takes part in applications like computation optimization and control performance ^[8]. As a result, the fractional order with neural networks has been taken into consideration systematically.

As a consequence, various studies have been performed using fractional order neural networks, revealing a variety of intriguing phenomena such as synchronization, chimera, and clustering, among others ^[9,10,11]. For example, synchronization behavior has been discovered in a non-integer order delayed neural network ^[12]. Adaptive synchronization emergence in a memristor neural network with different time delays has been demonstrated in ^[13]. Recently, synchronization and its transitions were reported in a fractional Hindmarsh-Rose (HR) neuron model ^[14]. Furthermore, the emergence of chimera states has also been observed in fractional order systems, such as fractional HR neurons and the magnetic HR model with a fractional derivative ^[15,16,17]. It was also discovered that strong coupling with decreasing derivative order reduces the level of synchronization with the emergence of chimera ^[18]. Recently studies have been showing the quasi-synchronization multi-layer fractional order networks with mismatched system parameters, as well as delay-dependent impulsive feedback ^[19] and fractional neural networks with leakage and discrete delays ^[20]. Controlling the dynamics of fractional order systems has been of interest in recent years, for example, the exponential stability problem for fractional order networks with intermittent control ^[21,22].

Among the various spatiotemporal patterns^[23], spiral waves are an intriguing phenomenon reported in many excitable chemical and biological systems ^[24,25,26]. Originally, it was reported in the Belousov-Zhabotinsky reaction ^[26]. Later, the spiral wave pattern was discovered in various systems including map lattices, neuronal systems, chaotic systems and so on. For instance, the occurrence of diverse patterns such as spiral waves, periodic, domain patterns and banded spiral waves are observed with different shapes such as polygonal shapes and domain patterns in a two-dimensional coupled map lattice of a sine circle ^[27]. Different types of wave propagation in a heterogeneous $2d$ lattice of coupled Nekorkin maps were investigated and similar phenomenon were observed in other systems such as reaction-diffusion models, FitzHugh-Nagumo models, phase oscillators and chemical oscillators ^[28]. It was also discovered in multilayer networks, where the two-layer Fitzhugh-Nagumo neural network exhibits high excitability in one layer and low excitability in the other layer because of the inter-layer coupling interaction ^[29]. In ^[30], the authors have also shown that the coupling parameters, external current, and amplitude of external stimuli are key parameters for the generation of spiral waves in the fractional order memristive synapse HR model. Very recently, the existence and the suppression of spiral waves were manifested through the $2d$ lattice of the fractional order neuron model.

With the knowledge of the above works, this work investigates the spatiotemporal pattern in fractional order flux coupled discrete neuron map (FDN) lattices. We believe that the present system is the simplest one-dimensional neuron model which shows all of the neuronal activities like spiking and bursting. The present model also has the merits of fractional order which has the memory function and flux coupling. To the best of our knowledge, no other model has shown a low-dimensional neuron map with memristive function and flux couplings. Primarily, the dynamical transitions are examined through a bifurcation analysis. Further, the observed bifurcation transitions are validated through the Lyapunov spectrum and entropy analysis. Followed by single system analysis, the spatiotemporal dynamics are investigated for the network of map lattices. In particular, we investigate the spatiotemporal behavior by fixing the different values of flux coupling $(k)$, fractional order $(q)$ and coupling strength $(D)$. We discover that increasing the flux coupling strength $k$ leads to multiple spiral wave zones as well as the spiral wave becomes highly turbulent at higher flux coupling values. When the fractional order $q$ is increased, the entire network loses its excitability and becomes turbulent, which results in distributed spiral waves. Similarly, spiral waves appear in the entire network only when the coupling strength $D$ increases to a sufficient range, and after that it becomes unstable.

The structure of the article is as follows: In Section 2, we introduce our considered model of FDN. The corresponding dynamical transitions are investigated through bifurcation analysis, Lyapunov spectrum analysis, and entropy analysis in Section 3. The existence of spiral waves is demonstrated for the $N\times N$ lattice network via spatiotemporal analysis in Section 4. Finally, we present our summary of the work and conclusions in Section 5.

2.   The Model

In ^[31], a simple modified Nagumo-Sato neuron model ^[31,32] was discussed the mathematical model is defined by a one-dimensional discrete map

The parameters $a$, $b$ and $c$ are the dissipation, weight and bias values, respectively. The value of the parameter $b$ is calculated by combining the initial bias current and the external stimuli current. The model in (2.1) exhibits rich dynamics as discussed in ^[31] and it is designed such that the discontinuity caused by the Heaviside function in the original model is now eliminated. Recently in ^[33] the authors have discussed the fractional order model of (2.1) and the modified mathematical model is defined as

where the function $M(\phi_n)$ is the memductance of the discrete memristor ^[34,35] and it represents the internal electric charge of the discrete memristor, which can be expressed as

Further expanding the above equation using the backward difference method, we end up with the following equations:

The memductance function can be defined as $M(\phi_n) = \alpha+\beta \phi_n^2$ where $\alpha$ and $\beta$ are the memristor internal parameters. The fractional order model of (2.2) can be mathematically expressed as

Here, the $\Delta$ represents the fractional order operator, $q$ refers to the order of the fractional derivative and $c$ is for the Caputo method. The initial value problem of the model (2.5) is solved using the Caputo delta difference method which can be expressed as

where $n \in N_{1-q}$. The following equation describes the discrete form of (2.6).

where $\Gamma$ is the Euler-gamma function. After several steps, we may reduce (2.7) to

The fractional order discrete neuron model with flux coupling (2.8) is henceforth called as the FDN. We use the method used in (2.3) and (2.4) to solve the discontinuity issue by replacing

3.   Dynamical analysis of fractional order flux coupled discrete neuron map

3.1.   Bifurcation and frequency spectrum analysis as a function of the order $q$ and parameter $k$

The dynamical behavior of the FDN can be investigated by using bifurcation plots using the local maxima and frequency, and via two-parameter bifurcations using Lyapunov exponents. We also use sample entropy and permutation entropy to analyze the complexity of the FDN model. We fixed the system parameters to $a = 0.8, b = 4, c = -16, \alpha = 0.8$ and $\beta = 0.01$. First, we investigate the dynamics of the FDN for various $q$ values by fixing $k = 0.1$. The bifurcation plot for $q$ is shown in Figure 1 and we can see two larger distinct regions showing chaos. The chaotic dynamics in the system disappear through the period by using the halving bifurcation (also known as inverse period route) route. We have also used the frequency bifurcation plot (Figure 1 (b)) to show the existence of chaotic and periodic regimes. The dominant frequencies (which are the Dirac-delta peaks or discrete peaks present in the power spectrum) in the periodic domains are identifiable compared to the periodic oscillations in their respective time-domain bifurcation plot (Figure 1(a)). The time series bifurcation is drawn by collecting the local peaks of the time series at every parameter on the x axis. While the chaotic dynamics of the system can be inferred from the continuous spectra (broad band spectra) from the fast Fourier transform (FFT) bifurcation diagram. In Figure 2, we have provided the bifurcation plot for the FDN with the control parameter of the flux coupling constant $k$ with the fixed value of the order of the fractional derivative $q = 0.01$. The dynamics show two regions of chaos and take an inverse period doubling to route from a chaotic state to periodic dynamics. The FFT bifurcation shown in Figure 2 (b) agrees with the corresponding time series bifurcation plot. Lyapunov exponents have been used to investigate the chaos diagram of FDN and which is shown in Figure 3. We can see the existence of positive Lyapunov exponents in two regions say the first region for $0.0001\leq q\leq 0.3, 0\leq k\leq 0.4$ and the second region for $0.4\leq q\leq 1, 0.2\leq k\leq 1$. The exact combination of $q$ and $k$ for chaotic regions can be seen in Figure 3. The color bar in the plot shows the value of the maximal Lyapunov exponent. The positive value in the figure corresponds to the chaotic dynamics of the FDN system. The stable oscillation of the system can be identified from the Figures 1 and 2, where we have discrete periodic oscillations. Similarly, from Figure 3, we can see that the value of the Lyapunov exponent to distinguish the periodic and chaotic dynamics, where the dark blue color shows the regular behavior.

3.2.   Entropy analysis

3.2.1.   Permutation entropy

Many works of literature ^[36,37,38] have shown that entropy analysis of time series is one of the best tools for investigating nonlinear system properties. Hence, we use permutation entropy ^[24,25] and sample entropy (SampEn) ^[36] as tools to discuss the complexity measure of the FDN model. We have subdivided the ways to calculate the permutation entropy from a time series data into several steps as below.

$\ast$ Step 1: Phase space reconstruction for a given time series ${w_i, \; i = 1, 2, \cdots, N}$,

where time delay and embedded dimensions are given by $\tau$ and $m$ respectively.

$\ast$ Step 2: Arrange the real values in each $W(i)$ in ascending order is

$\ast$ Step 3: Group the values in $W(i)$ as the symbols $S(l) = (k_1, k_2, …, k_m)$ and for example if $w(i+(k_1-1)\tau) = w(i+(k_2-1)\tau)$, then the sorting of these values can be in the order of $k_1\leq k_2, w(i+(k_1-1)\tau)\leq w(i+(k_2-1)\tau)$. The term $l = 1, 2, …, k$ where $k\leq m!$ with $m!$ being the substantial number of distinct symbols. $S(l)$ is one of the $ml$ symbol permutations and it is mapped onto the $m$ number symbols $(k_1, k_2, …, k_m)$ in $m$-dimensional embedding space.

$\ast$ Step 4: Denote the probability distribution function (PDF) of each symbol sequence as $P_1, P_2, \cdots, P_k$ such that $\sum_{(l = 1)}^kP_l = 1$, and the Shannon entropy for $k$ symbol sequences as

$\ast$ Step 5: Using the maximum value of $H_p(m)$ by calculating $ln(ml)$ such that all those symbol sequences have the same PDF as $P_l = 1⁄m!$; the permutation entropy can be defined as

3.2.2.   Sample entropy

Sample entropy has been used to categorize dynamics from a signal of a nonlinear time series. It is a modification of the version of approximate entropy. It has been used because of its advantages over the approximate entropy, which are (ⅰ) the data length independence compared with other methods and (ⅱ) it is a relatively trouble-free method of implementation. The interpretation of the outcome of the method is very straightforward. Larger values denote higher complexity while smaller values characterize more self-similar and regular signals. We used the time series analysis of the membrane potential to calculate the SampEn. We propose the following steps for the calculation of SampEn.

$\ast$ Step 1: Construction of phase space values from the time series of membrane potential

where $i = 1, 2, \cdots, N-m+1$. The relation used to calculate the Euclidean distance is given by

where $i, j = 1, 2, \cdots, N-m+1$ and $i\neq j$. The complexity measure can be calculated from a given threshold value $T$:

where $\psi$ denotes integer number that is $i = 1, 2, …, N-m+1$ given that $i\neq j$. The intermediate coefficient is calculated as

Considering $m = 1$, we redo the above procedure presented as (2.8)–(3.2); the SampEn is given by

The parameter $m$ represents the dimension and $T$ represents its tolerance, respectively. The tolerance is calculated using T$_t$ = T$\times$SD where SD being the standard deviation of the given time series.

In Figure 4, we have presented the entropies of the FDN for different values of $q$ and $k$. To calculate the SampEn, we used the time series of the same variable and we used the same set of initial conditions. Figure 4(a) is the plot of SampEn, of the time series $w$ as a function of order of the fractional derivative $q$ and Figure 4(b) shows the SampEn, as a function of the parameter $k$. We can compare Figure 4 with the plots shows in Figures 1 and 2. We could note that the regions showing chaos have larger entropies while the periodic regions have lesser entropies.

Next in Figure 5(a) we have calculated the permutation entropy of the FDN model and the larger entropy values confirm the existence of chaotic regions. In Figure 5(b) we have shown the entropy calculated using the SampEn and both Figure 5(b) and Figure 5(b) agree with each other; SampEn gives a clear view of the dynamical regions of the system. Comparing Figure 5 and Figure 3, we could observe that the Lyapunov exponent chaos diagram cannot identify certain regions of complexity shown in blue (Figure 5(a)) and cyan (Figure 5(b)). These regions may be the existence of strange nonchaotic attractors which have complexity more than periodic attractors and less than chaotic attractors and whose Lyapunov exponent is negative.

4.   Spatio temporal dynamics of the FDN network

In this section, we present the spatio-temporal dynamics of the proposed neuron map with electromagnetic flux coupling. We consider the same parameters given in Section 3.1 for simulation. The mathematical model of $N\times N$ lattice network is given by

In the $N\times N$ network of neurons, we would like to investigate the spatiotemporal dynamics with electromagnetic flux coupling. A $N\times N$ lattice array of the proposed neuron map is constant as in (4.1). A periodic stimulus is applied to the network from the center of the network by fixing $\theta_1 = \theta_2 = 50$.

In the first case, the control parameter electromagnetic flux coupling coefficient $k$ is used to capture the spatiotemporal dynamics. The observed dynamics are plotted in Figure 6. When $k = 0$, the network behaves like a one-dimensional neuron map and displays small regions of spiral waves which can be seen in Figure 6(a). When $k = 0.01$ the spatiotemporal dynamics of the network display much stronger spiral waves as shown in Figure 6(b). Increasing the flux coupling to $0.05$ creates many local pools of highly excitable nodes and thus creates multiple spiral wave regions in the network. As we keep increasing the flux coupling, the network becomes highly turbulent as seen in the case of $k = 0.1$.

In the second case, the control parameter is the order of the fractional derivative $q$ and the captured structure of the network is shown in Figure 7. For $q = 0.001$, the lattice array shows a spiral wave in the network with a symmetric pair formed around the corner of the lattice. The excitability of the network is heterogeneous which can be seen by the formation of spiral waves with different colors confirming different levels of excitability. Increasing the fractional order value from $0.001$ to $0.01$, the spiral waves increase in count whereas the radii of the spiral waves decrease. For $q > 0.01$, the entire network loses its excitability and exhibits turbulent and distributed spiral waves confirming that the nodes are having very close excitability voltages. This is due to the nodes extending their periodic oscillation and having unsynchronized behavior in the network. This can be due to the emergence of chimeras in the network with few nodes in synchronized states and this is highly excitable while the remaining in the unsynchronized state have low excitability.

Finally, we vary the coupling strength $D$ and the observed snapshots of the lattice array are given in Figure 8. The nodes are so weakly coupled for the coupling strength $D < 0.5$, making the entire media dissipate the applied stimuli because of the different excitabilities of the nodes. Hence, we made an effort to enhance the stimuli strength to explore whether the wave can travel to the boundaries but it was in vain-as shown in Figure 8(a) and (b). The wave irrespective of the stimuli amplitude cannot travel to the boundaries of the lattice network. Rising the coupling to $D = 0.5$, we could note the formation of multiple smaller spiral arms (Figure 8(c)) which consolidate large rotating spirals when increasing the coupling to $D = 0.6$ (Figure 8(c)). An increase in coupling strength further makes the network unstable.

5.   Conclusions

A fractional order discrete neuron map with electromagnetic flux coupling has been proposed. A fractional order discrete memristor was modeled to mimic the flux effects on the neuron model. Both the time domain and frequency domain bifurcation analysis were considered to study the dynamical properties of the proposed FDN. When considering the fractional order and the flux coupling coefficients as the control parameters, the FDN takes an inverse period-doubling exit from chaos. To show the different dynamical regions of the FDN in the $k$-$q$ plane, we have used the Lyapunov spectrum method. As an alternative, we have used the permutation and the sample entropy methods to derive the different dynamical regions in the same $k$-$q$ plane, and interestingly, we could identify several regions of complex behavior where the Lyapunov exponents are negative. Such regions confirm the existence of strange nonchaotic attractors which will be investigated in detail in our next work on the same FDN. To understand the wave propagation in the FDN network, we constructed a $N \times N$ lattice of FDN nodes and applied a periodic external stimulus to the center of the network.

The network of the system shows a rich variety of dynamics. In order to study the network, we considered the fractional order $q$, flux coupling coefficient $k$ and the network coupling strength $D$ as the control parameters. When considering the flux coupling coefficient, we could show that for lower values of $k$ the network supports spiral waves as the nodes are in high excitability. Increasing the $k$ values makes the network turbulent and for $k = 0.1$, the network becomes unstable. For the fractional order, the network supports spiral waves for $q < 0.1$ as because the nodes exhibit high excitability and heterogeneity in these regions. Interestingly, neither a strong coupling $(D > 0.7)$ nor a weak coupling $(D < 0.3)$ could trigger the formation of spiral waves in the network. This is one of the notable observations in the present study. For the coupling $0.5\leq D \leq 0.6$ the nodes are in high excitability confirming the formation of spiral waves. The present work gives a simple fractional order neuron model which exhibits almost all the dynamical behaviors observed in real neurons. The rich dynamics observed by the network of neurons show the excitability of the network with different parameter choices. The proposed system and results can be used to understand various aspects of network dynamics like chimera states, amplitude death, and short-time synchronizations.

Acknowledgments

The researchers would like to thank the Deanship of Scientific Research at Qassim University for funding the publication of this project.

Conflict of interest

The authors declare that there is no conflict of interest.

Ethics statement

This research did not involve human or animal subjects.

Data accessibility

The data that supports the findings of this study are available from the corresponding author upon reasonable request.

Authors contributions statement

All authors have contributed equally to the research and the writing up of the paper.