Firing mechanism based on single memristive neuron and double memristive coupled neurons

Memristive neurons and memristive neural networks constructed based on memristors have important research significance for revealing the mystery of the brain. This paper proposes a compound hyperbolic tangent cubic nonlinear memristor, which has nonvolatile memory characteristics and local activity characteristics. In particular, the memristor also has three stable pinched hysteresis loops under different initial values. The memristor is applied to Fitzhugh–Nagumo neuron and Hindmarsh–Rose neuron to establish five different memristive neural models, and a series of firing dynamics analysis are carried out on them. At the same time, we not only discuss multiple firing patterns on a single memristive neuron and double memristive coupled neurons, but also compare which neuron and which coupled neural network the proposed memristor is more suitable for, which is a lack of comprehensive investigation in the published research. Furthermore, digital circuit experiment is performed on the FPGA development board to verify the firing mechanism of these memristive neural models, which has potential application value for some practical projects.

According to the latest reports, it is aware that the dynamical investigations based on the memristive neuron model and the memristive neural network model are the most common, which mainly include two research aspects. On the one hand, various neuronal electrical activities are simulated by applying memristors on an individual neuron. For instance, Lv and Ma [32] emphasized that memristive electromagnetic radiation plays an important role in the mode transition of neuronal electrical activity in the threedimensional Hindmarsh-Rose neuron model. Bao et al. [33] modeled a simple memristor neuron model without equilibrium point, which can show the hidden coexistence behavior of asymmetric attractors. Kapetanovic et al. [34] demonstrated that memristive electromagnetic induction effect leads to changes in action potential dynamics of Hodgkin-Huxley neuron. Lin et al. [35] introduced a nonvolatile memory and local activity memristor with three steady-state pinched hysteresis loops, which can be used as the neuronal autapse to build a memristive HR neuron model. Kafraj et al. [36] describes the firing behavior of the improved Izhikevich neural network model under memristive electromagnetic induction and noise, where external memristive electromagnetic induction will increase the firing rate of the neuron. Zhang et al. [37] presented a non-equilibrium Hindmarsh-Rose neuron model with memristive electromagnetic radiation effect, which can generate multiscroll hidden attractors with complex topology. On the other hand, the neural network models composed of neurons coupled by memristor have more novel and intricate firing mechanisms. For example, Li et al. [38] constructed a locally active memristor coupled FitzHugh-Nagumo neuron model and calculated the Hamiltonian energy function, which can detect periodic burst and multiscroll chaotic burst under the excitation of multi-level pulse current. Tan et al. [39] employed a local active memristor to couple two HR neurons, which can produce complicated firing behaviors in the local active region. Xu et al. [40] observed extreme multistability with different initial conditions in the continuous non-autonomous memristive Rulkov model. Njitacke et al. [41] explained the dynamical influences of asymmetric memristive electromagnetic induction and electrical synaptic coupling on two neurons. Chen et al. [42] discovered the multistability phenomenon in the Hopfield neural network composed of the memristor synapse-coupled two neurons, and conducted a detailed stability analysis for the line equilibrium. Bao et al. [43] detected the coexisting and synchronous discharge behaviors induced by initial conditions in the memristive synapse-coupled Morris-Lecar double neurons. Li et al. [44] studied the firing phenomenon caused by the memristor coupling between different neurons, that is, the bistable local active memristor was used to couple the FN neuron model and the HR neuron model bidirectionally. This kind of interesting coupling method can bring great inspiration to readers.
There is no doubt that hardware implementation of these neuron models and neural network models is indispensable in practical engineering applications. The hardware implementation methods can be divided into analog circuit implementation and digital circuit implementation. At present, a large number of literatures have been published on the neuron and neural network models realized by analog circuits. Cai et al. [45] proposed a novel smooth nonlinear fitting scheme to realize HR neuron model, which greatly reduces the experimental cost of analog circuit and is conducive to the hardware implementation of large-scale neural network. Bao et al. [46] implemented the hyperbolictype memristive Hopfield neural network on breadboard and gave the same results as the numerical simulations. Chen et al. [47] verified the bursting dynamics characteristics in nonautonomous memristor Fitzhugh-Nagumo circuit through analog hardware experiment. On the contrary, it is relatively scarce for digital circuits to achieve these neuron models and neural network models. Zhu et al. [48] proved that bistable dynamics in Tabu learning neuron model by performing digital circuit experiments on FPGA. Tlelo-Cuautle et al. [49] utilized FPGA to realize Hopfield and Hindmarsh-Rose neurons, and applied them to pseudo-random number generator and image encryption. Yang et al. [50] established a Hopfield neural network system with complex chaotic phenomena and designed its digital circuit based on DSP platform. To sum up, these neuron models and neural network models have high requirements for the setting of initial conditions, and the implementation of digital circuit can accurately set the required initial conditions, so the experimental results are more ideal than those of analog circuit. Inspired by the above review, this paper designs a new nonvolatile memory and local activity memristor, whose pinched hysteresis loop has three stable states under different initial values. Based on the designed memristor, FN neuron model and HR neuron model, five memristive neural models are proposed, which are the memristive FN neuron model, the memristive HR neuron model, the memristive coupled FN-FN neurons model, the memristive coupled FN-HR neurons model and the memristive coupled HR-HR neurons model. The firing mechanisms of these models are orderly analyzed via phase portraits, timing waveform diagrams, bifurcation diagrams, Lyapunov exponents and attraction basins. At the same time, the FPGA implementation further demonstrated the results of the numerical simulation.

