# Diverse Dynamic Behaviors and Firing Activities of the Modified Fractional-Order Hindmarsh–Rose Neuronal Model Induced by Fractional-Order

Authors: Xin Yang; GuangJun Zhang; XueRen Li; Dong Wang
Journal: Complexity
 (2021)

---

## Abstract

It is important to investigate the firing activities of neurons, and previous experimental works have shown that fractional-order neuronal models depict the firing rate of neurons more verifiably. In this study, a modified fractional-order Hindmarsh–Rose neuronal model is proposed, and the dynamics and firing activities are investigated. Some novel phenomenon can be found. First, by analyzing numerically and theoretically, the Hopf bifurcation is found to occur when the external direct current stimulus is chosen appropriately. The effects of fractional-order on the bifurcation are also studied. Second, when injecting a direct current stimulus, compared with the integer-order model, the system has more varying dynamic behaviors and firing pattern transitions. Under different external current stimulus, periodic firing patterns and chaotic firing patterns occur when fractional-order changes, but the regions of chaotic firing patterns are different. In other words, the transition mode of periodic firing and chaotic firing induced by fractional-order is different under different external current stimulus. The two-dimensional colored diagram of firing patterns is also investigated. Finally, when injecting periodic current stimulus, regular/irregular bursting, multiple spiking, regular\irregular square wave bursting, and mixed firing mods are found by setting the appropriate fractional-order, amplitude, and frequency of the external current stimulus. Some firing patterns cannot be found in integer-order models. When the amplitude is chosen at appropriate values, the region of frequency when the system displays the mixed firing modes decreases with increasing fractional-order.


---

## Body

## 1. Introduction

The neuron was modeled and analyzed by a large amount of experimental data, and the results show that the firing behaviors of neurons are nonlinear processes. In 1952, Hodgkin and Huxley used equivalent circuits and large amounts of data from experiments to model and analyze the data. They constructed the Hodgkin–Huxley (HH) neuron model through theoretical derivation [1]. In 1961, FitzHugh simplified the variables in the HH model and constructed a low-dimensional model, the two-dimensional FitzHugh–Nagumo (FHN) model [2]. Morris and Lecar summarized the new neuron model (Morris–Lecar (ML) model) from the experimental data of the arctic goose muscle fiber, which is a further simplification of the HH model. In 1982, based on voltage clamp experiment data of snail nerve cells, Hindmarsh and Rose proposed the Hindmarsh–Rose (HR) model [3]. The modified HR model is called the polynomial model, which is the modified HR-type model [4]. This is also an important model for neurons, but research on this model’s firing activities is limited. The modified HR model has many dynamic characteristics and depicts bursting and spiking behaviors successfully. It is significant to investigate the dynamic behavior of the modified HR model. Firing activities are important phenomena in the neuronal system. A number of researchers have used different neuronal models to try to explain some of the phenomena observed in experiments. Studying the motor mechanism of the neuronal dynamical system helps us to understand relevant phenomena in the brain and contribute to the development of artificial intelligence.Neuronal models have been studied numerically and theoretically [5, 6]. In [7], the HR neuronal model was used to develop a new type of a central pattern generator. Neuronal models describe the firing patterns in real neurons, for example [8]. For neuronal model firing, the effects of electromagnetic radiation [9, 10] and extracellular alternating-current feld [11] have been investigated. In [9], chaotic bursting and periodic firing were found in the HR model under an electric field, and Hopf bifurcation was verified in this system. In [10], the mode transition of electrical activities induced by electromagnetic radiation was studied. The two-dimensional discrete HR model was investigated in [12]. Diverse bifurcation points and dynamic behaviors were found. The neuronal system’s synchronization behaviors were investigated in [9, 13]. For collective behaviors in neuronal networks, [14] depicted the dynamics of the neurons and neuronal network, such as the mode transition in electrical activity, pattern formation, and selection in networks of neurons. In recent years, some researchers have investigated modified HR models by analyzing bifurcation points theoretically and firing pattern transitions [15–17]. Reference [15] analyzed the dynamic behaviors of a modified HR model with induced electromagnetic radiation and injected current stimulus, and diverse bifurcation points and firing activities were found. The firing pattern transitions with direct current stimulus and without current stimulus were investigated in [16, 17], respectively. At present, there are few studies on the modified HR model, but this model’s dynamic behaviors can be important to understand the experimental results.The above studies studied the integer-order model only. In results of [18], after analyzing the dynamics of the firing rate with a range of stimulus dynamics, showed that the multiple time scale adaptation is consistent with fractional-order differentiation. Fractional computation is related to the intrinsic properties of dynamical systems [19]. Local memory, nonuniform diffusion effects, and many other aspects require fractional computation [19, 20]. Fractional-order differentiation is a fundamental and general computation that can contribute to efficient information processing, stimulus anticipation, and assessment of frequency-independent phase shifts of oscillatory neuronal firing. In this way, fractional-order differentiation has advantages over integer-order differentiation in the depiction of the firing characteristics of some types of neurons. The fractional-order model can give a more complete picture of neurons’ dynamic characteristics than the integer-order model. In the past, scholars have performed much research on fraction-order dynamical systems and applied them in many fields, such as financial systems [21], biomedicine [22, 23], and the spread of infectious diseases [24]. In terms of the fractional-order neuronal model, the single fractional-order HR neuronal model [25] was investigated, and the transitions of chaotic firing and periodical firing, spike firing, and bursting firing were observed. Other fractional-order models such as the fractional-order Lorenz [26, 27], piecewise-linear (PWL) hyperchaotic system [28], and Rossler system [27, 29] were also investigated, and diverse dynamic behaviors were found. The firing activities of neuronal models need to be further investigated, and modified fractional-order HR model may have more diverse firing patterns. However, there is no research studying the modified fractional-order HR model under external direct and periodic current stimulus. Like [21, 23], the Adomian decomposition methods (ADM) algorithm is used in the present study because of its great computational efficiency.From this perspective, this study is organized as follows: first, the modified fractional-order HR model under an external current stimulus is proposed. The bifurcation analysis and influence of fractional-order on the bifurcation are studied numerically and theoretically. Then, the firing patterns, firing pattern transition, and bifurcation diagrams with a bifurcation parameter of fractional-order when injected with a direct current stimulus were investigated. Finally, when a periodic current stimulus is injected, diverse firing activities are studied.


