Coexisting Firing Patterns in an Improved Memristive Hindmarsh–Rose Neuron Model with Multi-Frequency Alternating Current Injection

With the development of memristor theory, the application of memristor in the field of the nervous system has achieved remarkable results and has bright development prospects. Flux-controlled memristor can be used to describe the magnetic induction effect of the neuron. Based on the Hindmarsh–Rose (HR) neuron model, a new HR neuron model is proposed by introducing a flux-controlled memristor and a multi-frequency excitation with high–low frequency current superimposed. Various firing patterns under single and multiple stimuli are investigated. The model can exhibit different coexisting firing patterns. In addition, when the memristor coupling strength changes, the multiple stability of the model is eliminated, which is a rare phenomenon. Moreover, an analog circuit is built to verify the numerical simulation results.


Introduction
Neurons are the basic unit of the nervous system structure, and information is transmitted in the nervous system in the form of electrical pulses.The earliest neuron model is the Hodgkin-Huxley (HH) neuron proposed in 1952 [1].Subsequently, the FitzHugh-Nagumo (FHN) neuron model [2], Morris-Lecar (ML) neuron model [3], HR neuron [4], and Hopfield neural network [5], which are derived from the HH neuron model, were gradually proposed.In recent decades, the dynamic behavior based on various neuron and neural network models has been studied.Among them, HR neurons can effectively simulate the electrical activities in the brain [6].Many biological neural electrical activity phenomena such as silence, spiking firing, bursting firing, periodic oscillation, and even simple chaos have been verified in the neural model.To study the influence of parameter variation on the firing patterns of neural nonlinear dynamical system, an improved Hindmarsh-Rose neural nonlinear dynamical system model was proposed in [7].Various firing patterns of neurons can also be observed by changing the initial conditions as discussed in [8], which indicates the memory effect of the neuronal system.The study of firing patterns and chaotic dynamics in neurons has gradually become a popular issue in the international academic community [9].
Furthermore, with the continuous development of memristor theory, the application of memristor in many fields has achieved remarkable results and has bright development prospects [10][11][12][13][14].It has developed rapidly in the fields of memristor neurons and neural network dynamics.When a neuron is stimulated by electromagnetic radiation, it will accumulate magnetic flux on the cell membrane to generate a magnetically induced current, and the memristor correlates the relationship between magnetic flux and electric charge.
Therefore, the flux-controlled memristor can be used to describe the magnetic induction effect, in order to establish a neuron or neural network model affected by electromagnetic radiation.At present, the research on using flux-controlled memristor to describe the electromagnetic induction phenomenon on the cell membrane has gradually emerged.A new three-dimensional memristive HR neuron model is reported in [15].The HR neuron model can show a discrete memristor initial offset boosting behavior, with an infinite number of discontinuous attractors.By changing the initial value of memristor, the generated firing patterns can be boosted to different discrete levels.In [16], the memristive HR neuron model operates in hidden bursting firing patterns when using an electromagnetic induction current generated by the threshold memristor to replace the external current.The application of locally active memristor in a neuron model has gradually emerged [17][18][19][20].Ref. [21] proposed a new type of locally active and non-volatile memristor, and it was used as a connecting autapse, so that both the firing pattern and multistability can be observed.Ref. [22] proposed a new neuron model based on a simple structure and excellent characteristics with a locally active memristor, which can exhibit more complex firing behaviors within the locally active region.Ref. [23] coupled 2D HR neurons to a 3D Hopfield neural network using a memristor as a synapse, in which multiple firing patterns and multistability were observed.Additionally, multistable coexistence is common in chaotic systems, and it needs to be eliminated in certain situations [24,25].For this reason, by introducing a memristor into the HR model, an interesting phenomenon of multistability disappearance can be found when the memristor coupling strength changes, which is rarely reported in previous studies.
For various changes in the internal and external environments of the body, neurons encode, transmit, and decode the information in different firing patterns (the presence or absence of action potential, frequency, peak, etc.).Different external stimuli cause different neural firing patterns and lead to different physiological effects [26][27][28][29].In fact, a neuronal cell can be thought of as a signal processor that can receive multiple signal inputs at the same time.Strong stimulation can input enough energy to induce mode transition, and the angular frequency can slightly modulate the firing rhythm at a certain intensity.Multiple signal inputs can be equivalent to coherence resonance [30,31].In addition, the high-low frequency signal is widely used in neuron models.In [32], the mode transition of electrical activities in neurons can be investigated when an electrical stimulation with a high-low frequency current is injected into a new HR neuron model.As discussed in [33], an improved HR neuron model with the periodic high and low frequency electromagnetic radiation and the Gaussian white noise is constructed to explore the modes transition in the electrical activities of neurons.Based on the above ideas, this paper proposed a memristive HR model with multi-frequency AC injection.Under the effect of both the single and double excitation signals, we find that the improved model possesses a property of the coexistence of different firing patterns.In addition, the multistability of the model will be eliminated when the memristor coupling strength k changes.
The rest of this paper is organized as follows.An improved memristive HR neuron model with multi-frequency AC injection is constructed in Section 2, and the stability of the equilibrium points is studied.Section 3 discusses the dynamical behaviors under single excitation, and the influence of the memristor coupling strength is also illustrated.Different firing patterns under high-low frequency current are determined in Section 4. The analog circuit is designed in Section 5, and a summary of the full paper is given in Section 6.

