Firing Activities in Fractional-Order Hindmarsh-Rose Neuron with Multistable Memristor as Autapse

: Compared with integer order neurons, fractional-order neuron model can more accurately describe the firing behavior of biological neurons. Considering the fact that memristors have the characteristics similar to biological synapses, a fractional-order multistable memristor is firstly proposed in this study. It is verified that the fractional-order memristor has multiple local active regions and multiple stable hysteresis loops, and the influence of fractional order on its nonvolatility is also revealed. Then by considering the fractional-order memristor as an autapse of HR neuron model, a fractional-order memristive neuron model is developed. The effects of the initial value, external excitation current, coupling strength and fractional order on the firing behavior are discussed by time series, phase diagrams, Lyapunov exponents and inter spike interval (ISI) bifurcation diagrams. Three coexisting firing patterns, including irregulate A-periodic bursting, A-periodic bursting and chaotic bursting, dependent on the memristor initial values are observed. It is also revealed that t he fractional order can not only induce the transition of firing patterns, but also change the firing frequency of the neuron. Finally, a neuron circuit with variable fractional order is designed to verify the numerical simulations.


Introduction
As an important basic unit of the nervous system, biological neuron has rich and complex firing activities. In order to study the firing patterns and operating mechanisms of biological [21]. The dynamics and chaotic synchronization of two fractional coupled HR neurons with time delay under electromagnetic radiation were explored by Meng et al [22]. Their research results show that the fractional order has an important influence on the firing patterns. Similarly, the effects of fractional order on synchronization of neuronal networks under different external conditions were reported in [23][24][25].
Local activity is the origin of complexity [26]. In 2014, Chua proposed the first locally active memristor [27]. Compared with ideal memristors, local active memristors have more complex nonlinear dynamic characteristics [28][29][30]. The multi-stability characteristic of local active memristor [31] make it a good means to study the conversion and coexistence of various firing patterns in neuron models. Thus, local active memristor based neurons and neuronal networks have been extensively studied in recent years. For example, Li et al. constructed a bistable neural network based on local active memristor, and revealed the switching and bifurcation behaviors of firing patterns in the neural network [32]. Lin et al. constructed an HR neural network based on a tristable locally active memristor, and found that the proposed neural network produces three multi-stable phenomena with four different firing patterns simultaneously [33]. A model of HR neuron with four-stable locally active memristor as selfsynapses was proposed by Li et al. and they found that four stable firing patterns in the proposed neuron can be switched freely by selecting appropriate initial conditions [34]. With further development of memristor research, locally active memristor has been extended to fractionalorder form, and the related research has been carried out. Xie et al. proposed a fractional-order memristor with infinite locally active interval and coupled it to Chen's chaotic system. They found that the system presents different states under different fractional orders [35]. However, so far, there is no relevant research on fractional-order neuron model with locally active memristor.
Inspired by the above studies, a novel fractional multistable locally active memristor is proposed in this paper, and a new neuron model is constructed by coupling it with fractionalorder HR neuron. The main contributions of this work are summarized as follows. Firstly, a novel fractional-order locally active memristor is developed. Its nonvolatile memory and locally active characteristics, initial value-dependent coexisting pinched hysteresis loops and fractional order-dependent electrical behavior are analyzed. Secondly, firing activities of the fractional-order memristive neuron model under different external conditions are studied.
Finally, physical implementation circuit of the fractional-order neurons is proposed and the circuit simulation verification is performed.
The layout of this paper is as follows. In Sec. 2, a new fractional tristable locally active memristor is introduced and its related characteristics are studied in detail. In Sec. 3, a 4D fractional-order neuron model, by considering the fractional-order multistable locally active as an autapse of 3D-HR neuron, is constructed, and the stability of the neuron model is analyzed.
In Sec. 4, the influences of coupling strength, external excitation current and fractional order on the fractional-order neuron model are discussed. In Sec. 5, the circuit simulation of the proposed fractional-order neuron model is carried out, and the conclusion is given in the last section.

Model and characteristics of fractional-order multistable memristor
It is well known that there are three definitions for fractional calculus, i.e. Grünwald-Letnikov (GL) definition, Riemann-Liouville (RL) definition and Caputo definition. Among them, Caputo definition is very suitable for real world applications due to its initials with the same form as those of its integer order counterparts. In this paper, the Caputo definition, thus, is used to characterize the neuron model and the Caputo definition of f (t) is described as: where, 0 < α <1 is the fractional order, n is the first integer which is not less than α and Г() is Euler gamma function.

