Hidden firing patterns and memristor initial condition-offset boosting behavior in a memristive Hindmarsh-Rose neuron model

Electromagnetic induction can effectively induce abundant firing patterns in neurons. In modeling a neuron model with the electromagnetic induction effect, an electromagnetic induction current is frequently added to the state equation of membrane potential. To more properly reflect the non-uniform distribution of the ions inside and outside the neuron membrane, an ideal flux-controlled memristor with sinusoidal memductance function and non-linearly modulated input is raised to depict an electromagnetic induction effect on a Hindmarsh–Rose neuron model, and thereby, a three-dimensional (3D) memristive Hindmarsh–Rose (mHR) neuron model is built in this paper. The proposed mHR neuron model possesses no equilibrium point since the involvement of the ideal flux-controlled memristor, which induces the generation of hidden dynamics. Numerical results declare that the mHR neuron model can generate abundant hidden dynamics, i.e., periodic spiking, chaotic spiking, period-doubling bifurcation route, tangent bifurcation, and chaos crisis. These hidden dynamics are much related to the memristor coupling strength and externally applied stimulus. Afterward, the memristor initial condition-offset boosting behavior is revealed. This can trigger the generation of infinite multiple coexisting firing patterns along the memristor variable coordinate. These coexisting firing patterns have identical attractor topology but different locations in the phase plane. Finally, an analog circuit is designed for implementing the mHR neuron model, and PSIM-based circuit simulation is executed. The circuit-simulated results perfectly verify the generation of hidden infinite multiple coexisting initial condition-offset boosting firing patterns in the proposed mHR neuron model.


Introduction
The nervous system contains a huge number of biological neurons, which are the basic information handling and integrating units of a biological nervous system [1]. The dynamical properties of these biological neurons are crucial for determining the behaviors of the nervous systems [2,3]. Thus, modeling of the biological neuron and exploring its dynamical behaviors are research hotspots and attract many researchers' attention. Up to date, numerous neuron models have been constructed to depict different kinds of biological neurons, and they can roughly be divided into two categories, i.e., the continuous-time neuron model [4][5][6][7][8][9][10] and discrete-time map [11][12][13]. In the literature, some of the continuous-time neuron models were built based on the electrophysiological ion transport mechanism [4][5][6][7]. In addition, some continuous-time neuron models [8][9][10] and discrete-time maps [11,13] were built on dynamical assumptions to reproduce electrical activities without regard to the neuron structure [14]. No matter which category they are, all these neuron models can effectively reproduce abundant firing patterns in response to the change of the electrophysiological environment. Recently, the electromagnetic induction effect has immersed great scientists' concern, which can greatly affect the neuron dynamics [15][16][17] and neural network behaviors [18,19].
Actually, the media of a biological neuron can be magnetized during its polarizing and depolarizing processes [20]. On one hand, periodic firing can control the transition, pumping, and distribution of calcium, potassium, chloride, etc., ions in these processes. On the other hand, the distribution, pumping, and transition of these ions can induce fluctuation of the membrane potential. Meanwhile, the electromagnetic induction effect is induced when the ions pass through the neuron membrane. Thus, ion channel currents and the electromagnetic induction current simultaneously affect the membrane potential. In the literature, flux-controlled memristors were used in various neuron models to depict the dynamic relation between the membrane potential and magnetic flux [21]. In other words, the flux-controlled memristors were used in the Hodgkin-Huxley neuron model [22][23][24], Izhikevich neuron model [25], FitzHugh-Nagumo neuron model [26], and threedimensional (3D) [27][28][29][30][31]/two dimensional (2D) [32][33][34] Hindmarsh-Rose neuron model to depict the electromagnetic induction effect. These memristive neuron models can generate abundant firing patterns since the involvement of the flux-controlled memristor. To explain in detail, these flux-controlled memristors are non-ideal with a quadratic polynomial memductance function [22-30, 32, 33] and ideal with hypertangent/sinusoidal/cosinoidal memductance functions [20,31,33]. It is worth noting that the state equations of these memristors are linearly controlled by a membrane potential. Actually, the membrane potential possesses non-linear regulation on magnetization since the non-uniform distribution of the ions inside and outside the neuron membrane. To stress this issue, an ideal flux-controlled memristor with a sinusoidal memductance function and non-linear modulation on a memristor magnetic flux is raised to availably depict the electromagnetic induction effect in this paper.
The Hindmarsh-Rose neuron model is a simple kind of neuron model built on dynamical assumptions, which can reproduce main firing patterns of the biological neuron [9]. In the literature, the memristive Hindmarsh-Rose (mHR) neuron models with ideal memristors can easily generate multistability with coexisting firing patterns. In this case, the mHR neuron model has no equilibrium point, which induces the occurrence of hidden dynamics [20]. In particular, the initial condition-offset boosting behavior is triggered since the involved ideal flux-controlled memristor possessing sinusoidal/cosinoidal memductance functions [20,34,35], which is very different from the parameter-offset boosting behavior [36][37][38]. This induces the occurrence of extreme multistability with infinite multiple coexisting firing patterns. These coexisting firing patterns own attractors having identical topology and boosting along the memristor variable coordinate [20,34]. Herein, a mHR neuron model with our proposed memristor is tamed for simplicity but without losing generality. The hidden dynamics and initial condition-offset boosting behavior of the mHR neuron model are investigated by numerical simulation and PSIM-based circuit simulation in this paper. A brief comparison between some aforementioned mHR neuron models and the model reported in this paper is demonstrated in Table 1. It is demonstrated that the memristor employed in building the mHR neuron model in this paper is different from the memristors reported in the aforementioned literature works. The electromagnetic induction effect characterized by the memristor is established by considering the periodic magnetization processing and nonuniform distribution of the ions inside and outside the neuron membrane.
The remainder of this paper is formulated as follows: Section 2 explains the building of a mHR neuron model with hidden dynamics. Section 3 explains memristor parameter-and stimulus parameter-related dynamical distributions and bifurcation behaviors by numerical simulations. Section 4 explains the memristor initial condition-offset boosting behavior and infinite multiple coexisting firing patterns. Section 5 explains the analog circuit design and PSIM-based circuit simulation. Finally, Section 6 briefly concludes the main results of this paper.