Model Description
In recent years, memristors have developed rapidly in the fields of neurons and neural network dynamics.Memristor neuron network dynamics is a new interdisciplinary subject.When neuron membrane potential is affected by the external magnetic field, induced current will react on the neuron membrane potential.This kind of influence can be described by a flux-controlled memristor [34].
The HR model is a neuron model [4] described by a set of ordinary differential equations with cubic nonlinear terms, and it can be described as where x and y denote the membrane potential, and the recovery variable of the neuron, respectively.a, b, c, and d are positive parameters.I represents the external input current.
A simple generalized flux-controlled memristor model is described as where v and i are the input voltage and output currents of the memristor, respectively, ϕ is the magnetic flux indicated as the state variable, W(ϕ) is the flux-controlled conductance, and α and β are constant.In this paper, α and β are set to 0 and 0.01, respectively.By introducing the flux-controlled memristor and injecting two external alternating currents (ACs), an improved memristive HR neuron model with multi-frequency AC injection was constructed as follows where and k is the memristor coupling strength.x and y denote the membrane potential and the recovery variable of the neuron, respectively.
The equilibrium point of the proposed improved memristive HR neuron model can be obtained by setting the left side of the Equation (3) as equal to zero From 1 − 5x 2 − y = 0 and x − ϕ = 0, we know that y = 1 − 5x 2 , x = ϕ.Next, let x = η and the equilibrium point be set to E. Thus, the equilibrium point can be given by Then, we can substitute E, α = 0, and β = 0.01 into y + 3x 2 − x 3 + kW(ϕ)x + I 1 +I 2 = 0, and the value of η is determined by the real root of the equation According to the Cardan discriminant, if there is an equation of f (x) = ax 3 + bx 2 + cx + d = 0 and a = 0, then p, q, and ∆ can be defined as Then, for the AC-induced HR neuron model proposed in this paper, Generally, there are three solutions for Equation (6), which can be given by The discussion is divided into three situations.If ∆ > 0, there is only one real root foot for Equation (6).If ∆ < 0, there are three real root feet for Equation (6).If ∆ = 0, there are two real root feet for Equation ( 6) [31].
The Jacobian matrix at the equilibrium point E can be deduced as The eigenvalues can be obtained by solving the characteristic polynomial where The stability of the equilibrium point can be judged by the magnitude of λ.When the real parts of all eigenvalues λ are negative, the equilibrium point of the system is stable, and the solution of the system tends to the equilibrium point.When the real parts of all eigenvalues λ are positive, the equilibrium point is unstable and the solution of the system is far from the equilibrium point.When the real part of the eigenvalue λ is positive or negative, the equilibrium point is called the saddle point and is unstable.When k is set to 1, A 1 and A 2 are set to 1, and f 1 and f 2 are set to 0.05; Figure 1 shows the trajectory of the parameter η, the real part of the eigenvalues λ 1 , λ 2 , and the partial real part of the eigenvalues λ 3 changing with time.It is obvious that, with the evolution of time, λ 1 jumps between positive and negative, resulting in the changing stability of the equilibrium.The stability of the equilibrium under these parameters during the whole time interval (5,25) of a period of I 1 + I 2 is discussed in Table 1 in detail.The type of equilibrium will change over time, which shows that the injected AC in the HR model can generate complex dynamic behaviors.

AC-Induced Complex Dynamical Behaviors under Single Excitation
In this section, when a single excitation I 1 is considered and I 2 is set to 0, the proposed model can exhibit complex dynamic behaviors which can be discussed in different aspects.It is worth mentioning that the MATLAB ODE45 algorithm is used to draw the phase diagrams, maximum Lyapunov exponent, bifurcation diagrams, and time evolutions.
Phase diagrams, bifurcation diagrams, and Lyapunov exponent diagrams are important tools in the study of chaos.The phase space structure of a chaotic system can be observed in the phase diagram, which can help acquire an in-depth understanding of the dynamic behavior of the system.In addition, the structure of the chaotic system is unstable, and any small appropriate disturbance will cause a sudden change in the topological structure of the system.This sudden change is called bifurcation, and we can observe the bifurcation situation of the system from the bifurcation diagram.The positive Lyapunov exponent is an important characteristic of a chaotic system.When the Lyapunov exponent is greater than 0, the system is chaotic.