A locally active memristor model
In this section, we introduce a new nonvolatile memory and local activity memristor with tristable state, which is missing from existing memristor models. The pinched hysteresis loops of the memristor are observed at different amplitudes, frequencies and initial values. The nonvolatile memory and local activity of the memristor are demonstrated by the power-off plot and the DC V − I plot. And the advantages of the designed memristor are highlighted by comparing with other memristors of the same type.

Model description
According to Chua's memristor principle [51], some memristor models can show nonvolatile memory characteristics and local active characteristics, which are generally found in some memristor models containing specific functions. In this study, we use a composite cubic hyperbolic tangent function nonlinear term to construct a tristable memristor model with nonvolatile memory and local activity, which is expressed as where γ, δ and ε are three positive memristor parameters that can be adjusted independently.
As shown in Fig. 1, the pinched hysteresis loop of the memristor model is numerically simulated by applying a sinusoidal voltage excitation source v m = A sin(2π Ft) and maintaining memristor parameters γ = 1, δ = 0.5 and ε = 1. Obviously, the pinched hysteresis loop of the memristor model can describe three essential characteristics of the memristor. When the initial value φ(0) = 0 and frequency F = 1 are fixed and the amplitudes A = 0.6, A = 1 and A = 1.5 are set, respectively, the side lobe area of pinched hysteresis loop increases with the increase in amplitude. When the initial value φ(0) = 0 and amplitude A = 1 are fixed and the frequencies F = 0.6, F = 1 and F = 1.5 are set, respectively, the pinched hysteresis loop shrinks at the origin, and its side lobe area decreases monotonically with the increase of the frequency, and when the frequency tends to infinity, the pinched hysteresis loop shrinks into a single-valued function.
In addition, coexisting pinched hysteresis loops can be observed in the proposed locally active memristor model under appropriate amplitude, frequency and initial value, as shown in Fig

Nonvolatile memory
As pointed out in [51], not all memristors are provided with nonvolatile memory, which needs to be judged by observing the power-off plot of memristors. If there are two or more negative slopes on the power-off plot of the memristor, it means that the memristor is a nonvolatile memristor. When the power is turned off, that is to say, v m = 0, and the state equation of the memristor is simplified tȯ From Eq. (2), the power-off plot of the memristor can be obtained as shown in Fig. 3. Its dynamical path is a smooth curve, and the additional arrow indicates the evolution direction of the state variable φ. Since the dynamical path of the nonvolatile memristor follows the no backtracking rule, for any initial state φ(0) on the curve, ifφ > 0, it must evolve to the right along the curve, otherwiseφ < 0, it must evolve to the left along the curve. Moreover, wheṅ φ = 0, the memristor has five equilibrium points and Q 4 (−2, 0) respectively, in which the three equilibrium points Q 0 (0, 0), Q 3 (2, 0) and Q 4 (−2, 0) have negative slopes and are asymptotically stable, while the two equilibrium points Q 1 (0.7237, 0) and Q 2 (−0.7237, 0) are unstable. Hence, the three stable equilibrium states of the nonvolatile memristor under different initial values can be obtained as Further, the corresponding stable memductance can be worked out It can be clearly seen from Eqs. (3) and (4) that the statė φ is diverse from the initial value φ(0). The above analysis results just explain why the nonvolatile memristor has tristable pinched hysteresis loops at different initial values. Of course, not all nonvolatile memristors possess local activity [51]. It is usually necessary to determine whether the memristor has local active characteristics with the help of DC V − I plot. Unlike the simulation of power-off plot of memristor, the simulation of DC V − I plot of memristor needs to setφ = 0 to calculate the equilibrium equation, as shown below where V is DC voltage, I is DC current and Φ is a variable equilibrium state. From Eq. (5), the DC V − I plot of the memristor can be drawn as shown in Fig. 4 , which are consistent with the corresponding five segment curves in Fig. 3. As shown in Fig. 4, the blue curve segment, the yellow curve segment and the green curve segment all have negative slopes. According to the local activity theory, the memristor can be considered as a local active memristor. Furthermore, some proposed local active memristors are investigated and summarized in Table 1. Unfortunately, most of the local active memristors reported in Table 1 are bistable, while the tristable local active memristors have received little attention.