Memristive Hindmarsh-Rose neuron model
Considering the periodic magnetization process and non-linear regulation of electromagnetic induction, a memristor with sinusoidal memductance function and non-linearly modulated input is raised, which is mathematically expressed as follows: where v M and i M represent the terminal voltage and current, respectively. W(φ) = sin(φ) is the periodic memductance function, and magnetic flux φ is the memristor inner state variable. Different from the memristor reported in [20], this memristor has a non-linearly modulated input, i.e., a hypertangent function, to reflect the non-uniform distribution of the ions inside and outside the neuron membrane. The hypertangent function is continuously derivable and bounded above and below.
To investigate this kind of electromagnetic induction effect on a neuron, the memristor is introduced into the existing 2D Hindmarsh-Rose neuron model [9]; thereby, a 3D mHR neuron model is built as follows: where x is the membrane potential and y is the recovery variable. a, b, c, and d are four controllable parameters in the original model [39]. I is the externally applied stimulus, and k is the coupling strength of electromagnetic induction. We mainly consider the dynamical effect of the externally applied stimulus I and coupling strength k on the mHR neuron model in the following sections. Therefore, the four controllable parameters, namely, a = 1, b = 3, c = 1, and d = 5, are assigned as the original parameters in [39,40]; I and k are adjustable parameters with positive values, and their typical values are preset to I = 1.5 and k = 2. By setting the left sides of (2) equal to 0, one can obtain the following equation: Evidently, there is no solution of (3) since I has a positive value. In other words, model (2) has no equilibrium point. Therefore, the dynamical behaviors and firing patterns generated by the 3D mHR neuron model (2) are hidden [41].

Parameter-related hidden dynamics
In this section, we mainly focus our concern on the parameterrelated hidden dynamics with the two adjustable parameters of the coupling strength k and externally applied stimulus I. The initial conditions [x(0), y(0), φ(0)] = [0, 0, 0] are utilized. The MATLABbased ODE45 algorithm with a fixed time-step duration of 10 −2 s and time-end duration of 800 s is utilized to calculate the bifurcation diagram, and the Jacobi matrix-based Wolf's method with a timestep duration of 0.1 s and time-end duration of 4,000 s is employed to calculate Lyapunov exponent spectra [42].