## 2. Model Description and Bifurcation Analysis

### 2.1. Model Description

After Hindmarsh and Rose proposed the Hindmarsh–Rose model [3], a large number of studies have been constructed on the HR neuronal model. The HR neuronal model is described as follows:(1)x˙=y−ax3+bx2−z+I,y˙=c−dx2−y,z˙=rsx−x¯−z,where x is the membrane action potential, y is a recovery variable, z is a slow adaption current, and I is the external stimulus current. a,b,c,d,r,s,x¯ are the parameters in the system, and they are usually set as a=1,b=3,c=1,d=5,r=0.006,s=4,x¯=−1.56. From previous studies [15–17], the integer-order modified HR model under a current stimulus is described as follows:(2)x˙=−s−ax3+x2−y−bz+I,y˙=φx2−y,z˙=εsa1x+b1−kz,where x represents the membrane voltage of the model; y and z measure the gating dynamics of the ion channels and the dynamics of the cytosolic channels, respectively; a,b,a1,b1,k,s,φ are the parameters of the modified HR model, and in this study, they are set as a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I is a variable.There are many definitions of fractional derivatives, and in practice, the most frequently used definitions are three: the Grunwald–Letnikov derivative, Riemann–Liouville derivative, and Caputo derivative. According to [25], under some conditions, the three definitions are equivalent and can be intertranslated. The Caputo derivative makes the Laplace transformation more concise, so it is simpler to solve the fractional-order derivative. Similar to most literature, the Caputo derivative is adopted in this study.Definition 1.
The Caputo derivative of functionfx is defined as(3)D0Ctqfx=1Γn−q∫0tfnτt−τq−n+1dτ,where n−1<q<n, Γ• is the gamma function which is defined as(4)Γz=∫0∞tz−1e−1dt.
Especially, when0<q<1,(5)D0Ctqfx=1Γ1−q∫0tf′τt−τqdτ.
The fractional-order modified HR neuronal model can be proposed as(6)Dtqx=−s−ax3+x2−y−bz+I,Dtqy=φx2−y,Dtqz=εsa1x+b1−kz,where Dtq is the differential operator defined by Caputo, q is the fractional-order operator, and the other variables and parameters are the same as in the integer-order model. In this study, similar to the integer-order model, the parameters are set as follows: a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I and q are the variables.