Fractional-order multistable memristor model
In this paper a new fractional-order multistable local active memristor is proposed which is defined by the state-dependent Ohm's law where i and v represent the input current and input voltage, φ the is internal state variable and G(φ) = nφ is the defined memductance function. The corresponding state equation can be described as where  is the memristive parameter. It is worth noting that the fractional-order state equation is constructed with sign functions which are related to the internal state φ and the terminal voltage v. Compared with the existing locally active memristors whose state equations are implemented with polynomial functions [28,[30][31], the proposed fractional-order locally active memristor is more convenient for circuit implementation with low cost. In the following, we will study the basic characteristics of the memristor from the aspects of pinched hysteresis loop, nonvolatile memory and locally active characteristics.

Pinched hysteresis loop
Pinched hysteresis loop is one of the important characteristics of memristor. When the parameter β = 0.2 and the order α = 0.9, a sinusoidal voltage v = Asin(2πFt) with amplitude A = 2 is applied to the memristor, the pinched hysteresis loops with different excited frequencies are shown in Fig.1(a). It can be seen from Fig.1(a) that when the input frequency F increases In addition, when A = 2 and F = 0.314 remain unchanged, three stable pinched hysteresis loops can be obtained under different initial states φ0 = 4, 0, and 4, as shown in Fig.1(b). In order to study the influence of fractional order α on the pinched hysteresis loop, keep A = 2, F = 3.14, and memristor initial value φ0= 0, the fractional order α is selected as 0.98, 0.9, 0.8, 0.7 and 0.5, respectively, the hysteresis loops of the memristor are shown in Fig.1(c). It is observed that the smaller the fractional order is, the larger the area of the hysteresis loop is, and thus the better memory of the memristor is [36].

Nonvolatile memory and locally active characteristics
Non-volatility is also a key characteristic of memristor, which means memristor can keep its latest memductance value when the power is switched off. Usually, the POP (Power-Off Plot) is used to judge whether an electronic device is nonvolatile [37]. Let the input voltage v = 0, and the state equation (3) becomes: When order α = 0.9, the POP of the fractional memristor is shown in Fig.2(a). It can be seen are unstable equilibrium points. For arbitrary initial value (0), the memristor will eventually tend to one of the three stable equilibrium points, Namely Thus, the proposed memristor is nonvolatile and it acts as a discrete memory device, which can encode three memory states.
DC V-I curve is an important tool to judge whether the memristor has local activity. When the DC V-I curve of the memristor has a negative slope, it is locally active [38]. Let Dt α φ = 0, Eqs. (2) and (3) can be written as (sgn( 2) sgn( 2) 0.5 ) where V, I denote DC voltage and current, respectively, and Φ is a variable equilibrium state which satisfies = ( = ) = 0. By virtue of Eqs. (6) and (7), the DC V-I curve of the memristor is obtained and shown in Fig.2 (b). Observation from Fig. 2(b) that the slopes of the red, green and blue curves are negative, so the fractional memristor is locally active and suitable to realize electronic synapses.

Fractional-order memristive HR model
According to the Caputo definition, the fractional-order HR neuron model can be described as the following differential equations  (9) where x represents the membrane potential, y is the recovery or spiking variable, z is the adaptation variable, I is an external excitation current, α is the system order and a, b, c, d, r, s are system parameters. In this study, we consider the proposed locally-active memristor governed by (2) and (3) as an autapse of the fractional-order HR neuron (9), then a new 4D neuron model is developed, which is formulated as where k is coupling strength, I represent the external excitation current. In this paper, the parameters are determined as a = 1, b = 3, c = 1, d = 5, r = 0.0021, s = 4, and β = 0.2. The effects of external excitation current I, fractional order α and coupling strength k on the firing dynamics of the memristive neuron model are discussed in detail in the Section 4.

Stability analysis of equilibrium points
In order to analyze the stability of the equilibrium point of the fractional-order neuron, the following theorem is introduced.
Lemma 1 [39]: For a fractional-order system with order 0< α ≤1: is fractional-order system equilibrium point, if all the eigenvalues ( = 1,2,3, … , ) of the Jacobian matrix satisfy the following condition: the equilibrium point is asymptotically stable.
Let the right side of Eq. (10) equal to 0, the equilibrium point of the system can be written as E (x0, 15x0 2 ,4(x0+1.6), φ0), where x0 and φ0 are determined by Eq. (13) It can be seen from Eq. (13) that the equilibrium point of the system is related to the external excitation current I and the coupling strength k. The Jacobian matrix of system (10) at equilibrium point E is expressed as 0 0 1 1 where The characteristic equation of system (10) at equilibrium point E is a B x ab With k = 0.1, I = 2.7, and other parameters remain unchanged, the two expressions in Eq. (13) can be graphically expressed, as displayed in Fig.3. As we can see from are not equilibrium points. The system, thus, has three equilibrium points E1, E2 and E3.
In order to discuss the influence of external excitation current I and the fractional order α on the stability of the equilibrium points, we list the eigenvalues of the three equilibrium points  According to Lemma 1 and Table 1, we can find when the real parts of the eigenvalues of the equilibrium points are all less than 0, the neuron is stable and independent of the fractional order α. When the eigenvalues of the equilibrium point have both positive real part and imaginary parts, the stability of the equilibrium point is decided by the fractional order α. When the eigenvalues of the equilibrium point of the neuron have pure positive real roots, the system is unstable and independent of the fractional order α.