Dynamical distribution
When the two adjustable parameters k and I are varied in 0.5 ≤ k ≤ 3 and 0 ≤ I ≤ 3, dynamical distributions of the bifurcation diagram and dynamical map in the k-I parameter plane are simulated, as shown in Figure 1. The 2D bifurcation diagram is depicted by checking the periodicities of the membrane potential x, as shown in Figure 1A, that is, the trajectories with different periodicities are painted by different colors. The red marked by CH represents chaos, and the other colors represent period-1 to period-8 marked by P1 to P8, respectively. One can see that the 2D bifurcation diagram possesses a ribbon structure in some regions and the ribbons marked by P1, P2, P4, and P8 appear in sequence.

Paper
Dimensionality Memristor state equation Implementation Frontiers in Physics frontiersin.org Also, numerous ribbons marked by P3, P5, and P7 are embedded in the CH (red) region or near the neighborhood of the CH region. These declare that numerous periodic windows generated via tangent bifurcations [42] and period-doubling bifurcations [29] are triggered by varying the two parameters. In addition, the 2D dynamical map described by the largest Lyapunov exponent (LLE) is employed to depict the parameter-related dynamical distribution in the k-I parameter plane, as shown in Figure 1B. The colorized domains are painted with different colors according to the values of LLE: red for chaos with positive LLE and other colors for a period with negative LLE. It is demonstrated that the dynamical behaviors depicted by dynamical distributions of the 2D bifurcation diagram and 2D dynamical map are completely identical. These numerical results demonstrate that the coupling strength and externally applied stimulus can induce abundant dynamical behaviors on neuron properties of the 3D mHR neuron model.

Bifurcation behavior
To more clearly demonstrate the bifurcation behaviors with the coupling strength k and externally applied stimulus I, the onedimensional (1D) bifurcation diagram and Lyapunov exponent spectra (LEs) are numerically simulated with the variations of k and I, as shown in Figures 2A, B, respectively. The representations at the top of Figures 2A, B display the 1D bifurcation diagrams of the membrane potential x, while the representations at the bottom exhibit LEs.
Herein, the externally applied stimulus I = 1.5 is fixed, and the coupling strength k is varied in 0.5 ≤ k ≤ 3. The 1D bifurcation diagram for the maximum value of the membrane potential x (marked as x max ) is depicted in Figure 2A. When k increases from 0.5, the trajectory of mHR neuron model (2) starts from period-1, then enters chaos via the forward perioddoubling bifurcation route [29], returns to period-2 via tangent bifurcation [42], and ends up to period-1. It is worth noting that the period-doubling bifurcation route demonstrates the transition of P1-P2-P4-P8-CH. In Figure 2A, only the first Lyapunov exponent LE 1 and partial second Lyapunov exponent LE 2 are shown for better visualization since the third exponent is very small. The LE 1 exponent has a zero value for periodic states with different periodicities and a positive value for the chaotic state. The LE 2 exponent increases to zero and immediately returns to the negative value along with the occurrence of period-doubling bifurcations. It is observed that the bifurcation behaviors revealed by the 1D bifurcation diagram (up) are effectively verified by the LEs (bottom) in Figure 2A.
Then, we fix coupling strength k = 2 and change the externally applied stimulus I in 0 ≤ I ≤ 3. In Figure 2B, one can see that with the increase of I, the mHR neuron model undergoes period-1 to chaos via chaos crisis [43], to period-2 via tangent bifurcation, to chaos via forward period-doubling bifurcations, to period-3 via reverse period-doubling bifurcations, to chaos via chaos crisis, and finally to period-1 via reverse period-doubling bifurcations. It is worth noting that there exists a small periodic window near the neighborhood of I = 1.95. The LE 1 exponent drops to zero and maintains in the narrow parameter range and then returns to a positive value with the appearance of a periodic window. Obviously, the evolution of LEs confirms the occurrence of bifurcation behaviors.
Consequently, these bifurcation behaviors can lead to the occurrence of abundant periodic firing patterns with different periodicities and a chaotic firing pattern in the 3D mHR neuron model and can be regulated by the memristor coupling strength and externally applied stimulus.