### 2.2. Bifurcation Analysis

The equilibrium point can be obtained by solving the equation systemx˙=y˙=z˙=0, which is(7)x˙=−s−ax3+x2−y−bz+I=0,y˙=φx2−y=0,z˙=εsa1x+b1−kz=0.In this part, similar to the integer-order model in the above part,I is variable, and the other parameters are set as a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. The equilibrium point Sexe,ye,ze can be described as follows after calculation:(8)Sexe,xe2,1ksa1xe+b1.From equation system (7) and the equilibrium point Sexe,ye,ze, the equation with one unknown quantity is as follows:(9)−saxe3+s+1xe2+bsa1kxe+bb1k−I=0.To solve the equation, the Cardan formula [5] is introduced. Now, the conversions are m1=−sa, m2=s+1, m3=bsa1/k, m4=bb1/k−I, xe=ye−m2/3m1. Equation (9) can be simplified as follows:(10)ye3+pye+q=0,where p=3m1m3−m22/3m12, q=27m12m4−9m1m2m3+3m23/27m13, and the Cardan discriminant is Δ=q/22+p/32. After calculation, the Cardan discriminant Δ>0, so there exists a single real root which is defined by(11)xe=−q2+Δ3+−q2−Δ3−m23m1.Now, we consider the existence of Hopf bifurcation. Suppose that the systemx˙=fμx,x∈ℝn,μ∈ℝ has an equilibrium x0,μ0 at which the following properties are satisfied:(1)
Dxfμ0x0 has a simple pair of purely imaginary eigenvalues and other eigenvalues with negative real parts(2)
d/dμRλμ|μ=μ0≠0Then, the system has a Hopf bifurcation at the equilibriumx0,μ0. The eigenvalues can be analyzed by determining the Jacobian matrix associated with system (7), and the characteristic equation detJe−λL=0 can be obtained as equation (12), where L is an identity matrix of the same size as Je.(12)Je=3saxe2−2sxe−1−b2φxe−φ0εsa10−kε,Pλ=λ3+F1λ2+F2λ+F3,in which(13)F1=φ−G+kε,F2=2φxe−φG+kεφ−G+bεsa1,F3=kε2φxe−φG+bεsa1φ,G=3saxe2−2sxe.Using the same method as in equation (9), equation (9) can be simplified as follows:(14)λ'3+pλ'2+q=0,in which p1=3F2−F12/3, q1=27F3−9F1F2+3F13/27. The three roots of equation (12) are(15)λ1=−q12+Δ3+−q12−Δ3−F13,λ2=−1+j32−q12+Δ3+−1−j32−q12−Δ3−F13,λ3=−1−j32−q12+Δ3+−1+j32−q12−Δ3−F13.From the above analysis, whenI=Icr=−0.253, there exists a simple pair of purely imaginary eigenvalues and a negative real root at equilibrium x0,μ0. Then, the results can be calculated as(16)dλ1dI|I=Icr≠0.The conclusion that there exists a Hopf bifurcation whenI=Icr can be drawn. For the fractional-order system, the stability can be studied using the following lemma:Lemma 1.
Consider the linear autonomous systemDqX=AX, X0=X0, where X∈Rnn∈N and A∈Rn×n. The system is asymptotically stable when and only when argλ>qπ/2 is satisfied for an arbitrary eigenvalue λ. The system is stable when and only when argλ>qπ/2 is satisfied for an arbitrary eigenvalue λ. The stable area of the q-order linear system is shown in Figure 1.According to Lemma1, the critical values of the external current stimulus Icr for the fractional-order model to start firing at different fractional-orders can be calculated, and the results are shown in Figure 2.Figure 1
Stable area of theq-order linear system.Figure 2
TheIcr∼q curve.


### 2.3. Firing Behavior under Different Direct Current Stimulus and Fractional-Orders

