Multiple Firing Patterns in Coupled Hindmarsh-Rose Neurons with a Nonsmooth Memristor

A model is introduced by coupling two three-dimensional Hindmarsh-Rose models with the help of a nonsmooth memristor. The firing patterns dependent on the external forcing current are explored, which undergo a process from adding-period to chaos. The stability of equilibrium points of the considered model is investigated via qualitative analysis, from which it can be gained that the model has diversity in the number and stability of equilibrium points for different coupling coefficients. The coexistence of multiple firing patterns relative to initial values is revealed, which means that the referred model can appear various firing patterns with the change of the initial value. Multiple firing patterns of the addressed neuron model induced by different scales are uncovered, which suggests that the discussed model has a multiscale effect for the nonzero initial value.


Introduction
As the building elements of the nervous system, the neuron is the most fundamental unit in neural processing. Research on the nonlinear dynamics of neurons is crucial for revealing the mechanisms underlying perception and transmission of neural information.
In 1952, to describe the ionic conductance dynamics of the giant axon, Hodgkin and Huxley established the Hodgkin-Huxley (HH) neuron model [1], which started the research on the neuron model. From then on, other simplified neuron models were proposed successively, which mainly explained lots of ion channels, various synapses, and spatial geometry of individual cells, such as the FitzHugh-Nagumo (FHN) model [2] depicting a prototype of a neuron, the Hindmarsh-Rose (HR) model [3] simulating the characteristics of neurons in the hippocampus of the brain, the Morris-Lecar (ML) model [4] obtained in the research on muscle fiber of Arctic goose, the Chay model [5] as a new theoretical model with unity based on many different types of excitable cells, and the Izhikevich neuron model [6] regarded as a mathematical simplification of the HH model using a binary tree. These models represented a variety of neurons. These neuron models have demonstrated different electrical activities and attracted many researchers' attention. For example, the FHN neuron model exhibits discontinuous transition between different oscillations [7] and double coherence resonance induced by phase noise [8]; the HH neuron model displays evoking spiking caused by enough noise intensity [9], chaotic resonance dependent on current intensity [10], and extrinsic stochastic resonance caused by ion shot noise [11]; in the presence of periodic input, the HR neuron model can show nonlinear resonance behavior [12], periodic and chaotic firing patterns [13], transition between chaotic firing and periodic firing [14], and bursting phenomenon [15]; the Izhikevich neuron model can appear chaotic resonance [16,17]; the ML neuron model can exhibit mono-and bistable dynamic regimes [18] and responses to two temperature-sensitive ion channels, calcium and leak current, respectively [19]. These classical models and their dynamical analysis are motivating researchers to develop more realistic or refined neuron models.
In recent years, improving the classical neuron models by different schemes, such as by introducing the electrical magnetic effect, has received intense attention. For example, considering magnetic flux and Gaussian white noise, a modified Izhikevich model was set up and diverse firing patterns for a different set of parameters were brought to light [20]. An improved HH neuron model was raised, and the firing frequency dependent on the external forcing current was analyzed [21]. A memristor-based HH neuron model was constructed, and multiple patterns of electrical activities as well as stochastic resonance were detected [22]. A stochastic HH model was brought up, and the effects of ionic channel blockage on the electrical activities of it were studied [23]. Kinds of improved FHN models [24][25][26] came up, and the complicated dynamics were revealed, such as noiseenhanced stability and resonance [24], hidden extreme multistability [25], and subcritical Hopf bifurcation and stochastic resonance [26]. The ML neuron model was improved [27][28][29] along with the dynamical properties being explored, like the memristor synapse-based ML model with abundant periodic and chaotic bursting [27], the random dynamical behavior driven by Gaussian white noise [28], and the stochastic hybrid ML model with the extended parameter regime of oscillations related to noise [29]. In particular, many researchers paid much attention to the well-known HR model and revised it utilizing various methods [30][31][32][33][34][35][36][37][38], including impulsive effect and period-adding bifurcation in a hybrid HR model [30]; alternating current-induced coexisting behaviors of asymmetric busters in an external alternating current injected HR model [31]; diversity of firing patterns in a memristor-coupled HR model [32]; multiple firing patterns dependent on the complex electrophysiological condition found in a memristor HR model [33]; coherence resonance induced by phase noise in a 3D memristor-based HR model [34]; the relationship between the energy and the firing pattern in a type of memristor-based HR model [35]; the coexisting phenomenon of diverse firing patterns in a 5D memristor-coupled HR neuron model [36]; hyperchaotic behavior in a kind of neuron model with discontinuous magnetic induction [37]; and different types of firing patterns in a modified Hindmarsh-Rose model, such as square-wave bursting, chattering, fast spiking, periodic spiking, and mixed-mode oscillations [38]. Additionally, a novel neuron model, the Wang-Zhang model, built directly at the neuron level was brought forth [39]. Research results [40,41] suggested that the Wang-Zhang model is equivalent to the HH neuron model from the perspective of neural energy calculation and energy coding. The aforementioned results indicate that to better show electrical activity of the neuron system, many improved neurons were brought up from different perspectives and various dynamical phenomena were gained. It is worth mentioning that among the referred methods to improve the neuron model, memristor coupling was an important approach when considering the effect of the electric field. However, the memristor coupling was mainly about a smooth memristor. Nonsmooth memristor-coupled neurons as well as the dynamics were hardly reported.
Actually, there is often a nonsmooth memristor, which can be used and make the system appear complicated dynamics. Inspired by this idea, in this paper, a nonsmooth memristor is in consideration in the HR neuron model and firing patterns of it are to be discovered. Other parts of this paper are arranged as follows. In Section 2, a novel neuron model is described by coupling two HR neurons using a nonsmooth memristor and different firing patterns are given when changing the external forcing current. In Section 3, the equilibrium points along with their stabilities of the introduced model are analyzed quantitatively. In Section 4, the coexistence of different firing patterns dependent on the initial values is exhibited via numerical simulations. In Section 5, the multiscale feature of the neuron model is discussed. Some conclusions are drawn in Section 6.