Hidden firing patterns
In this section, five values of coupling strength k with I = 1.5 are selected from Figure 2A to partially display the firing patterns emerged from the 3D mHR neuron model, as shown in Figure 3. Time-domain waveforms of the membrane potential x (left) and corresponding phase portraits in the φ-x phase plane (right) are demonstrated; they are period-1 spiking, period-2 spiking, period-4 spiking, period-8 spiking, and chaotic spiking for k = 1, 1.5, 1.6, 1.65, and 2, respectively. These numerical results further demonstrate that the 3D mHR neuron model can generate hidden firing patterns of periodic spiking behaviors with different periodicities and a chaotic spiking behavior. Moreover, the state transition (P1-P2-P4-P8-CH) of firing patterns confirms the generation of the period-doubling bifurcation route.

Initial condition-related dynamics
Of particular interest, the 3D mHR neuron model (2) can show the initial condition-offset boosting behavior since the involvement of sinusoidal memductance function [20,44]. This can trigger infinite multiple coexisting firing patterns for a fixed set of model parameters. In this section, we mainly focus our attention on this issue. Herein, we only consider the memristor initial conditioninduced dynamical effect and set the initial conditions as [x(0), y(0), φ(0)] = [0, 0, φ(0)]. The numerical simulation settings are identical with those employed in Section 3.

Initial condition-offset boosting
Herein, two sets of model parameters, i.e., I = 1.5, k = 1.5 and I = 1.5, k = 2, are selected as paradigms to demonstrate the memristor initial condition-offset boosting behavior. The memristor initial condition is adjusted in the region [-20, 20]. The bifurcation plots of the 1D bifurcation diagram and Lyapunov exponent spectra for the two sets of model parameters are shown in Figures 4A, B, respectively.

Frontiers in Physics
frontiersin.org Figure 4A shows the bifurcation plots for I = 1.5 and k = 1.5, which demonstrates period-2 spiking firing patterns with the increase of the initial condition φ(0). One can see that the locations of these firing patterns possess a step-by-step structure and their dynamic amplitude is identical. In addition, the memristor initial condition applies at 2π initial condition-offset, which is the period of the sinusoidal memductance function sin(φ). In this case, the step change happens periodically. In addition, the LE 1 and LE 2 exponents demonstrate constant Lyapunov exponent spectra and are not related to the memristor initial condition, i.e., LE 1 = 0 and LE 2 = −0.15.
As shown in Figure 4B, the bifurcation plots for I = 1.5 and k = 2 are elaborated. The results also display the step-by-step memristor initial condition-offset boosting behavior; thereby, multiple coexisting chaotic firing patterns are generated. Also, the memristor initial condition applies at 2π initial condition-offset, and these steps periodically occur. Identically, the Lyapunov exponent spectra display constant values of LE 1 and LE 2 and possess LE 1 = 0.21 and LE 2 = 0 for the chaotic firing patterns. It is worth noting that the initial condition-offset boosting behavior exists for other firing patterns under different fixed model parameters. In addition, the memristor initial condition-offset boosting behavior can trigger infinite multiple coexisting firing patterns. We only demonstrate this behavior in a finite range of memristor initial conditions in this study.

Infinite multiple coexisting firing patterns
In the previous bifurcation analysis for the memristor initial condition-offset boosting behavior, the results display the occurrence of hidden infinite multiple coexisting firing patterns in the 3D mHR neuron model. To further demonstrate this striking memristor initial condition-offset boosting behavior, the coexisting firing patterns are displayed by the phase portrait in the φ-x phase plane, as shown in Figure 5. The coexisting period-2 spiking firing patterns for I = 1.5 and k = 1.5 under the initial conditions φ(0) = −18, −12, −6, 0, 6, 12, and 18 are shown in Figure 5A. These initial conditions are selected in each step from the bifurcation diagram in Figure 4A. The offsets among each of the two adjacent attractors of   Figure 5B. It is worth noting that we do not select the initial conditions with the interval 2π for their convenient setting in PSIM-based circuit simulation. These results demonstrate the generation of hidden infinite multiple coexisting firing patterns in our proposed 3D mHR neuron model.