In this section, the parameters are also set asa=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I and q are the variables. In [17], the firing activity of integer-order and the influence of parameters on the firing patterns were studied. In this section, the influence of fractional-order and direct current stimulus are investigated. In [17], when I=0.06, the integer-order modified HR model shows a spiking firing pattern. Figure 3 shows the bifurcation diagram varying with the fractional-order when I=0.06. Figure 3 shows that the firing patterns change with varying fractional-order.Figure 3
Bifurcation diagram with a bifurcation parameter of fractional-order whenI=0.06.Figure4 shows the corresponding time series of x and the phase diagram of x,y. When q=0.2, as shown in Figures 4(a) and 4(b), the neuronal model displays periodic-6 bursting. From Figures 4(c) and 4(d), the neuronal model displays periodic-5 bursting when q = 0.25. When q = 0.3 and q = 0.39, the neuronal model displays chaotic bursting (Figures 4(e), 4(f), 4(i), and 4(j)). The neuronal model displays periodic-4 bursting when q = 0.35 (Figures 4(g) and 4(h)), periodic-2 bursting when q = 0.415 (Figures 4(k) and 4(l)), and spiking when q = 0.55 (Figures 4(m) and 4(n)). From the above analysis, under the influence of the fractional-order, more complex dynamic behaviors occur. The fractional-order can induce a periodic firing pattern for the integer-order model transiting to the chaotic firing pattern.Figure 4
Corresponding time series ofx and the phase diagram of x,y when (a)-(b) q=0.2, (c)-(d) q=0.25, (e)-(f) q=0.3, (g)-(h) q=0.35, (i)-(j) q=0.39, (k)-(l) q=0.415, and (m)-(n) q=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)To prove that if fractional-order can induce a chaotic firing pattern for an integer-order model transitioning to a periodic firing pattern, the firing pattern varies with fractional-order whenI=0.015 is set as an example. In [17], the integer-order neuronal model displays chaotic firing. As shown in Figure 5 (the bifurcation diagram when I=0.015), periodic bursting occurs when the fractional-order decreases. There exists a periodic window around q=0.65.Figure 5
Bifurcation diagram with a bifurcation parameter of fractional-order whenI=0.015.WhenI=0.015, the corresponding time series of x and the phase diagram of x,y are shown in Figure 6. When q=0.8, as shown in Figures 6(a) and 6(b), the neuronal model displays chaotic firing. When q = 0.4, q = 0.3, q = 0.2, the neuronal model displays periodic-3, periodic-4, and periodic-5 bursting (Figures 6(c)–6(h)). Meanwhile, fractional-order can induce the chaotic firing pattern to transition to a periodic firing pattern.Figure 6
Corresponding time series ofx and the phase diagram of x,y when (a)-(b) q=0.8, (c)-(d) q=0.4, (e)-(f) q=0.3, and (g)-(h) q=0.2.
(a)(b)(c)(d)(e)(f)(g)(h)From Figure7, the two-dimensional colored diagram of firing patterns, we can easily find the chaotic firing range and periodic firing range. We find that the chaotic range is larger when the direct current stimulus is smaller.Figure 7
Two-dimensional colored diagram of firing patterns.


### 2.4. Firing Behavior under External Periodic Current Stimulus