Dynamic analysis
In order to study the dynamic behaviors of the neuron, nonlinear dynamical analysis tools such as time series, phase diagram, bifurcation diagram, Lyapunov exponent and ISI bifurcation diagram are used to analyze the firing activities. In the process of numerical analysis, the system parameters are fixed as a = 1, b = 3, c = 1, d = 5, r = 0.0021, s = 4, β = 0.2. It should be noted that Adomian fractional-order algorithm [40] is adopted in this paper, the step length of simulation h = 0.01, the computing software used is Matlab2019b, and the computer system is based on win10 home edition.
What needs special attention is that many studies have shown that autonomous and nonautonomous fractional-order systems have not any exact periodic solution [41][42]. However, many studies have shown that there exist asymptotically periodic oscillations in fraction-order systems, which are named "Asymptotically T-periodic oscillations" [43][44][45][46]. In this paper, we also find the existence of asymptotically periodic oscillations in the proposed fractional-order neuronal model. To facilitate recording, the term "Asymptotically periodic" is written as "Aperiod" in the following.

Coexisting firing behavior
In order to explore the coexisting firing behavior of the fractional-order neuron model In order to verify the above three different firing patterns, the initial value of the memristor, φ0, is considered as the controlled parameter and the other parameters remain unchanged, the bifurcation diagram and Lyapunov exponents with respect to the initial value φ0 are obtained in Fig.5 (a) and (b), respectively. For clarity, the other two Lyapunov exponents LE3 and LE4 in Fig.5 (b) are omitted since their values are much less than 0. When the initial value φ0 locates at the region of (6, -2) and (2, 2), the system is A-periodic, while when the initial value φ0 is in the region of (2, 6), the system presents a chaotic state. From Fig.5 (b), one can find that the maximum Lyapunov exponent is approximately equal to 0 when the initial value φ0 locates at (6, 2) and (2, 2), meaning the neuron is in a periodic firing state; while the maximum Lyapunov exponent is greater than 0 when the initial value range φ0 is in (2,6), showing that the neuron produces chaotic firing behavior. So, the bifurcation diagram and the Lyapunov exponents are well consistent.  (1, 0, 0, -4), (1, 0, 0, 0) and (1, 0, 0, 4)).

Influence of excitation current
In the real world, neurons respond differently to different external stimuli. In this part, we In a fractional-order system, if the trajectory of the system is A-periodic oscillation in some intervals, while the trajectory of other adjacent regions is chaotic, this phenomenon is called "local chaos" [35]. From Fig.6(a), as the external excitation current I increases from 1.5 to 2.6, the fractional-order neuron gradually increases from A-periodic-1 to A-periodic-4. When I increases from 2.6 to 3.1, the neuron appears an A-periodic-chaos-A-periodic switching process.
Namely, the neuron presents "local chaos". The Lyapunov exponents of the system are shown  Fig.6(b), which are well consistent with the ISI bifurcation diagram.

Influence of coupling strength
Set the fractional order α = 0.9,  = 0.2, I = 2.7, and the initial value of the system is (1, 0,0,0), the coupling strength k is considered as the controlled parameter to explore its influence on the firing dynamics of the neuron model. When the coupling strength k increases from 0.1 to 1.4, the ISI bifurcation diagram and the corresponding Lyapunov exponents are illustrated in Figs.7 (a) and (b), respectively. From Fig.7(a), when coupling strength k (0,0.5), the system is in an A-periodic state, when coupling strength k (0.5,0.95), the system alternates between A-periodic and chaotic states, and as the coupling strength k increases from 0.95 to 1.4, the neuron undergoes a series of period doubling bifurcations and finally shows A-period 1 oscillations. Fig.7(b) illustrates the corresponding Lyapunov exponents.