Nonsmooth Memristor-Coupled
Neuron Model 2.1. Model Descriptions. In this section, a nonsmooth memductance function [42] is considered to be which is an ideal and active flux-controlled memristor with absolute value nonlinearity, where α and β are two memristor parameters with positive values. When there exists a membrane potential difference between two neurons, electromagnetic induction current will be sensed. For this reason, memristor synapse can be applied to characterizing the dynamics caused by the membrane potential difference. Consequently, on the basis of the 3D HR model [43], a model involving two neurons coupled with a nonsmooth memristor can be presented as where x and u represent the membrane potentials in coupled neurons, y and v mean the exchanges of ions in coupled neuron membranes, and z and w denote the adaptation currents. I ext is the external forcing current. k is the coupling strength.
x − u expresses the difference of membrane potential. System parameters a and b are the fitting coefficients of a cubic function which is used to describe the rate of change of membrane potential. The values of c and d are required to ensure the time course of adaptation current in voltage-clamp condition. r and S mean the corresponding parameters relative to a short depolarizing current and are utilized to depict the change of adaptation current.

Firing Patterns Affected by the External Forcing Current.
From Section 2.1, it can be known that by adjusting the rate of change of membrane potential, voltage-clamp condition, and change of adaptation current to an appropriate state, the system parameters can be chosen as a = 1:0, b = 3:0, 2 Neural Plasticity c = 1:0, d = 5:0, r = 0:006, and S = 4:0. With these parameter values, the 3D HR model [43] can demonstrate complicated dynamics. To explore the firing patterns of model (2), the coupling strength is taken as k = 0:1 and the initial value is considered x 0 = ð0, 0, 0, 0, 0, 0, 0Þ, and then the sampled time series of neuron model (2) are calculated and drawn in Figures 1-5, from which it can be seen that the addressed model (2) appears multiple firing patterns with the change of the external forcing current I ext while other parameters are kept unchanged. Specifically, for I ext = 1:4, 1:8, 2:5, 2:8, neuron model (2) displays firing patterns of period-1, period-2, period-3, and period-4, respectively. Namely, it shows an adding-period phenomenon. But for I ext = 3:3, neuron model (2) comes into a chaotic state. To better characterize this process of change, Figure 6 draws the bifurcation of membrane potential x with bifurcation parameter I ext , which confirms the results of Figures 1-5. As has been noted, the firing patterns of model (2) are similar to the characteristics of the HR neuron model coupled with a smooth memristor [44].