Coexisting Asymmetric Firings When A 1 Changes
Setting the parameter k as 1, and f 1 as 0.5, take A 1 as a bifurcation parameter within the range of (0, 8).The bifurcation diagram with respect to A 1 is numerically simulated by the MATLAB ODE45 algorithm in Figure 2a, in which (5, 0, 0) and (−5, 0, 0) are the initial values of the red and blue trajectories, respectively.And, the corresponding maximum Lyapunov exponent graph is shown in Figure 2b, which is consistent with the bifurcation diagram.As can be seen from Figure 2, the presented neuron model (3) shows a reverse period-doubling bifurcation under the initial condition of (−5, 0, 0), while the neuron model is always in a stable periodic state under the initial condition of (5, 0, 0).Choose several representative values of A 1 to observe the coexisting behaviors of asymmetric firings.Then, the phase plane orbits corresponding to different values of A 1 are shown in Figure 3.When A 1 = 0.1, the time sequence exhibits a coexistence of fixed-point and chaotic spiking firing under the symmetric initial conditions of (−5, 0, 0) and (5, 0, 0).The corresponding phase diagram and time series are drawn in blue and red, which are shown in Figure 3a and Figure 3b, respectively.When A 1 = 3, as shown in Figure 3c,d, the coexistence of chaotic bursting firing with multiple spikes and periodic limit cycles is found in this improved HR model.It can be seen from Figure 3e,f that, when A 1 = 5.5, the coexistence of the periodic spiking firing and periodic limit cycles is found.

Coexisting Firing Patterns When f 1 Changes
Set the parameter k as 1, A 1 as 3, take f 1 as a bifurcation parameter within the range of (0.01, 0.1), and the initial conditions are (−5, 0, 0).Then, the bifurcation diagram and maximum Lyapunov exponent are depicted in Figure 4a and Figure 4b, respectively.It can be seen from Figure 4 that the improved HR model has rich dynamic behaviors including chaos, period, periodic windows, crisis scenario, and so on.The maximum Lyapunov exponent is greater than zero.which means that the model is in a chaotic state.Obviously, the chaotic region is small in this range, and there is no coexistence phenomenon in this range.Then, a larger range is chosen, and when f 1 is in the range of (0.1, 1), the coexistence phenomenon appears in the partial region of (0.438, 1), and the bifurcation diagrams under different initial conditions are depicted in Figure 5a.In contrast to Figure 4, there is a larger chaotic range in Figure 5, and the maximum Lyapunov exponent spectra shown in Figure 5b is consistent with the bifurcation diagram.

Influence of Coupling Strength k on Dynamics
As we all know, multistability depends on the initial condition of the systems, and is a common phenomenon in dynamical systems.However, sometimes it is necessary to avoid multistability when designing commercial devices with certain special characteristics, which often brings an inconvenience to the design.In this study, through numerical simulation, it can be found that the coexistence phenomenon can be eliminated when k is less than −6.9.This is an interesting phenomenon which is rarely reported before.Details are discussed below.
When A 1 is equal to 3 and f 1 is equal to 0.5, set k = −7, for which the phase diagrams and corresponding time series are shown in Figure 6, where the initial condition of the blue track is (−5, 0, 0) and the initial condition of the red track is (5, 0, 0).It is evident that the blue trajectory is in a state of chaotic bursting fire during the whole time range, while the red trajectory enters into chaotic bursting fire from the periodic orbit when t is approximately equal to 270.Similarly, when k = −10, and f 1 = 0.5, A 1 = 0.1, as shown in Figure 7, the blue trajectory is always in a state of spiking while the red trajectory enters into the spiking state after t = 200.It can be seen from these two sets of parameters that the memristor coupling strength k can eliminate the coexistence phenomenon with the time evolution when k is less than −6.9 and the proposed neuron model will eventually enter into a stable state under any initial conditions.This phenomenon makes a great contribution to multistability control.