Memristive neuron model
In this section, the same locally active memristor is connected to the two-dimensional FN neuron model and the two-dimensional HR neuron model, respectively, so that two three-dimensional locally active memristive neuron models with electromagnetic induction effect can be established. The dynamics of neuronal electrical activity of these two new locally active memristive neuron models are discussed through appropriate parameters and initial conditions. In contrast, the proposed local active memristor is more suitable for HR neuron model, which can stimulate it to produce different firing patterns.

Memristive FN neuron model
The two-dimensional FN neuron model proposed by Fitzhugh and Nagumo can be expressed as where x is the membrane potential, y is the ion current, I is the external stimulus, and a, b, c and d are the control parameters [52,53]. For this two-dimensional FN neuron model (6), we introduce a local active memristor in its first equation, and select appropriate model parameters a = where the new term ρ 1 φx denotes the induced current applied externally to FN neurons, ρ 1 indicates the electromagnetic induction intensity and the new variable φ represents the magnetic flux.
When the memristor parameters are respectively selected as ρ 1 = 0.1, γ = 1, δ = 0.5 and ε = 0.2, the three-dimensional locally active FN neuron model displays the periodic firing pattern in Fig. 5a, b, in which the green and red trajectories are respectively derived from the initial conditions (0, 0, 2) and (0, 0, −2). Obviously, the amplitude of periodic firing attractor and corresponding time-domain waveform will change with the different initial conditions. In addition, under different initial conditions (0, 0, 2) and (0, 0, −2), the bifurcation diagram and the first two Lyapunov exponents (calculated using Wolf's algorithm [54], same below.) with varying parameter ρ 1 are simulated in

Memristive HR neuron model
The two-dimensional HR neuron model proposed by Hindmarsh and Rose can be described as where x is the membrane potential, y is the recovery variable, I is the external stimulus, and a, b, c and d are the control parameters [55].
In the same way, for this two-dimensional HR neuron model (8), we also apply a local active memristor to its first equation, and choose suitable model parameters a = 1, b = 3, c = 1, d = 5 and external stimuli I = 0. Therefore, a three-dimensional local active HR neuron model with electromagnetic induction effect can where the new term ρ 2 φx stands for the induced current added externally to HR neurons, ρ 2 signifies the electromagnetic induction intensity and the new variable φ means the magnetic flux.
When the memristor parameters are respectively chosen as ρ 2 = 0.9, γ = 0.1, δ = 0.1 and ε = 0.5, the electrical activity of the three-dimensional locally active HR neuron model is investigated in Fig.  6a, b. As shown in Fig. 6a, b, the model exhibits chaotic firing pattern with the initial conditions (0, 0, 2) and is marked with green trajectory, while the model reveals periodic firing pattern with the initial conditions (0, 0, −2) and is marked with red trajectory. Accordingly, the change of electrical activity in the parameter range ρ 2 ∈ [0, 2] of three-dimensional locally active HR neuron model is analyzed by bifurcation diagrams and Lyapunov exponents, as shown in Fig. 6c, d. Evidently, the bifurcation diagrams and Lyapunov exponents are completely consistent under the two sets of initial conditions (0, 0, 2) and (0, 0, −2), but there is no significant difference in the change of electrical activity.

Memristor-coupled neurons model
In this section, two FN neurons, two HR neurons, and FN neuron and HR neuron are bidirectionally cou-pled through the local active memristor. The electromagnetic induced current generated by the local active memristor can describe the dynamical influence of the membrane potential difference between two neurons. The firing mechanisms of the three new memristive coupled neural network models are respectively studied by numerical simulation tools. The final results show that the firing patterns of HR neuron model coupled by the local active memristor are more abundant and interesting.

Memristive coupled FN-FN neurons model
Firstly, we utilize the local active memristor to couple two FN neurons in two directions to create a neural network model, which can be described as By setting external stimulus I 1 = 0.38, I 2 = 0.45 and memristor parameters ρ 3 = 0.03, γ = 2, δ = 1 and ε = 1, the firing behaviors of memristor coupled FN-FN neurons model (10) with two different initial conditions (0, 0, 0, 0, 0) (green represents chaos) and (−1, 0, 1, 0, 0) (red represents periodic) are investigated in Fig. 7. In particular, the occurrence of early after/depolarization (EAD) prearrhythmia behavior, a phenomenon of heart plateau spiking associated with arrhythmias, is shown in Fig. 7b. EAD and triggering activity in the autonomic nervous system were studied by Charpentier et al. [59] as early as 1993. In 2018, Bortolotto et al. [60] first reported the presence of cardiac spikes in the HR model. It should be noted that when the initial conditions are (0, 0, 0, 0, 0) and the memristor parameters remain unchanged, the two FN neurons adding the same external stimulus I 1 = I 2 = 0.38 can present phase synchronization and complete overlap of the corresponding timing waveforms, as shown in Fig.  8. At the same time, the bifurcation diagrams and Lyapunov exponents of the memristor coupled FN neural network model with parameter ρ 3 are respectively plotted in Fig. 9. It can be seen that the adjustable range of firing behavior dependent on parameter ρ 3 is relatively small, less than 0.05.

Memristive coupled FN-HR neurons model
Secondly, we adopt the local active memristor to interconnect one FN neuron and one HR neuron in two ways to compose a neural network model, which can be expressed as The firing mechanisms of the memristor coupled neurons model (11) composed of the FN neuron and HR neuron are discussed by using the memristor parameter γ = 1, δ = 0.5, ε = 0.5 and the coupling strength ρ 4 as the main firing mode control parameter. Considering that the initial conditions are (0, 0, 0, 0, 0), the phase portraits and corresponding timing waveforms of the model with three different coupling strengths are analyzed in Fig. 10. When the coupling strength ρ 4 is equal to 0.1, the model shows chaotic spike, when the coupling strength ρ 4 is equal to 0.5, the model displays chaotic burst, when the coupling strength ρ 4 is equal to 1, the model exhibits periodic burst. Moreover, we conducted in-depth research using bifurcation diagrams and Lyapunov exponents, as shown in Fig.  11, to clarify the dynamics of electrical activity of the memristor coupled neurons model. It is obvious that the spike firing and burst firing of the memristor coupled neurons model appear mode transition with the increase in coupling strength ρ 4 .

Memristive coupled HR-HR neurons model
Finally, we employ the local active memristor to connect two HR neurons bidirectionally to form a neural network model, which can be depicted as Under the initial conditions (−1, 0, 1, 0, 0) and the memristor parameters γ = 2, δ = 1, ε = 1, the firing behaviors of the memristor coupled HR-HR neurons model with the variation of coupling strength are examined in Fig. 12. When the coupling strength ρ 5 is 0.15, the chaotic firing and its corresponding time-domain waveform of the model are marked with green track. When the coupling strength ρ 5 is 0.25, the periodic firing and its corresponding time-domain waveform of the model are marked with red track. In addition, as shown in Fig. 13, when the initial conditions are (0, 0, 0, 0, 0) and the memristor parameters remain unchanged, the model can also achieve phase synchronization, while when the initial conditions are (−1, 0, 1, 0, 0) and the memristor parameters keep invariant, the model cannot be synchronized. Therefore, as long as the parameters of the two HR neuron models and the memristor models are adjusted to the same, its initial conditions will affect the synchronization phenomenon.
It is very interesting that the firing multistability are observed in the memristor coupled HR-HR neurons model, which is not found in the memristor coupled FN-FN neurons model and memristor coupled FN-HR neurons model. In order to further explore the firing multistability of the memristor coupled HR-HR neurons model, the diagrams with attractor changes for different initial conditions and the bifurcation diagram depending on the coupling strength are respectively drawn in Fig. 14. For example, when verifying the initial value x 1 (0) in the [−10, 10] region, other initial values should be set to 0 and the model parameters should be kept unchanged, as shown in Fig. 14a. In a similar way, the diagrams of changes under different initial conditions from Fig. 14b, e can be obtained one after another. For the bifurcation diagram with the change of coupling strength ρ 5 , namely Fig. 14f, the initial conditions need to be fixed as (−1, 0, 1, 0, 0) and other model parameters remain unchanged. In addition,     Fig. 16 under different initial conditions. It is obviously that the coexisting multiple attractors are exhibited in the memristor coupled HR-HR neurons model, that is, it generates the phenomenon of multistability.

Digital hardware experiment
In this section, the firing behaviors of the proposed five memristive neural models are effectively verified in digital hardware experiment. The process and note of FPGA implementation are described in detail, and the experimental results are displayed one by one on the oscilloscope. As far as the experimental results are concerned, it is not difficult to see that the results achieved by the digital circuit are quite ideal, which are highly consistent with the numerical simulation results.

FPGA implementation
For real biological neurons and nervous systems, it is vital to restore their functions through hardware experiments. As we all know, field programmable gate array is a reconfigurable electronic platform, which has the advantages of rich wiring resources, repeatable programming and high integration. It is very suitable to realize the hardware standard of these memristive neural models. Hence, we designed the digital hardware circuit of memristive neuron and memristive neural network based on FPGA on Vivado 2018.3. In this experiment, we chose the Xilinx xc7z100ffg900-2 chip as an validation tool for digital circuit realization, which has abundant programmable logic units inside and powerful functions. Then, according to the basic operation symbols contained in the equations of the five memristive neural models, such as addition, subtraction, multiplication, division and hyperbolic tangent function, to complete the design of each module and the required IP core by using hardware description language Verilog. And these IP cores and modules all follow the 32-bit IEEE754-1985 floating point standard. Furthermore, the classical fourth-order Runge-Kutta algorithm is used to numerically solve these equations, which is the most important part for FPGA implementation. The formula of the classic fourth-order Runge-Kutta algorithm is given as where where y n+1 is the next iteration value, y n is the current iteration value, h is the time interval, and k 1 , k 2 , k 3 and k 4 represent different slopes, respectively. And the hyperbolic tangent function contained in these equations can be approximately described by a simple piecewise function, which can be expressed as It should be noted that because the dimensions of five memristive neural models proposed in this work are different, including two three-dimensional memristive neuron models and three five-dimensional memristive neural network models, the program of the RK4 algorithm module needs to be modified accordingly. This Xilinx xc7z100ffg900-2 chip uses a differential clock which is slightly different from ordinary clocks in other chips. Set the number of iterations to 1000 and the sampling interval to 0.01, as well as assign the correct pins in the constraint file. Then, the compiled Verilog program is synthesized and implemented in Vivado 2018.3, and the generated bitstream file is downloaded to the Xilinx xc7z100ffg900-2 chip. Ultimately, it is connected to the oscilloscope through the digital to analog converter to observe the experimental results.  Fig. 19. Obviously, the experimental results indicate that the digital circuit implementation is a wonderful choice. The FPGA implementation can determine the precise initial conditions of these memristive neural models, so the obtained experimental results can highly reproduce the results of the numerical simulation, and the error cannot be distinguished by the naked eye.

Conclusion
In this work, we not only propose a new memristor model, but also apply it to FN neuron and HR neuron to construct five memristive neural models. The electrical activity dynamics of these memristive neuron models and memristive coupled neural network models are studied numerically and experimentally. The results show that each model exhibits different shapes of attractors and its own unique firing mechanism. Even the same memristor may not be able to stimulate rich firing modes for each model. For instance, the memristor FN neuron model can only show a single periodic firing mode, while the memristor HR neuron model can show both periodic firing mode and chaotic firing mode. There is also firing multistability in the memristive coupled HR-HR neurons model, However, it is not found in the memristive coupled FN-FN neurons model and memristive coupled FN-HR neurons model. This full investigation provides readers with a new insight that will help to emancipate the firing mechanism of the brain. In future research work, we will propose more novel memristive neural models and explore their fir-ing dynamics, so as to contribute our meager efforts to the development of neuroscience and neuromedicine. CX20200884, and by the Postgraduate Scientific Research Innovation project of Changsha University of Science and Technology under Grants CX2021SS72 and CX2021SS69.

Declaration
Conflict of interest The authors declare that they have no conflict of interest. The authors have no relevant financial or nonfinancial interests to disclose. All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by [Hui Shen], [Fei Yu] and [Shou Cai].