Superextreme spiking oscillations and multistability in a memristor-based Hindmarsh–Rose neuron model

In this paper, we investigate the occurrence of superextreme spiking (SES) oscillations and multistability behavior in a memristor-based Hindmarsh–Rose neuron model. The presence of SES oscillations has been identified as arising due to the occurrence of an interior crisis. As the membrane current I(t), considered as the control parameter is varied, the system transits from bounded chaotic spiking (BCS) oscillations to SES oscillations. These transitions are captured numerically using geometrical representations like time series plots, phase portraits and inter-spikes interval return maps. The characterization of SES from the BCS oscillations is made using statistical tools such as phase shift analysis and probability density distribution function. The multistability nature has been observed using bifurcation analysis and confirmed by the Lyapunov exponents for two different sets of initial conditions. The numerical simulations are substantiated through real-time hardware experiments realized through a nonlinear circuit constructed using an analog model of the memristor.


Introduction
In the past few decades, a large number of researchers have started investigating neuronal models because of the fact that neuronal dynamics has impacted greatly the developments in artificial intelligence and biological systems. In particular, they have turned their attention towards spiking and bursting oscillations (BOs) in neuronal systems. Many mathematical models have been developed by different researchers to study these spiking and BOs. Examples of these are the Hodgkin-Huxley (HH) [1], Morris-Lecar (ML) [2], FitzHugh-Nagumo (FHN) [3], Hindmarsh-Rose (HR) [4,5] and Izhikevich (IZ) [6]. Though these models have been successful to a certain extent in explaining many aspects of neuronal behavior, they have many limitations as well. For instance, the HH neuron model, which explains with considerable depth, the ionic processes and electrical currents on the membrane surfaces, it is not ideal for specific applications owing to difficul-ties in calculations [1]. The FHN model is a simple model, in fact it is a simplified version of the HH neuron model. In spite of its simplicity, it fails to explain bursting patterns in the neurons [3]. On the other hand, the HR neuron model, which is based on the FHN model, is fairly successful in explaining a wide range of dynamical events and neuronal activities [5]. For example, using the HR neuron model, noise-induced torus bursting oscillations [7], noise-induced spikingbursting transition [8] and complex dynamics in discontinuous magnetic induction [9] have been studied. Very recently, a memristor-based HR model showing up rich dynamical behaviors such as coexistence attractors [10], multi-scroll hidden attractors [11] and hidden coexisting asymmetric behaviors [12] has been reported. In addition, the electrical activity of timedelay memristive neurons studied with Gaussian white noise disrupts [13], phase noise-induced coherence resonance [14] and hidden dynamics in a fractional-order memristive HR model [15] have also been reported.
From a different perspective, the spiking oscillations can be considered as sort of extreme events exhibited by neurons. Extreme events (EEs) are found to occur in dynamical systems, when any one of the observables or dynamical state variables takes on an abnormally large variations in magnitudes when compared to the normal states [16][17][18][19][20]. They have been reported in coupled system of bursting neurons [21], in intermittent large amplitude bursting in HR model [22]. Motivated by the aforementioned studies, we report a memristive Hindmarsh-Rose (HR) system, that exhibits superextreme spiking (SES) oscillations and multistability states which can model the spiking activities of the neurons. We have carried out numerical simulations and real-time hardware experiments to support our claim.
The rest of this work is structured as follows. In Sect. 2, the proposed memristor-based HR neuron model and its linear stability analysis are given. In Sect. 3, the emergence of periodic and chaotic spiking oscillations, are reported. While the mechanism of the existence of the extreme spikes is discussed in Sect. 4, the multistability of spikes is briefly discussed in Sect. 5. The real time hardware experimental observations using an analog circuits are detailed in Sect. 6. The experimental results are found to be in complete conformity to those obtained from numerical studies. Finally in Sect. 7, the conclusions of our investigations are summarized.