In the above analysis, various firing activities are investigated when the modified fractional-order HR model is under direct external current. When the direct external current is replaced by a periodic external current stimulus, there are various firing activities. New firing patterns are found when the fractional-order changes. Spiking, bursting, and other firing patterns can occur at an appropriate fractional-order and by injecting a periodic current stimulus. The modified fractional-order HR model under periodic external current stimulus is as follows:(17)Dtqx=−s−ax3+x2−y−bz+Asinηt,Dtqy=φx2−y,Dtqz=εsa1x+b1−kz,where A is the amplitude of the external current stimulus and η is the angular frequency. In this section, similar to the above section, the parameters are set as a=0.5,b=1, φ=1, ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2, and the fractional-order q, the amplitude A, and angular frequency η are considered as variables.We setA=0.01, and regular bursting and irregular bursting can be found. When q=0.3,η=0.01, the system is in regular bursting, as shown in Figure 8(a). Irregular bursting can be found when q=0.6,η=0.01 (Figure 8(b)). After further investigating the firing activities, the critical value of fractional-order exists between regular bursting and irregular bursting. When the frequency η is changed, we can also find irregular bursting and regular bursting. For example, when η=0.1, the system is in regular bursting when q=0.3 (Figure 8(c)) and irregular bursting when q=0.6 (Figure 8(d)).Figure 8
Firing patterns whenA=0.01, η=0.01, (a)q=0.3, (b)q=0.6, η=0.1, (c)q=0.3, and (d)q=0.6.
(a)(b)(c)(d)The amplitude of the external current stimulus is increased asA=0.05. Multiple spiking consisting of different spiking frequencies, spiking, and square wave bursting can be found. Set η=0.01, and the system is in regular bursting for q=0.3 (Figure 9(a)), multiple spiking when q=0.7 (Figure 9(b)), and square wave bursting when q=0.9 (Figure 9(c)). Setting η=0.05, the system displays mixed firing modes that include bursting and spiking when q=0.4 (Figure 9(d)), irregular bursting when q=0.8 (Figure 9(e)), and another irregular bursting that is different from the above bursting activities when q=0.95 (Figure 9(f)). Furthermore, setting η=0.1, the system displays regular spiking when q=0.9 (Figure 9(g)), mixed mode oscillations that consist of multiple numbers of spiking and subthreshold oscillations when q=0.95 (Figure 9(h)), and the system can also display regular bursting and mixed firing modes when q=0.4 (Figure 9(i)).Figure 9
Firing patterns whenA=0.05, η=0.01, (a)q=0.3, (b)q=0.7, and (c)q=0.9; A=0.05, η=0.05, (d)q=0.4, (e)q=0.8, and (f)q=0.95; A=0.05, η=0.1 (g)q=0.9, (h)q=0.95, and (i)q=0.4.
(a)(b)(c)(d)(e)(f)(g)(h)(i)WhenA=0.1, the firing activities are similar to the above firing patterns. However, a novel phenomenon can be found. First, setting q=0.5, the system is in mixed firing modes when η=0.01 (Figure 10(a)), irregular square wave bursting when η=0.05 (Figure 10(b)), and regular square wave bursting when η=0.1 (Figure 10(c)). Then, setting q=0.4, the system is in mixed firing modes when η=0.01 and η=0.05 (Figures 10(d) and 10(e)), and regular square wave bursting when η=0.1 (Figure 10(f)). When q=0.3, we can find that the system displays mixed firing modes when η=0.01, η=0.05, and η=0.1. From the above analysis, we can conclude that the region of frequency η when the system displays the mixed firing modes decreases with increasing fractional-order.Figure 10
Firing patterns whenA=0.1, q=0.5 (a)η=0.01, (b)η=0.05, and (c)η=0.1; A=0.1, q=0.4 , (d)η=0.01, (e)η=0.05, and (f)η=0.1.
(a)(b)(c)(d)(e)(f)


## 2.1. Model Description

After Hindmarsh and Rose proposed the Hindmarsh–Rose model [3], a large number of studies have been constructed on the HR neuronal model. The HR neuronal model is described as follows:(1)x˙=y−ax3+bx2−z+I,y˙=c−dx2−y,z˙=rsx−x¯−z,where x is the membrane action potential, y is a recovery variable, z is a slow adaption current, and I is the external stimulus current. a,b,c,d,r,s,x¯ are the parameters in the system, and they are usually set as a=1,b=3,c=1,d=5,r=0.006,s=4,x¯=−1.56. From previous studies [15–17], the integer-order modified HR model under a current stimulus is described as follows:(2)x˙=−s−ax3+x2−y−bz+I,y˙=φx2−y,z˙=εsa1x+b1−kz,where x represents the membrane voltage of the model; y and z measure the gating dynamics of the ion channels and the dynamics of the cytosolic channels, respectively; a,b,a1,b1,k,s,φ are the parameters of the modified HR model, and in this study, they are set as a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I is a variable.There are many definitions of fractional derivatives, and in practice, the most frequently used definitions are three: the Grunwald–Letnikov derivative, Riemann–Liouville derivative, and Caputo derivative. According to [25], under some conditions, the three definitions are equivalent and can be intertranslated. The Caputo derivative makes the Laplace transformation more concise, so it is simpler to solve the fractional-order derivative. Similar to most literature, the Caputo derivative is adopted in this study.Definition 1.
The Caputo derivative of functionfx is defined as(3)D0Ctqfx=1Γn−q∫0tfnτt−τq−n+1dτ,where n−1<q<n, Γ• is the gamma function which is defined as(4)Γz=∫0∞tz−1e−1dt.
Especially, when0<q<1,(5)D0Ctqfx=1Γ1−q∫0tf′τt−τqdτ.
The fractional-order modified HR neuronal model can be proposed as(6)Dtqx=−s−ax3+x2−y−bz+I,Dtqy=φx2−y,Dtqz=εsa1x+b1−kz,where Dtq is the differential operator defined by Caputo, q is the fractional-order operator, and the other variables and parameters are the same as in the integer-order model. In this study, similar to the integer-order model, the parameters are set as follows: a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I and q are the variables.