Influence of fractional order
The firing frequency of neuron has great influence on its performance. The more spikes in per unit of time, the faster the action per unit time, and the faster the discharge frequency of neuron. In order to study the effect of fractional order α on the firing frequency of the fractionalorder neuron model, set k = 0.1,  = 0.2, I = 2.7, and the fractional order α varies from 0.55 to 1, the ISI bifurcation diagrams of the fractional neuron with the initial value of (1,0,0,4) is show in the Fig.8.
From Fig.8, when α increases from 0.55 to 0.7, the neuron exhibits A-periodic firing and (a) the ISI of the neuron increases gradually, meaning the firing frequency decreases gradually.
Figs.9(a) and (b) illustrate the time series of the membrane potential x with the fractional order α = 0.6 and 0.68, respectively. In the duration of t = (1000, 1200), Fig.9(a) contains 6 complete A-periodic burster, while there are just 4 complete A-periodic burster in Fig.9(b). So, we can intuitively observe that the firing frequency decreases as the increase of the fractional order α.
However, from Fig.8, as the fractional order α increases from 0.7 to 1, the firing pattern of the neuron changes from A-periodic to chaotic firing. In order to verify the transition of firing pattern, the time series of the membrane potential x with α = 0.7 and 0.8 are provided in Figs. 10(a) and (b), respectively. And the corresponding phase diagrams are shown in Fig. 10(c) and (d). It is easy to find when the fractional order is 0.7, the neuron exhibits A-periodic 5 firing, while the neuron is in chaotic firing state when the fractional order is 0.8.
Comparing the dynamical Eq. (16) with the neuron model (10), when the circuit parameters and system parameters satisfy the following relations Then the circuit equations (16) are equivalent to those of (10), thus, the proposed circuit can simulate the firing dynamics of the fractional-order neuron. According to Ref. [47], the equivalent expression of α-order capacitor in complex frequency domain can be written as: The corresponding resistance and capacitance values with different fractional order α are   Table 2. In the fractional-order capacitor module, the switching elements are used to realize different fractional orders. For example, when the switch K1 is closed and the other switches is turned off, M is a fractional-order capacitor with α = 0.98. In this way, we can obtain the other three fractional-order capacitors.

Circuit simulation results
Referring to the circuit schematic diagram provided in Fig.11 Firstly, we simulate the electrical behavior of fractional-order memristor by using Multisim 14.0. Closing K2 so that the fractional order α = 0.9, applying a sinusoidal voltage with amplitude A = 2V and frequency F = 3.14 Hz to the fractional-order memristor, we obtain three stable pinched hysteresis loops with different initial values. As illustrated in Fig.12 (a), the red, green and blue hysteresis loops are related to VC1 = 4V, 0V and 4V, respectively. Thus, the circuit simulation results effectively verify the numerical simulation results shown in Fig.1(b), indicating that the proposed fractional memristor can indeed exhibit the coexisting multiple stable hysteresis loops. Let VC1 = 0V, and close K4, K3, K2 and K1 sequentially, the corresponding pinched hysteresis loops corresponding to different fractional orders are provided in Fig.11(b), which has the same results as the numerical simulation results as shown in Fig. 1(c).
Secondly, the coexisting firing patterns are also verified by Multisim circuit simulation.
close K2 so that the fractional order α = 0.9, let VC1 = 4V, 0V and 4V, respectively, the corresponding time series of the membrane potential x with different initial values are shown in Figs.13(a)-(c) and the corresponding phase portraits in z-x plane are illustrated in Fig.13(d).
From Fig.13, one can find that the neuron circuit exhibits three coexisting firing patterns, namely, irregular A-periodic bursting, A-periodic bursting and chaotic firing, which are well consistent with the numerical simulation results illustrated in Fig.4.  Finally, the firing activities of the neuron circuit dependent on different fractional orders are plotted in Fig.14. For simplicity, we only consider two cases, namely, α = 0.7 and α = 0.8.
When only the switch K4 is closed, the corresponding fractional order is α = 0.7, we find the neuron circuit exhibits A-periodic bursting with 5 spikes, the corresponding time series and phase diagram are plotted in Fig.14 (b) and (a). While K3 is on, which means α = 0.8, we can observe complex chaotic firing activities, as illustrated in Fig.14 (c) and (d). From the numerical simulation results in Fig. 10 and the circuit simulation results in Fig. 14, it is easy to find that the fractional order has an important effect on the firing activities of the neuron.

Conclusions
In this paper, a multistable fractional-order locally active memristor is proposed, and its nonvolatile and locally active characteristics are studied and analyzed. Then, a 4D fractional order HR neuron model is developed by introducing the designed memristor as an autapse of a fractional order 3D HR neuron model. The effects of initial condition, coupling strength, external excitation current and fractional-order α on the firing activities of the fractional-order neuron model are discussed in detail by numerical and circuit simulations. Some complex firing Furthermore, we find that the fractional order α can not only affect the firing patterns but also change the firing frequency of the neuron. It is expected that this study will be helpful to the theoretical analysis and understanding of fractional-order local active memristor and the firing activities of fractional-order neuron.