Different Firing Patterns Are Driven by High-Low-Frequency Current
The application of two frequency signals is rife in neural systems, and it plays a significant role in the biological field, such as that the neural cell behaviors under high-lowfrequency ultrasonic irradiation can be examined.Then, the dynamical behaviors under high-low-frequency current are worth discussing.In this section, f 1 is regarded as a highfrequency electromagnetic radiation, and f 2 is regarded as a low-frequency electromagnetic radiation.When k = 1, A 1 = A 2 = 3, and f 1 = 0.5, the initial condition is (−5, 0, 0), with the increase in low frequency f 2 , the electrical activities of the HR neuron are investigated in Table 2.When f 2 = 0.002, the model presents a periodic bursting firing pattern with chaotic spiking per bursting, and when f 2 = 0.02, there is a periodic bursting firing pattern with three spikes per bursting.When f 2 = 0.04, a period-2 spiking pattern is found, and the model presents a chaotic spiking pattern when f 2 = 0.07.The three corresponding three Lyapunov exponents are also calculated in Table 2. Obviously, the quiescent state between the two bursts becomes shorter as the frequency increases.These different mode transitions can be achieved by adjusting the value of f 2 .

Circuit Implementation
To prove that multiple firing patterns of the HR neuron model under the multifrequency AC injection can be realized in hardware circuits, we carry out circuit implementation.In this section, the circuit schematic will be analyzed and designed.The amplifiers and the multipliers are selected as UA741 and AD633JN.They are both under the supply voltage ±15 V. And, the gain of the multiplier is 0.1.In order to complete the circuit implementation, the following transformation is taken out.Let τ = τ 0 t where the integral constant τ 0 is the time-scaling factor and τ 0 = 1 RC = 10,000.Set C = 50 nF, then R can be solved as 2 kΩ.To verify the MATLAB simulation results, the circuit schematic was designed in Figure 8.When k = 1, the circuit equivalent equation derivation process is: 1.
Variable-scale reduction transformation.Since the range of the attractor does not exceed the dynamic range of ±13 V, variable-scale reduction transformation of variety is not required.

2.
Time-scale transformation: 3. Differential-integral conversion: Because the inverse addition proportional arithmetic unit is used in the circuit, the Equation ( 14) is normalized to: ϕ]dt (15) where the memductance function equivalent circuit is In order to make the values of resistances as integral as possible, set R 7 = 300 kΩ and R 11 = R 14 = 100 kΩ.Then, by comparing the corresponding coefficients in Equation (15) with Equation ( 3), all the values of the circuit parameters can be obtained in Table 3 in detail.Finally, the multisim simulation is constructed based on the component parameters above.When only the external current I 1 is connected to the circuit, and the amplitude A 1 is 3 and 5.5, respectively, the frequency f 1 is equal to 5 kHz.Then the corresponding phase diagrams in the oscilloscope are shown in Figure 9.When A 1 is equal to 3, a chaotic bursting firing pattern is obtained, and when A 1 is equal to 5.5, a periodic spiking firing pattern is obtained which is completely consistent with the numerical simulation in Figure 3.The circuit of two external currents is also constructed.Set the amplitudes A 1 = A 2 = 3 and f 1 as 5 kHz, then adjust the value of f 2 to 20 Hz, 200 Hz, 400 Hz, and 700 Hz.The time series in the oscilloscope are shown in Figure 10.These are a periodic bursting fire pattern with chaotic spiking per bursting; a periodic bursting fire pattern with periodic-3 spiking per bursting; a periodic-2 spiking firing pattern; and a chaotic spiking firing pattern.The numerical simulation results are verified accurately through analog circuit simulation.

Conclusions
In this paper, we constructed an improved memristive HR neuron model with multifrequency AC injection.Then, the equilibrium of the model is analyzed.When the model is under single excitation, it can produce the coexistence of different firing patterns.With the introduction of the memristor, it can be found that the multistability of the model is eliminated when the memristor coupling strength k is less than −6.9.Additionally, when the external current is a high-low frequency excitation, the system can exhibit different kinds of firing patterns under the changing low frequency.Finally, the circuit experiment is carried out to verify the validity of the numerical simulation results.

Figure 1 .
Figure 1.The η-values of the equilibrium point E and the corresponding eigenvalues λ 1 , λ 2 , λ 3 , with the time evolutions.

Figure 4 .
Figure 4. Bifurcation diagram and maximum Lyapunov exponent spectra with respect to f 1 in the range of (0.01, 0.1): (a) Bifurcation diagram; and (b) Maximum Lyapunov exponent.

Figure 8 .
Figure 8. Circuit implementation of the improved memristive HR neuron model.

Author Contributions:
Conceptualization, M.W. and J.D.; methodology, M.W. and B.D.; writing-original draft, M.W. and J.D.; investigation, M.W., B.D. and J.D.; writing-review and editing, S.H., H.H.-C.I. and J.D.; supervision, M.W.; project administration, M.W.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the (National Natural Science Foundation of China) grant number (62071411), and the (Research Foundation of Education Department of Hunan Province) grant number [20B567]).Data Availability Statement: Data are contained within the article.

Table 2 .
Different firing patterns under low frequency f 2 .