## 2.2. Bifurcation Analysis

The equilibrium point can be obtained by solving the equation systemx˙=y˙=z˙=0, which is(7)x˙=−s−ax3+x2−y−bz+I=0,y˙=φx2−y=0,z˙=εsa1x+b1−kz=0.In this part, similar to the integer-order model in the above part,I is variable, and the other parameters are set as a=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. The equilibrium point Sexe,ye,ze can be described as follows after calculation:(8)Sexe,xe2,1ksa1xe+b1.From equation system (7) and the equilibrium point Sexe,ye,ze, the equation with one unknown quantity is as follows:(9)−saxe3+s+1xe2+bsa1kxe+bb1k−I=0.To solve the equation, the Cardan formula [5] is introduced. Now, the conversions are m1=−sa, m2=s+1, m3=bsa1/k, m4=bb1/k−I, xe=ye−m2/3m1. Equation (9) can be simplified as follows:(10)ye3+pye+q=0,where p=3m1m3−m22/3m12, q=27m12m4−9m1m2m3+3m23/27m13, and the Cardan discriminant is Δ=q/22+p/32. After calculation, the Cardan discriminant Δ>0, so there exists a single real root which is defined by(11)xe=−q2+Δ3+−q2−Δ3−m23m1.Now, we consider the existence of Hopf bifurcation. Suppose that the systemx˙=fμx,x∈ℝn,μ∈ℝ has an equilibrium x0,μ0 at which the following properties are satisfied:(1)
Dxfμ0x0 has a simple pair of purely imaginary eigenvalues and other eigenvalues with negative real parts(2)
d/dμRλμ|μ=μ0≠0Then, the system has a Hopf bifurcation at the equilibriumx0,μ0. The eigenvalues can be analyzed by determining the Jacobian matrix associated with system (7), and the characteristic equation detJe−λL=0 can be obtained as equation (12), where L is an identity matrix of the same size as Je.(12)Je=3saxe2−2sxe−1−b2φxe−φ0εsa10−kε,Pλ=λ3+F1λ2+F2λ+F3,in which(13)F1=φ−G+kε,F2=2φxe−φG+kεφ−G+bεsa1,F3=kε2φxe−φG+bεsa1φ,G=3saxe2−2sxe.Using the same method as in equation (9), equation (9) can be simplified as follows:(14)λ'3+pλ'2+q=0,in which p1=3F2−F12/3, q1=27F3−9F1F2+3F13/27. The three roots of equation (12) are(15)λ1=−q12+Δ3+−q12−Δ3−F13,λ2=−1+j32−q12+Δ3+−1−j32−q12−Δ3−F13,λ3=−1−j32−q12+Δ3+−1+j32−q12−Δ3−F13.From the above analysis, whenI=Icr=−0.253, there exists a simple pair of purely imaginary eigenvalues and a negative real root at equilibrium x0,μ0. Then, the results can be calculated as(16)dλ1dI|I=Icr≠0.The conclusion that there exists a Hopf bifurcation whenI=Icr can be drawn. For the fractional-order system, the stability can be studied using the following lemma:Lemma 1.
Consider the linear autonomous systemDqX=AX, X0=X0, where X∈Rnn∈N and A∈Rn×n. The system is asymptotically stable when and only when argλ>qπ/2 is satisfied for an arbitrary eigenvalue λ. The system is stable when and only when argλ>qπ/2 is satisfied for an arbitrary eigenvalue λ. The stable area of the q-order linear system is shown in Figure 1.According to Lemma1, the critical values of the external current stimulus Icr for the fractional-order model to start firing at different fractional-orders can be calculated, and the results are shown in Figure 2.Figure 1
Stable area of theq-order linear system.Figure 2
TheIcr∼q curve.


## 2.3. Firing Behavior under Different Direct Current Stimulus and Fractional-Orders