Equilibrium Points and Their Stability Analysis
In this section, system parameters of neuron model (2)  (2) is provided with a chaotic feature. The equilibrium points of (2) can be determined as where δ 1 and δ 2 can be regarded as the intersection points of function curves It can be solved via the graphic analytic method, which is a method to solve the problem in geometry. According to the conditions and conclusions of the problem to be solved, one or more basic graphics of it are found by analyzing. Then, the properties of these basic graphics are applied to solve the problem. Because ðδ 1 , δ 2 Þ is believed as the intersection point of function curves (4) and (5), it can be obtained by picturing the curves of (4) and (5). Therefore, the intersection coordinate ðδ 1 , δ 2 Þ can be achieved. Substitute the values of δ 1 and δ 2 into (3), and   (4) and (5) can be calculated numerically and pictured in Figures 7 and 8, respectively, which indicate that for k = 0:1 and k = 0:5, the number of intersection points of (4) and (5) is 5 and 3, respectively. It means that the value of coupling coefficients has an effect on the number of intersection points. That is to say, the number of nonzero equilibrium points of model (2) is relative to the value of the coupling coefficient k.
To judge the stabilities of the equilibrium points, the Jacobian matrix of model (2) at A is yielded as where For the given coupling strengths k = 0:1 and 0.5, the values ðδ 1 , δ 2 Þ of the equilibrium points are calculated and listed in Tables 1 and 2. Meanwhile, utilizing the Jacobian matrix (6), the corresponding eigenvalues at A are computed and also given in Tables 1 and 2. By doing so, the stabilities of the equilibrium points can be asserted (see Tables 1 and 2), from which we can conclude that for k = 0:1, the equilibrium points A 1 , A 2 , A 3 , A 4 , and A 5 are all unstable saddle points and for k = 0:5, A 1 and A 3 are stable center points while A 2 is an unstable saddle point. It is known that equilibrium points have complicated distributions, which can result in the multistability of a chaotic system [45].