Memristive Hindmarsh-Rose neuron model
The Hindmarsh-Rose (HR) model of neurons is a twodimensional one, reported by Hindmarsh and Rose in the year 1982 [4]. It simulates the variations in membrane potential, recovery time and spike length activities in neuronal dynamics. In the year 1984, the Hindmarsh and Rose reported an improved threedimensional HR neuron model [5] for investigating the bursting activity and large amplitude spikes in the neurons. In this paper, we present a memristor-based Hindmarsh-Rose model for the neurons. Memristors are two terminal devices that exhibit a pinched hysteresis characteristics in the charge-flux plane when excited by an alternating excitations. Its existence was first postulated by the circuit theorist Leon O. Chua based on axiomatic considerations in 1971 [23] and was first discovered in nano-scale devices in 2008 by scientists at the Hewlett-Packard labs [24]. Ever since this discovery, different analog models of memristors [25,[28][29][30] and many fabrications of it in the nano-scale [28] have been reported.
The memristive HR neuron model discussed here uses a flux-controlled memristor [25]. A flux-controlled memristor can be described by where W ( ) = k denotes the memductance. The memristive HR model [12] can be described by a set of three coupled first-order differential equations, where x, y and are the dynamical state variables representing the membrane potential, spiking amplitude, and magnetic flux, respectively. k x is the external electromagnetic induction, in which k is the electromagnetic strength. The constants a, b, c, d and k are system parameters and I is the membrane input current. While all the constants are kept fixed, the current I is taken as the system control parameter to study the dynamics.
For linear stability analysis, we set up the Jacobian matrix from the Eq. (2) and find the characteristic equation as, where p = −1, From the Eq. (4), the eigenvalues are given as, To study the dynamics of the system defined by Eq. (2), we fix the values for parameters as a = 1.0, b = 3.13, c = 1.0, d = 5.0 and k = 1.0. Pertaining to these parametric values, we obtain three different eigenvalues for Eq. (4), namely (λ 1 , λ 2,3 ) as given in Eq. (5). The positive and complex conjugate eigenvalues (λ 2,3 ) show the presence of an unstable equilibrium point for the system of Eq. (2). This is found using the XPPAUT-AUTO simulation [29]. Here, we have plotted nullcline curves in the (x − y) plane for typical membrane input current I = 1.013 as shown in Fig. 1. In Fig. 1, the red and green color curves represent the x nullcline and y nullcline, respectively. Here, the intersection of x and y nullclines indicates the unstable (US) equilibrium point, (black dot online).

Numerical and statistical studies
The numerical investigations on the memristive HR neuron model of Eq. (2) are examined in this section. The system parameters are set as mentioned in Sect. 2. The system of Eq. (2) is integrated by the well-known RK 4 method for the initial conditions (x, y, ) = (−1.0, −2.0, −3.0) with step size 0.01. By decreasing the control parameter, I in the range of I ∈ (0.1, 2.5), we find the system transits from periodic behavior to chaos to large amplitude spiking behavior and then back to chaos and periodicity. This is seen clearly in the one parameter bifurcation diagram ( Fig. 2a) in the (I − y min ) plane and its corresponding Lyapunov exponents spectra (Fig. 2b) in the (I −L E 1,2 ) plane.

Periodic spiking oscillations
In our investigations, we found that the system exhibits spiking oscillations of different periods as the control parameter, I is decreased. Looking at the bifurcation diagram ( Fig. 2a) in the (I − y min ) plane, we find that the system exhibits period-1 spiking oscillations (P 1 S) for I ∈ (2.5, 2.29), period-2 spiking oscillations (P 2 S) for I ∈ (2.28, 2.0), period-3 spiking oscillations (P 3 S)

Bounded chaotic spiking (BCS) oscillations
When the membrane input current (I ) is decreased in the range of I ∈(1.37, 1.02), the system exhibits bounded chaotic spiking (BCS) oscillations. This is seen in Fig. 2a The time series y(t) for the BCS oscillations is plotted for the control parameter I = 1.2 as shown in Fig. 4a(i). The chaoticity is reflected in the unpredictable amplitudes of the spikes. To differentiate spiking oscillations as being bounded chaotic or extreme or superextreme, the threshold value for the amplitude of the state variable y(t) is determined by the relation where P n is the time-average of negative peaks of y(t), σ is its standard deviation and n is an arbitrary integer [16,30,31]. Generally, for extreme spiking oscillations the integer can be any value in the range 4 < n ≤ 8 while for superextreme spiking oscillations the integer is specified as n > 8. For our system of Eq.
(2), we have assumed n = 4 for extreme spiking and calculated the threshold value as H s = −21.31. This is indicated by the (red color dashed line) in the time plot of y(t) in Fig. 4a(i). As the chaotic spiking amplitudes do not exceed this threshold limit, namely H s = −21.31, we can say that the spiking is a bounded chaotic one. As an additional measure of confirmation, we have plotted the probability distribution function (PDF) in the (P n −PDF) plane for the same variable y(t) for the same value of the control parameter, I = 1.2 in Fig. 4a(ii). From this figure, we find that there are no events exceeding the threshold amplitude (H s ), shown in this case by the red color vertical line. This also confirms that the oscillations are bounded chaotic spiking (BCS) oscillations.  Fig. 4b(i). We find the superextreme spiking oscillations occurring for this control parameter. To identify superextreme spiking oscillations, we set the integer value as n = 9 in Eq. (6) and calculate the threshold as H s = −43.15. This is indicated by pink color horizontal line in the time series plot of y(t) variable in Fig. 4b(i). As the spiking amplitudes exceed this threshold value, we confirm that the oscillations are superextreme spiking oscillations. Also, we identify the extreme spiking oscillations and we set the integer value as n = 4 in Eq. (6) and calculate the threshold as H s = −31.17. This is indicated by red color horizontal line in the time series plot of y(t) variable in Fig. 4b(i). Similar to the case of bounded chaotic spiking oscillations, we provide an additional confirmation using the probability distribution function (PDF) plot in the Fig 4b(ii). Here, the threshold for superextreme spiking oscillations is indicated by the vertical dashed line in pink color and we find a definite probability of spiking events beyond this threshold, thereby lending credit to our assertion that the spiking is superextreme oscillations.

Phase slips analysis
The bounded chaotic spiking (BCS) and superextreme spiking (SES) oscillations can be differentiated using phase slips analysis. For that, we have calculated the instantaneous phase difference ( φ) between the phase values of the variable y(t) and the electromagnetic induction ( ). The time series for the variable y(t) for the control parameter (I =1.2) is plotted in Fig. 5a and its corresponding phase slips diagram is plotted in Fig. 5b. While the time series, y(t) shows the absence of large amplitude spikes, the phase slips diagram shows changes in phase only in small integral multiples of π . On the other hand, similar time series and phase slips diagram plotted for superextreme spiking oscillations for (I =1.013) in Figs. 6a and b, respectively, shows the large amplitude spiking oscillations as well as large values of phase slips of the order of a few hundreds ( φ = n100, n = 1, 2, 3..).

Return maps
Further, we have calculated the return map as an additional tool for differentiating bounded chaotic spiking (BCS) from the superextreme (SES) oscillations. For this, we first calculate the inter-spike interval (ISI), that is, the time that has elapsed between two successive spikes. Then, using this ISI value, we calculate the corresponding discrete times, say (n) and (n + 1).
Using the values of y(t) variable at these instances, namely (y n ) and (y n+1 ) we have plotted the return maps in the (y n − y n+1 ) plane. The phase portraits in the (x − y) plane and the corresponding return map in the (y n − y n+1 ) plane for the BCS attractor are plotted in Fig. 7a(i) and a(ii) for the value of control parameter value I = 1.2. Similarly, the phase portrait of SES attractor in the (x − y) plane and its corresponding return map in the (y n − y n+1 ) plane for I = 1.013 are plotted in Figs. 7b(i) and b(ii), respectively. The random and long range distribution of points representing the ISI of the SES attractor in Fig. 7b(ii) when compared to the short range and apparently regular distribution of ISI points in Fig. 7a(ii), once again confirms the existence of SES oscillations in the memristor-based HR neuron model.
In addition, by decreasing the value of control parameter, I further, we have also observed the following transitions after the SES oscillations in the bifurcation diagram of Fig. 2. In the range of control parameter, I ∈(0.863, 0.788), the system exhibits periodic window (PW) spiking oscillations. It then transit into

Mechanism of extreme spikes
The presence of superextreme spiking in this memristive HR neuron model of Eq. (2), can be attributed to the occurrence of an interior crisis. This can be picturized by plotting the phase portrait in the (x − y) plane superimposed on the nullcline plots as shown in Fig. 8. For the values of the current higher than a critical value, say I =1.370, the system exhibits bounded chaotic spikings. If the control parameter is reduced below this critical value to say I =1.013, an unstable periodic orbits collides with a chaotic spiking orbits, causing a huge blowup in the amplitude of the oscillations. These points of intersection are called as the interior crisis points (ICPs) and are shown in Fig. 8. The unstable fixed point is shown by the intersection of x and y nullclines. When the ICPs lie close to the y− nullcline, the amplitude of variable x(t) will be blown up, while the amplitude of y(t) variable will be blown up, if the ICPs lie close to the x− nullcline. From the Fig. 8, we find that the ICPs lie close to the x− nullcline and consequently the amplitude of y(t) variable is blown up as the phase trajectory winds its way in the phase space.

Multistability
Multistability refers to the coexistence of multiple asymptotic stable states of dynamical system for a given set of parameters. It offers a significant deal of diversity in terms of dynamics, as each attractor reflects a different dynamical state of the system. Here, the initial conditions play a significant role in the choice of the ultimate asymptotic state that a multistable system can take [32]. In particularly, the multistability in memristor-based neuron models were reported in [11,12]. In this section, we explore the multistability nature of our system, namely the memristive HR neuron model, Eq. (2). For studying this, we fix the system parameters as (a = 1.0, b = 3.13, c = 1.0, d = 5.0, k = 1.0) and vary the control parameter I in the range of I ∈ (0.1, 2.5) and construct the bifurcation diagram and its Lyapunov exponents, shown in (Figs. 9a and b), respectively, for two different initial conditions. The bifurcation diagram of Fig. 9a, shows two distinct bifurcation curves, one represented by red color dots and the other by blue color dots pertaining to the two distinct initial conditions (−1, −2, −3) and (1, 2, 3), respectively. We have identified two distinct behaviors, namely the bounded chaotic spiking oscillations (BCS) oscillations for the initial conditions (−1, −2, −3) and period-3 spiking oscillations for the initial conditions (1, 2, 3) from the one parameter bifurcation diagram of Fig. 9a.
The corresponding phase portraits in the (x − y) plane is shown in Fig. 10 a(i) and time series of y(t) variable for these two cases for the control parameter, I = 1.31 are shown in Fig. 10 a(i), 10 a(ii), respectively. The BCS attractor represented by blue color

Experimental observations
Finally, we present the experimental investigations of the memristor-based Hindmarsh−−Rose (HR) neuron model of Eq. (2). The experimental circuit shown in Fig. 11 contains operational amplifiers (UA 741), analog multipliers (AD 633), resistors, and capacitors. The DC ±12 V programmable power supply (U8031A) is used for this circuit operation and the output voltages v 1 , v 2 , v 3 are the voltages of capacitors C 1 , C 2 and C 3 , respectively. For further details of the circuit (Fig. 11), see in Ref. [12]. The data acquisition of the output voltages was made using the mixed signal oscilloscope (MSO−X3014A) and the resolution is 100 MHz, 4 GSa/s (Fig. 12). In our experimental studies, we had fixed the circuit parameters as C 1 = C 2 = C 3 = 10 nF, R 1 , R 2 , R 3 , R 5 , R 6 , R 8 , R 9 , R 10 =10 K , R 4 = 243 , R 7 = 1.5 K , R 11 = 106 , R 12 = 55 , R 13 = 100 , R C = 97.8 K , and had kept the variable resistor R I as a control parameter. When this control parameter R I was varied in the range of R I ∈(1 K , 100 K ), we had observed periodic spiking and superextreme spiking oscillations, namely period 1 spikes for R I ∈(5.08 K , 7.77 K ), period 2 spikes for R I ∈(7.78 K , 21.54 K ), period 3 spikes for R I ∈(21.55 K , 23.62 K ) and period 6 spikes for R I ∈(23.63 K , 41.5 K ). These experimental observations are summarized in Fig. 13. On varying the value of control parameter R I further, we observed the bounded chaotic spiking oscillations for R I ∈(41.6 K , 42.2 K ) and superextreme spiking oscillations for R I ∈(42.3 K , 44.38 K ) as are shown in Fig. 14.
To characterize the behavior of the system and differentiate the BCS oscillations from the SES oscillations, we used the same statistical tool employed in our numerical studies, namely the probability distribution function using time series data of the variable v 2 (t) acquired using the mixed signal oscilloscope (MSO). The calculated probability distribution func- t) (x-axis: 500 ms/div, y-axis: 100 mV/div) across the capacitor C 2 and (ii) the corresponding the typical phase portraits in the (v 1 − v 2 ) plane (x-axis: 60 mV/div, y-axis: 100 mV/div) capturing a bounded chaotic spiking (BCS) oscillations for R I = 41.6 K and b superextreme spiking (SES) oscillations for R I = 42.36 K , respectively tion (PDF) is plotted in the (P n −PDF) plane as shown in Fig. 15 for both BCS and SES oscillations. From the Fig. 15a, the red color dashed line represents the threshold value H s = (μ + nσ ) = −3.768 with n = 4. We find that the amplitude of the bounded chaotic spiking (BCS) oscillations are well below this threshold for superextreme spiking (SES) oscillations. Naturally for the amplitudes exceeding these levels, we find the extreme spiking ES oscillations and the superextreme spiking oscillations SES occurring in the system.

Conclusion
In our work, we have investigated numerically the dynamics of the memristor-based Hindmarsh-Rose (HR) neuron model and have identified the presence of superextreme oscillations and multistable states in it.
We have found out that the superextreme spiking (SES) oscillations arise due to an interior crisis. We have observed the SES oscillations geometrically using time series plots, phase portraits and return maps. Also, confirmed the presence of SES oscillations and differentiated them from bounded chaotic spiking (BCS) oscillations using the statistical tools namely, phase slips analysis and probability distribution functions (PDFs). That the system exhibits multistability has been observed by drawing the bifurcation diagrams, Lyapunov exponent spectrum, time series plots and phase portraits for two different sets of initial conditions. Further, we have substantiated our numerical results with experimental investigations using an analog circuit model of the memristor. We have used mixed storage oscilloscope (MSO−X3014A) for data acquisition. We hope that the methodologies used in our studies and the expertise gained in our experimental investigations can be profitably used to develop mathematical models and an analog nonlinear circuits for many neuronal systems and understand the underlying causes for their varied dynamics.