In this section, the parameters are also set asa=0.5,b=1,φ=1,ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2. I and q are the variables. In [17], the firing activity of integer-order and the influence of parameters on the firing patterns were studied. In this section, the influence of fractional-order and direct current stimulus are investigated. In [17], when I=0.06, the integer-order modified HR model shows a spiking firing pattern. Figure 3 shows the bifurcation diagram varying with the fractional-order when I=0.06. Figure 3 shows that the firing patterns change with varying fractional-order.Figure 3
Bifurcation diagram with a bifurcation parameter of fractional-order whenI=0.06.Figure4 shows the corresponding time series of x and the phase diagram of x,y. When q=0.2, as shown in Figures 4(a) and 4(b), the neuronal model displays periodic-6 bursting. From Figures 4(c) and 4(d), the neuronal model displays periodic-5 bursting when q = 0.25. When q = 0.3 and q = 0.39, the neuronal model displays chaotic bursting (Figures 4(e), 4(f), 4(i), and 4(j)). The neuronal model displays periodic-4 bursting when q = 0.35 (Figures 4(g) and 4(h)), periodic-2 bursting when q = 0.415 (Figures 4(k) and 4(l)), and spiking when q = 0.55 (Figures 4(m) and 4(n)). From the above analysis, under the influence of the fractional-order, more complex dynamic behaviors occur. The fractional-order can induce a periodic firing pattern for the integer-order model transiting to the chaotic firing pattern.Figure 4
Corresponding time series ofx and the phase diagram of x,y when (a)-(b) q=0.2, (c)-(d) q=0.25, (e)-(f) q=0.3, (g)-(h) q=0.35, (i)-(j) q=0.39, (k)-(l) q=0.415, and (m)-(n) q=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)To prove that if fractional-order can induce a chaotic firing pattern for an integer-order model transitioning to a periodic firing pattern, the firing pattern varies with fractional-order whenI=0.015 is set as an example. In [17], the integer-order neuronal model displays chaotic firing. As shown in Figure 5 (the bifurcation diagram when I=0.015), periodic bursting occurs when the fractional-order decreases. There exists a periodic window around q=0.65.Figure 5
Bifurcation diagram with a bifurcation parameter of fractional-order whenI=0.015.WhenI=0.015, the corresponding time series of x and the phase diagram of x,y are shown in Figure 6. When q=0.8, as shown in Figures 6(a) and 6(b), the neuronal model displays chaotic firing. When q = 0.4, q = 0.3, q = 0.2, the neuronal model displays periodic-3, periodic-4, and periodic-5 bursting (Figures 6(c)–6(h)). Meanwhile, fractional-order can induce the chaotic firing pattern to transition to a periodic firing pattern.Figure 6
Corresponding time series ofx and the phase diagram of x,y when (a)-(b) q=0.8, (c)-(d) q=0.4, (e)-(f) q=0.3, and (g)-(h) q=0.2.
(a)(b)(c)(d)(e)(f)(g)(h)From Figure7, the two-dimensional colored diagram of firing patterns, we can easily find the chaotic firing range and periodic firing range. We find that the chaotic range is larger when the direct current stimulus is smaller.Figure 7
Two-dimensional colored diagram of firing patterns.


## 2.4. Firing Behavior under External Periodic Current Stimulus