Coexistence of Multiple Firing Patterns
As it is well known, a chaotic system is sensitive to the initial value. In terms of the neuron system, the change of the initial value can make neurons demonstrate different electrical    Figure 7: Function curves of (4) and (5) as well as the intersection points when k = 0:1. The red curve corresponds to function (4), and the blue one is for function (5). A 1 , A 2 , A 3 , A 4 , and A 5 are five intersection points. 6 Neural Plasticity appears periodic firing and chaotic firing for I ext = 1:4 and I ext = 3:3, respectively. Namely, in model (2), various external forcing currents can lead to multiple firing patterns. In this part, the existence of multiple firing patterns of model (2) dependent on initial values is demonstrated. For this purpose, considering I ext = 1:4, another initial value is chosen as ð−2, 0, 0, 2, 0, 0, 0Þ and the sampled time series for membrane potential in model (2) are computed and pictured in Figure 9. Comparing Figure 1 with Figure 9, it is obvious to see that model (2) shows the period-1 firing pattern for initial value ð0, 0, 0, 0, 0, 0, 0Þ and becomes close to being stable with small tiny vibration around equilibrium points for initial value ð−2, 0, 0, 2, 0, 0, 0Þ. The phase portrait in Figure 10 verifies the coexistence of two kinds of firing patterns. Take I ext = 3:3 and initial value ð−2, 0, 0, 2, 0, 0, 0Þ; the sampled time series for membrane potential in model (2) are computed and drawn in Figure 11. By analyzing Figures 5 and 11, it can be found that model (2) presents a chaotic phenomenon for initial value ð0 , 0, 0, 0, 0, 0, 0Þ and tends to be stable with small tiny vibration around equilibrium points for initial value ð−2, 0, 0, 2, 0, 0, 0Þ. Figure 12 confirms the result of Figures 5 and 11. Equally important, when k = 0:5, I ext = 1:4, and initial values are selected as ð0, 0, 0, 0, 0, 0, 0Þ and ð−2, 0, 0, 2, 0, 0, 0Þ, corresponding time series of membrane potential in model (2) are computed and pictured in Figure 13, which means that model (2) can illustrate different firing patterns for various initial values. That is to say, the coexistence of multiple firing patterns can also be disclosed, which can be tested in Figure 14. When k = 0:5, I ext = 3:3, initial values are also selected as ð0, 0, 0, 0, 0, 0, 0Þ and ð−2, 0, 0, 2, 0, 0, 0Þ, sampled time series for membrane potential in model (2) are counted and drawn in Figure 15, which shows the chaotic characteristic for initial value ð0, 0, 0, 0, 0, 0, 0Þ, while appears to be stable with small tiny vibration around equilibrium points for initial value ð−2, 0, 0, 2, 0, 0, 0Þ. It suggests that model (2) also can exhibit different firing patterns dependent on the initial values. Phase portraits in Figure 16 verify this result.
The above results indicate that when coupling coefficients are taken as k = 0:1 and k = 0:5, whether the firing pattern is periodic or chaotic, it can be changed into an approximately stable state by selecting the initial value. Namely, the coexistence of multiple firing patterns dependent on the initial values can be uncovered, while in the coupled HR model with a smooth memristor, it is hard to detect this coexistence of firing patterns.

Multiscale Effect of the Coupled HR Neuron Model by a Nonsmooth Memristor
Traditional approaches to establish a model often focus on one scale, but multiscale phenomena are part of our daily life whether we explicitly recognize it or not. For instance, as a result of the multiscale dynamics of the solar system, our time can be organized in days, months, and years. Our society is in a hierarchical structure from towns to states, countries, and continents. That is to say, it is not an easy task to think of a situation which is not involved in any multiscale characteristics. Therefore, the multiscale model is crucial in precise modeling and can provide support for the dynamic analysis of systems. For example, considering the effects of different drugs on cardiac electrophysiological activity, the drug-induced multiscale model was addressed and the mechanism of drug-induced changes in cardiac behavior was studied [46]. A multiscale model linking the cell level and the subcellular level was proposed [47], which illustrated the prediction of cancer cell migration. In terms of the inertia and response time of the hydraulic servo system, multitime scales modeling of the hydroturbine governing system was put up and the effects of time scales on the dynamical behavior of the system were analyzed [48]. Existing results suggest that the dynamics of the system are often affected by different scales and have a multiscale effect. Encouraged by the above results, the multiscale effect of the coupled HR neuron model with a nonsmooth memristor is to be investigated. As a matter of fact, in a neuron system, electromagnetic induction current is often caused by the difference of membrane potential between two neurons. Therefore, there are different time scales in model (2). Suppose the original time scale is t, the fast time scale is T 1 , and the slow time scale is T 2 . Variables x, y, z, u, v, and w are related to T 1 , and φ is associated with T 2 . Let T 1 = t and T 2 = εt, and then model (2) can be rescaled as  (4) and (5) as well as the intersection points when k = 0:5. The red curve corresponds to function (4), and the blue one is for function (5). A 1 , A 2 , and A 3 are three intersection points.      System parameters are also chosen as a = 1:0, b = 3:0, c = 1:0, d = 5:0, r = 0:006, S = 4:0, α = 4:0, and β = 5:0, and the external forcing current is taken as I ext = 3:0. When the initial value is considered x 0 = ð0, 0, 0, 0, 0, 0, 0Þ, ε has little effect on the firing pattern of model (8). Therefore, in the following discussions, the nonzero initial value is considered, e.g., x 0 = ð1, 0, 0, 0, 0, 0, 0Þ. Two cases are studied with coupling coefficients k = 0:1 and k = 0:5. The effect of smallscale ε on the firing pattern of model (8) is to be studied.