Analog design and PSIM-based circuit simulation
The analog circuit design of neuron models is crucial for investigating the neuron dynamics and exploring neuron-based engineering applications [45,46]. The 3D mHR neuron model can be easily designed by utilizing passive circuit components of a capacitor and resistor and the integrated chips of operational amplifier, multiplier, and trigonometric circuit modules. The circuit schematic representation is well-designed and given in Figure 6. The memristor equivalent circuit contains a hypertangent circuit module [47], a sinusoidal function chip U 5 , an operational amplifier U 4 , a multiplier M 3 , a capacitor, and two resistors, as shown in the bottom part of Figure 6. The main circuit involves two integrators, two multipliers, and an inverter, as shown in the top part of Figure 6. Then, the circuit state equations can be correspondingly built as follows: where v x , v y , and v φ are three circuit state variables corresponding to the model variables x, y, and φ, respectively. We suppose C 1 = C 2 = C 3 = C and the integral time constant RC = 1 ms, i.e., C = 100 nF and R = 10 kΩ. In addition, the recovery variable v y is linearly transformed to reduce its dynamic amplitude as follows: The linear transformation is conducive to a hardware experiment since the value of the recovery variable approaches the saturation voltage of operational amplifiers without transformation. Thereby, the other circuit parameters are calculated as R 1 = 10 kΩ, R 2 = 1 kΩ, R 3 = 3.3 kΩ, R 4 = 10 kΩ, R 5 = 10 kΩ, R 6 = 20 kΩ, R 7 = 100 kΩ, R 8 = 10 kΩ, R 9 = 3.3 kΩ, g 1 = g 2 = g 3 = 1 V −1 , and V c = 1 V by comparing (2) with (4). It is worth noting that the two adjustable parameters can be regulated by I = V I and k = R/R 9 . The initial states of three capacitors are assigned as Employing the circuit schematic representation illustrated in Figure 6, a PSIM-based simulation circuit is built and circuit FIGURE 6 Circuit schematic representation of the 3D mHR neuron model. Frontiers in Physics frontiersin.org parameters are accurately set. First, the two adjustable circuit parameters are set to V I = 1.5 V and R 9 = 6.667 kΩ corresponding to I = 1.5 and k = 1.5. By respectively setting the memristor initial condition v φ (0) to −18 V, −12 V, −6 V, 0 V, 6 V, 12 V, and 18 V, PSIM-based circuit simulations are executed and phase trajectories are obtained in the v φ -v x phase plane, as shown in Figure 7A. Then, the two adjustable circuit parameters are set to V I = 1.5 V and R 9 = 5 kΩ corresponding to I = 1.5 and k = 2. The PSIMbased circuit simulations are illustrated in Figure 7B. One can see that the PSIM-based circuit simulations in Figure 7 are very consistent with the numerical results in Figure 5. These circuitsimulated results manifest the occurrence of hidden infinite multiple coexisting initial condition-offset boosting firing patterns in our proposed 3D mHR neuron model. It is worth noting that the power supplies for operational amplifiers and multipliers are, respectively, set to ±30 V and ±15 V in PSIM-based circuit simulation.

Conclusion
In this paper, an ideal flux-controlled memristor with sinusoidal memductance and non-linearly modulated input was presented to depict the electromagnetic induction effect in biological neurons. Then, the electromagnetic induction effect on an existing 2D Hindmarsh-Rose neuron model was elaborated. Theoretical analysis and numerical simulation demonstrated that the 3D mHR neuron model can generate the hidden memristor initial condition-offset boosting behavior with infinite multiple coexisting firing patterns. The attractors of these firing patterns have an offset along the memristor variable coordinate, and the offset is identical with the period of memductance function. In addition, PSIM-based circuit simulation further confirmed the validation of the analog circuit design and generation of the initial condition-offset boosting behavior. It is worth noting that the power supplies should be suitably set in PSIM-based circuit simulation to capture the offset boosting firing patterns. The settings can refer to the dynamic range for the attractors of these firing patterns along with the memristor variable coordinate. This hinders the hardware experimental measurement of initial condition-offset boosting firing patterns. In addition, it is not easy to accurately set the initial conditions in each step to acquire corresponding firing patterns from hardware experiments. However, the memristor initial condition-offset boosting firing patterns have potentiality in neuron-based engineering applications [37,48], i.e., the waveform bias of chaotic signal and random signal generation [49]. These deserve our future concern.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.