In the above analysis, various firing activities are investigated when the modified fractional-order HR model is under direct external current. When the direct external current is replaced by a periodic external current stimulus, there are various firing activities. New firing patterns are found when the fractional-order changes. Spiking, bursting, and other firing patterns can occur at an appropriate fractional-order and by injecting a periodic current stimulus. The modified fractional-order HR model under periodic external current stimulus is as follows:(17)Dtqx=−s−ax3+x2−y−bz+Asinηt,Dtqy=φx2−y,Dtqz=εsa1x+b1−kz,where A is the amplitude of the external current stimulus and η is the angular frequency. In this section, similar to the above section, the parameters are set as a=0.5,b=1, φ=1, ε=0.02,s=−1.61,a1=−0.1,b1=−0.045,k=0.2, and the fractional-order q, the amplitude A, and angular frequency η are considered as variables.We setA=0.01, and regular bursting and irregular bursting can be found. When q=0.3,η=0.01, the system is in regular bursting, as shown in Figure 8(a). Irregular bursting can be found when q=0.6,η=0.01 (Figure 8(b)). After further investigating the firing activities, the critical value of fractional-order exists between regular bursting and irregular bursting. When the frequency η is changed, we can also find irregular bursting and regular bursting. For example, when η=0.1, the system is in regular bursting when q=0.3 (Figure 8(c)) and irregular bursting when q=0.6 (Figure 8(d)).Figure 8
Firing patterns whenA=0.01, η=0.01, (a)q=0.3, (b)q=0.6, η=0.1, (c)q=0.3, and (d)q=0.6.
(a)(b)(c)(d)The amplitude of the external current stimulus is increased asA=0.05. Multiple spiking consisting of different spiking frequencies, spiking, and square wave bursting can be found. Set η=0.01, and the system is in regular bursting for q=0.3 (Figure 9(a)), multiple spiking when q=0.7 (Figure 9(b)), and square wave bursting when q=0.9 (Figure 9(c)). Setting η=0.05, the system displays mixed firing modes that include bursting and spiking when q=0.4 (Figure 9(d)), irregular bursting when q=0.8 (Figure 9(e)), and another irregular bursting that is different from the above bursting activities when q=0.95 (Figure 9(f)). Furthermore, setting η=0.1, the system displays regular spiking when q=0.9 (Figure 9(g)), mixed mode oscillations that consist of multiple numbers of spiking and subthreshold oscillations when q=0.95 (Figure 9(h)), and the system can also display regular bursting and mixed firing modes when q=0.4 (Figure 9(i)).Figure 9
Firing patterns whenA=0.05, η=0.01, (a)q=0.3, (b)q=0.7, and (c)q=0.9; A=0.05, η=0.05, (d)q=0.4, (e)q=0.8, and (f)q=0.95; A=0.05, η=0.1 (g)q=0.9, (h)q=0.95, and (i)q=0.4.
(a)(b)(c)(d)(e)(f)(g)(h)(i)WhenA=0.1, the firing activities are similar to the above firing patterns. However, a novel phenomenon can be found. First, setting q=0.5, the system is in mixed firing modes when η=0.01 (Figure 10(a)), irregular square wave bursting when η=0.05 (Figure 10(b)), and regular square wave bursting when η=0.1 (Figure 10(c)). Then, setting q=0.4, the system is in mixed firing modes when η=0.01 and η=0.05 (Figures 10(d) and 10(e)), and regular square wave bursting when η=0.1 (Figure 10(f)). When q=0.3, we can find that the system displays mixed firing modes when η=0.01, η=0.05, and η=0.1. From the above analysis, we can conclude that the region of frequency η when the system displays the mixed firing modes decreases with increasing fractional-order.Figure 10
Firing patterns whenA=0.1, q=0.5 (a)η=0.01, (b)η=0.05, and (c)η=0.1; A=0.1, q=0.4 , (d)η=0.01, (e)η=0.05, and (f)η=0.1.
(a)(b)(c)(d)(e)(f)


## 3. Conclusion

In this study, the dynamic behaviors of the modified HR model were investigated, and diverse firing activities and some novel phenomena were found. First, according to the integer-order model and Cardan formula, the existence of Hopf bifurcation can be proven. The critical values of the direct external current stimulus for each fractional-order are also obtained. Then, when direct current stimulus is injected, the fractional-order model exhibits more complex dynamic behaviors compared with the integer-order model. In some direct current stimulus, the fractional-order can induce periodic firing for the integer-order model transitioning to chaotic firing, vice versa. In an appropriate direct current stimulus, from the bifurcation diagram with a bifurcation parameter of fractional-order, there is some transition of firing patterns such as chaotic bursting, periodic bursting (including different bursting modes), and periodic spiking. In addition, from the two-dimensional colored diagram, we can clearly find the chaotic firing range and periodic range. Finally, when periodic current stimulus is injected, at some appropriate fractional-order and parameters of periodic current stimulus, the fractional-order model can display regular/irregular bursting, regular/irregular spiking, multiple spiking, regular/irregular square wave bursting, and mixed firing modes. The region of frequency when the system displays the mixed firing modes decreases with the increasing fractional-order at some amplitude values of the external current stimulus.



---

*Source: 8490695-2021-11-12.xml*