Neural Plasticity
For this purpose, ε is selected as different values. Sampled time series for membrane potential in model (8) are attained and given in Figures 17 and 18. From Figure 17, it can be known that for k = 0:1, when ε = 0:6, the firing patterns of membrane potentials in model (8) gradually converge to a state vibrating around a point with very low amplitude, which can be called a quasisteady state; when ε = 0:06, the firing patterns of membrane potentials in model (8) tend to be stable; when ε reduces to 0.006, the firing patterns of membrane potentials in model (8) appear a new chaotic phenomenon, in which there appears low-amplitude oscillation in the resting state between two adjacent bursting firings. It means     Neural Plasticity that low-amplitude oscillation and high-amplitude oscillation can take on alternatively. From Figure 18, it is obvious to see that for k = 0:5, with the change of ε, model (8) also presents various firing patterns, which is similar to that for k = 0:1. Figures 17 and 18 Figure 16: Coexistence of multiple firing patterns dependent on initial values for k = 0:5 and I ext = 3:3. The red phase portrait is for initial value ð0, 0, 0, 0, 0, 0, 0Þ while the blue one is for ð−2, 0, 0, 2, 0, 0, 0Þ.

11
Neural Plasticity nonzero initial value. Namely, model (8) has an obvious multiscale effect for the nonzero initial value. It is worth noting that compared with model (8), the firing pattern of the coupled HR model with a smooth memristor has little multiscale effect whether the initial value is zero or nonzero. This result may contribute to the optimization analysis and control of the neuron system.

Conclusions
In this paper, considering that the membrane potential difference between different neurons may be responsible for electromagnetic induction current, a neuron model coupled with a nonsmooth memristor is introduced and the multiple firing patterns are explored.
Firstly, multiple firing patterns induced by different external forcing currents are discussed. It is found that with the change of the external forcing current, the referred neuron model experiences a process of period-1, period-2, period-3, period-4, and chaos, namely, addingperiod bifurcation to chaos. Secondly, the diversity of the number and stability of equilibrium points caused by the coupling coefficient is studied and it is disclosed that different coupling coefficients can make the mentioned model be provided with different numbers of equilibrium points as well as various stable states. Thirdly, the 12 Neural Plasticity coexistence of multiple firing patterns depending on the initial value is revealed, which suggests that the dynamics of the addressed neuron model is sensitive to the initial value. Fourthly, the multiscale effect of the mentioned neuron model is uncovered and results suggest that a small scale has much effect on the firing pattern of the referred neuron model.
As has been noted, the electrical activities in coupled Hindmarsh-Rose neurons with a nonsmooth memristor display diversity. The firing pattern of it is relevant to many factors. By regulating certain factors, some required electrical activities can be achieved, which helps uncover the dynamics of neurons and then provides a basis for controlling the firing rhythm of the nervous system.

Neural Plasticity
Additionally, the results in this paper can be generalized to the cases of three neurons or even more neurons, because in the real nervous system, more than two neurons can be regarded as the coupling between another neuron and the coupled neurons.

Data Availability
The data supporting the results can be found in our paper.

Conflicts of Interest
The authors declare that they have no conflicts of interest in this work.