In recent years, the fractional differential equations have attracted more and more attention from researchers at home and abroad. Due to the wide application background, they have played an important role in mathematics [1], physics [2], biology [3], engineering [4] and other fields [5]. As a branch of calculus, fractional calculus is the extension of ordinary integer order differential and integral to arbitrary real order differential and integral. With the increasing maturity of fractional calculus theory and the rapid development of computer science and technology, it is found that using fractional calculus to build a model can accurately describe the properties and morphology of matter, and the fractional calculus model also has a larger stable region compared with the integer model. The theoretical study of fractional differential system not only has theoretical guiding significance, but also has a very wide application prospect. At present, the study of fractional neurodynamics has become an important branch in the field of biological nervous system.

The neuronal system connected by a considerable number of neuron cells is a very complex nonlinear system. Integer order neuronal model can not well describe the neuron movement process, since the fractional order neuronal model is established. We find that the fractional order neuronal model can more accurately describe some characteristics of neuron activities [6]. Hindmarsh-Rose(HR) neuronal model, as a classical neuronal model, was put forward by Hindmarsh and Rose in the 1980s [7]. It plays an important role in describing the physiological processes of neuronal system, neuron excitability and cardiac muscle fiber, and has a good effect in studying the firing behavior of neurons. The rich firing modes also contain abundant bifurcation structures.The study of bifurcation can not only reveal the dynamic mechanism of neuron firing and clarify the intrinsic nature of neural coding, but also reveal the clustering nature of neural network.

Hopf bifurcation is one of the important nonlinear characteristics. It describes the sudden situation of equilibrium stability and the dynamic behavior of nonlinear system such as periodic solution and chaos. By Hopf bifurcation analysis, more detailed information of periodic solutions near the equilibrium point can be obtained. Shi M, et al [8] considered the influence of time delay on fractional HR neurons. A qualitative bifurcation condition is given with time delay as a bifurcation parameter. In [9], based on the stability and bifurcation theory of fractional systems, a two-dimensional improved fractional HR neuron model is analyzed, and a state feedback method is proposed to control the Hopf bifurcation of the model. Huang C, et al [10] studied the dynamics of a class of high-dimensional fractional ring structured neural networks with multiple delays. Taking time delay as bifurcation parameter, the delay-dependent stability and explicit conditions for Hopf bifurcation are given. The results show that the stability and bifurcation of the proposed network depend heavily on the sum of delays, and careful selection of the sum of delays can significantly improve the stability of the network. It is further proved that the order and number of neurons have great influence on the stability and bifurcation of the network.

We arrange the paper as follows: In section 2, some knowledge of fractional differential is briefly introduced. The fractional-order HR neuronal model is given in Section 3. Stability and Hopf bifurcation analysis are investigated in Section 4. In Section 5, two examples are given to demonstrate the validity of the formula and to compare it with the integer-order model. Finally, some useful conclusions are obtained.

There are several definitions of fractional-order calculus, of which the most common are Gr$ \ddot{u} $nwald-Letnikov definition, Riemann-Liouville definition and Caputo definition. In practical application, Caputo definition is generally used to define the time fractional derivative, while Gr$ \ddot{u} $nwald-Letnikov definition or Riemann-Liouville definition is generally used to define the reciprocal of space fractional derivative. What's more, the physical meaning of the initial condition of the Caputo definition is clear, so we will adopt it.

Definition 2.1 [11] The q-derivative of a function of one variable is defined as follows:

$ \begin{equation} _{a}D^{q}_{t}f(t)=\frac{1}{\Gamma{(m-q)}}\int^{t}_{a}\frac{f^{(m)}(x)}{(t-x)^{q+1-m}}dx, (m-1<q\leq m) \end{equation} $

(2.1)

where $ m $ is the least integer that is not less than $ q(m-1<q\leq m) $, $ a $ and $ t $ are the upper and lower limits of the integral, respectively. $ \Gamma{(m-q)} $ is the Gamma function. $ f^{(m)}(x) $ is the m-th derivative of the function $ f(x) $.

Consider the fractional differential linear system as follows:

$ \begin{equation} \begin{cases} D^{q}x=Ax\\ x(0)=x_{0} \end{cases} \end{equation} $

(2.2)

where $ 0<q\leq 1 $, $ x\in R^{n} $, $ A\in R^{n\times n} $, $ D^{q} $ is the fractional differential operator of order $ q $.

The stability of fractional-order linear systems has been thoroughly studied, and the necessary and sufficient conditions for the stability of fractional-order linear systems have been obtained.

Lemma 2.1 [12] System(2.2) is asymptotically stable if and only if $ |\arg(\lambda)|>\frac{q\pi}{2} $; System (2.2) is stable if and only if $ |\arg(\lambda)|\geq\frac{q\pi}{2} $; System (2.2) is unstable if and only if $ |\arg(\lambda)|<\frac{q\pi}{2} $, where $ \lambda $ is the eigenvalue of the matrix $ A $.

Consider the fractional differential nonlinear system as follows:

$ \begin{equation} D^{q}x=g(x) \end{equation} $

(2.3)

where $ 0<q\leq 1 $, $ x\in R^{n} $ and $ g\in R^{n} $ is a nonlinear function, $ D^{q} $ is the fractional differential operator of order $ q $.

Lemma 2.2 [13] If $ |\arg(\lambda)|>\frac{q\pi}{2} $, $ \lambda $ is the eigenvalue of the Jacobian $ J=\frac{\partial g}{\partial x} $ of system (2.3) at the equilibrium $ x_{0} $, that is, $ x_{0} $ satisfies $ g(x_{0})=0 $, then system (2.3) is locally asymptotically stable at the equilibrium point $ x_{0} $.

Hindmarsh-Rose neuronal model is a mathematical expression of neuron cluster behaviors, which can imitate repeated peaks and irregular behaviors of mollusk neurons [14, 15]. The fractional Hindmarsh-Rose neuronal model used here has the following form:

$ \begin{equation} \begin{cases} D^{q}_{t}x(t)=y-ax^{3}+bx^{2}-z+I\\ D^{q}_{t}y(t)=c-dx^{2}-y\\ D^{q}_{t}z(t)=r[s(x-\bar x)-z] \end{cases} \end{equation} $

(3.1)

where $ x $, $ y $, $ z $ represent the membrane potential, the recovery variable and the slow adaption current, respectively. $ a $, $ b $, $ c $, $ d $, $ r $, $ s $, and $ \bar x $ are system parameters, $ I $ is the external stimulus. $ D^{q}_{t} $ is the differential operator by definition of Caputo, $ 0<q\leq 1 $ is the number of fractional order.

We adopt the method of Pyragas [16, 17], adding a time-delay difference feedback $ k[x(t-\tau)-x(t)] $ to the membrane potential of system (3.1), then the system is described by

$ \begin{equation} \begin{cases} D^{q}_{t}x(t)=y-ax^{3}+bx^{2}-z+I+k[x(t-\tau)-x(t)]\\ D^{q}_{t}y(t)=c-dx^{2}-y\\ D^{q}_{t}z(t)=r[s(x-\bar x)-z] \end{cases} \end{equation} $

(3.2)

where $ k $ is the feedback gain and $ \tau $ is the time delay.

As we know, the equilibrium of Eq.(3.1) is the same as Eq.(3.2).

Assume that $ E^{*}(x^{*}, y^{*}, z^{*}) $ is the equilibrium of Eq.(3.1) which satisfies the following equations:

$ \begin{equation} \begin{cases} y^{*}-ax^{*3}+bx^{*2}-z^{*}+I=0\\ c-dx^{*2}-y^{*}=0\\ r[s(x^{*}-\bar x)-z^{*}]=0 \end{cases} \end{equation} $

(4.1)

The Jacobian matrix of Eq.(3.1) near equilibrium $ E^{*} $ is:

$ \begin{equation} J=\left(\begin{matrix}-3ax^{*2}+2bx^{*}&1&-1\\ -2dx^{*}&-1&0\\ rs&0&-r\end{matrix}\right) \end{equation} $

(4.2)

In this paper, we always choose $ a=1 $, $ b=3 $, $ c=1 $, $ d=5 $, $ s=4 $, $ \bar x =-1.6 $, $ r=0.005 $, and $ I=1.7 $, so the single HR neuron without time-delay feedback control exhibits periodic firing with period-1. We can get the equilibrium $ (-1.2147, -6.3772, 1.5413) $ with eigenvalues: $ \lambda_1 =-12.7468 $, $ \lambda_2 =0.0137+0.0348i $, $ \lambda_3 =0.0137-0.0348i $.

From Lemma 2.2, we calculate that $ |\arg(\lambda_{2, 3})|=\frac{q_0 \pi}{2} $, then $ q_0=0.7612 $, it means the critical order of periodic firing in HR neuronal model. Thus, when $ 0<q<q_0=0.7612 $, the equilibrium $ (-1.2147, -6.3772, 1.5413) $ is stable, that means the neuron is in a rest state; when $ q_0=0.7612<q\leq1 $, the equilibrium $ (-1.2147, -6.3772, 1.5413) $ is unstable, and the neuron is in a firing state.

The following cases are given in Fig. 1 with $ q $ being $ 0.6 $, $ 0.7 $, $ 0.8 $, and $ 0.9 $, respectively.

From Fig. 1 (a) and (b), we can see that when $ q=0.6 $ and $ 0.7 $, the equilibrium $ (-1.2147, -6.3772, 1.5413) $ is stable, the neuron goes through a slight movement and gradually returns to a rest state; while $ q=0.8 $ and $ 0.9 $ shown in Fig. 1 (c) and (d), the equilibrium $ (-1.2147, -6.3772, 1.5413) $ is unstable, meanwhile, the neuron exhibits firing with period-1. Comparing these two situations from (c) and (d), we can also find that with the increase of fractional order $ q $, the firing frequency of neuron slows down, and the firing peak increases instead.

In this section, we will investigate dynamics of the Hopf bifurcation of the controlled fractional system (3.2).

By linear transformation, Eq.(3.2) can be linearized to be the follow equations:

$ \begin{equation} \begin{cases} \tilde{x}=x-x^{*}\\ \tilde{y}=y-y^{*}\\ \tilde{z}=z-z^{*} \end{cases} \end{equation} $

(5.1)

$ \begin{equation} \begin{cases} D^{q}_{t}\tilde{x}(t)=\tilde{y}+(-3ax^{*2}+2bx^{*})\tilde{x}-\tilde{z}+k[\tilde{x}(t-\tau)-\tilde{x}(t)]\\ D^{q}_{t}\tilde{y}(t)=-2dx^{*}\tilde{x}-\tilde{y}\\ D^{q}_{t}\tilde{z}(t)=r(s\tilde{x}-\tilde{z}) \end{cases} \end{equation} $

(5.2)

The equilibrium $ (x^{*}, y^{*}, z^{*}) $ of Eq.(3.2)becomes the equilibrium $ (0, 0, 0) $ of Eq.(5.2). Then Eq.(5.2) can be rewritten to be the following vector form:

$ \begin{equation} D^{q}_{t}X(t)=BX(t)+KX(t-\tau) \end{equation} $

(5.3)

where

$ X=(\tilde{x}, \tilde{y}, \tilde{z})^{T}, B=\left(\begin{matrix}-3ax^{*2}+2bx^{*}-k&1&-1\\ -2dx^{*}&-1&0\\ rs&0&-r\end{matrix}\right), K=\left(\begin{matrix}k&0&0\\ 0&0&0\\ 0&0&0\end{matrix}\right). $

The corresponding characteristic equation of Eq.(5.3) is

$ \begin{equation} h(s)=|s^{q}I-B-Ke^{-s\tau}|=s^{3q}+P_{1}s^{2q}+P_{2}s^{q}+P_{3}+(Q_{1}s^{2q}+Q_{2}s^{q}+Q_{3})e^{-s\tau}=0 \end{equation} $

(5.4)

where

$ \begin{equation} \begin{cases} P_{1}=3ax^{*2}-2bx^{*}+k+r+1, \\ P_{2}=(3ax^{*2}-2bx^{*})(r+1)+r(k+s+1)+k+2dx^{*}, \\ P_{3}=r(3ax^{*2}-2bx^{*}+k+s+2dx^{*}), \\ Q_{1}=-k, \\ Q_{2}=-k(r+1), \\ Q_{3}=-kr. \end{cases} \end{equation} $

(5.5)

Next, we consider the effect of time-delay $ \tau $ on the stability of the equilibrium $ (0, 0, 0) $ of Eq.(5.2).Since the roots of the characteristic equation (5.4) will be affected by time delay $ \tau $, a small change in time delay $ \tau $ will inevitably lead to a change in the roots of Eq.(5.4). If there is a critical time-delay, Eq.(5.4) has a pair of pure imaginary roots with it. Then at this critical value, the stability of the equilibrium $ (0, 0, 0) $ of Eq.(5.2) will change, and under certain conditions, a cluster of periodic solutions with small amplitude will bifurcate from the equilibrium $ (0, 0, 0) $. In other words, there is a Hopf bifurcation at the equilibrium $ (0, 0, 0) $.

Assume that the equation has a pure imaginary root $ s=i\omega=\omega(\cos{\frac{\pi}{2}}+i\sin{\frac{\pi}{2}}) $ with $ \omega>0 $, then it satisfies the characteristic equation (5.4), substitute $ i\omega $ into Eq.(5.4) and separate the real part and imaginary part, so we have

$ \begin{equation} \begin{cases} \bar{A}+Q_{1}\omega^{2q}\cos{(q\pi-\omega\tau)}+Q_{2}\omega^{q}\cos{(\frac{q\pi}{2}-\omega\tau)}+Q_{3}\cos{(\omega\tau)}=0\\ \bar{B}+Q_{1}\omega^{2q}\sin{(q\pi-\omega\tau)}+Q_{2}\omega^{q}\sin{(\frac{q\pi}{2}-\omega\tau)}-Q_{3}\sin{(\omega\tau)}=0 \end{cases} \end{equation} $

(5.6)

where

$ \begin{equation} \begin{cases} \bar{A}=\omega^{3q}\cos{(\frac{3q\pi}{2})}+P_{1}\omega^{2q}\cos{(q\pi)}+P_{2}\omega^{q}\cos{(\frac{q\pi}{2})}+P_{3}\\ \bar{B}=\omega^{3q}\sin{(\frac{3q\pi}{2})}+P_{1}\omega^{2q}\sin{(q\pi)}+P_{2}\omega^{q}\sin{(\frac{q\pi}{2})} \end{cases} \end{equation} $

(5.7)

Using the trigonometric function formula, Eq.(5.7) can be rewritten as

$ \begin{equation} \begin{cases} M_{1}\cos{(\omega\tau)}+M_{2}\sin{(\omega\tau)}=-\bar{A}\\ M_{2}\cos{(\omega\tau)}-M_{1}\sin{(\omega\tau)}=-\bar{B} \end{cases} \end{equation} $

(5.8)

where $ M_{1}=Q_{1}\omega^{2q}\cos{(q\pi)}+Q_{2}\omega^{q}\cos{(\frac{q\pi}{2})}+Q_{3}, M_{2}=Q_{1}\omega^{2q}\sin{(q\pi)}+Q_{2}\omega^{q}\sin{(\frac{q\pi}{2})} $. Then, we can get

$ \begin{equation} \begin{cases} \cos{(\omega\tau)}=-\frac{\bar{A}M_{1}+\bar{B}M_{2}}{M^{2}_{1}+M^{2}_{2}}\\ \sin{(\omega\tau)}=-\frac{\bar{A}M_{2}-\bar{B}M_{1}}{M^{2}_{1}+M^{2}_{2}} \end{cases} \end{equation} $

(5.9)

According to $ \cos^{2}{(\omega\tau)}+\sin^{2}{(\omega\tau)}=1 $, we can obtain that $ \bar{A}^2+\bar{B}^2=M_1^2+M_2^2 $. Suppose that

$ \begin{equation} F(\omega)=\bar{A}^2+\bar{B}^2-(M_1^2+M_2^2)=0 \end{equation} $

(5.10)

substitute $ \bar{A} $, $ \bar{B} $, $ M_{1} $ and $ M_{2} $ into Eq.(5.10), we obtain the following expression:

$ \begin{equation} F(\omega)=\omega^{6q}+\phi_{1}\omega^{5q}+\phi_{2}\omega^{4q}+\phi_{3}\omega^{3q}+\phi_{4}\omega^{2q}+\phi_{5}\omega^{q}+\phi_{6}=0, \end{equation} $

(5.11)

where

$ \begin{equation} \left\{ \begin{aligned} \phi_{1}=&2P_1\cos{(\frac{q\pi}{2})}, \\ \phi_{2}=&P_1^2-Q_1^2+2P_2\cos{(q\pi)}, \\ \phi_{3}=&2P_3\cos{(\frac{3q\pi}{2})}+2(P_1P_2-Q_1Q_2)\cos{(\frac{q\pi}{2})}, \\ \phi_{4}=&P_2^2-Q_2^2+2(P_1P_3-Q_1Q_3)\cos{(q\pi)}, \\ \phi_{5}=&2(P_2P_3-Q_2Q_3)\cos{(\frac{q\pi}{2})}, \\ \phi_{6}=&P_3^2-Q_3^2.\\ \end{aligned} \right. \end{equation} $

(5.12)

Lemma 5.1 [8] Assume that Eq.(5.11) has at least one positive real root $ \omega_0 $, then when $ \tau=\tau_j $, Eq.(5.4) exists a pair of purely imaginary roots $ \pm{i\omega_0} $, where

$ \begin{equation} \tau_j=\frac{1}{\omega_0}[\arctan{(\frac{\bar{A}M_2-\bar{B}M_1}{\bar{A}M_1+\bar{B}M_2})+j\pi}], j=0, 1, 2, \cdots. \end{equation} $

(5.13)

Next, we use time delay $ \tau $ as a variable parameter to investigate the Hopf bifurcation behavior of Eq.(3.2). Assume that the roots of Eq.(5.4) are $ s(\tau)=\psi(\tau)+i\omega(\tau) $ which is time-delay depended. We substitute $ s(\tau) $ into Eq.(5.4), and according to the implicit function theorem, we get

$ \begin{equation} \frac{ds}{d\tau}=\frac{s(Q_1s^{2q}+Q_2s^q+Q_3)}{(3qs^{3q-1}+2qP_1s^{2q-1}+qP_2s^{q-1})e^{s\tau}+(2qQ_1s^{2q-1}+qQ_2s^{q-1})-\tau(Q_1s^{2q}+Q_2s^q+Q_3)} \end{equation} $

(5.14)

From Lemma 5.1, we know that when $ \tau=\tau_j $, $ \psi(\tau_j)=0, \omega(\tau_j)=\omega_0 $. Define the bifurcation point $ \tau_0=min{\{\tau_j\}}, j=0, 1, 2, \cdots $, we have

$ \begin{equation} Re\left[\frac{ds}{d\tau}\right]|_{\omega=\omega_0, \tau=\tau_0}=\frac{N_1N_3+N_2N_4}{N_3^2+N_4^2} \end{equation} $

(5.15)

where

$ \begin{equation} \left\{ \begin{aligned} N_1=&-[Q_1\omega_0^{2q+1}\sin{(q\pi)}+Q_2\omega_0^{q+1}\sin{(\frac{q\pi}{2})}]\cos{(\omega_0\tau_0)}\\ &+[Q_1\omega_0^{2q+1}\cos{(q\pi)}+Q_2\omega_0^{q+1}\cos{(\frac{q\pi}{2})}+Q_3\omega_0]\sin{(\omega_0\tau_0)}\\ N_2=&[Q_1\omega_0^{2q+1}\cos{(q\pi)}+Q_2\omega_0^{q+1}\cos{(\frac{q\pi}{2})}+Q_3\omega_0]\cos{(\omega_0\tau_0)}\\ &+[Q_1\omega_0^{2q+1}\sin{(q\pi)}+Q_2\omega_0^{q+1}\sin{(\frac{q\pi}{2})}]\sin{(\omega_0\tau_0)}\\ N_3=&3q\omega_0^{3q-1}\sin{(\frac{3q\pi}{2})}+2qP_1\omega_0^{2q-1}\sin{(q\pi)}+qP_2\omega_0^{q-1}\sin{(\frac{q\pi}{2})}\\ &+[2qQ_1\omega_0^{2q-1}\sin{(q\pi)}+qQ_2\omega_0^{q-1}\sin{(\frac{q\pi}{2})}-Q_1\tau_0\omega_0^{2q}\cos{(q\pi)}\\ &-Q_2\tau_0\omega_0^q\cos{(\frac{q\pi}{2})}-Q_3\tau_0]\cos{(\omega_0\tau_0)}-[2qQ_1\omega_0^{2q-1}\cos{(q\pi)}\\ &+qQ_2\omega_0^{q-1}\cos{(\frac{q\pi}{2})}+Q_1\tau_0\omega_0^{2q}\sin{(q\pi)}+Q_2\tau_0\omega_0^q\sin{(\frac{q\pi}{2})}]\sin{(\omega_0\tau_0)}\\ N_4=&-3q\omega_0^{3q-1}\cos{(\frac{3q\pi}{2})}-2qP_1\omega_0^{2q-1}\cos{(q\pi)}-qP_2\omega_0^{q-1}\cos{(\frac{q\pi}{2})}\\ &-[2qQ_1\omega_0^{2q-1}\cos{(q\pi)}+qQ_2\omega_0^{q-1}\cos{(\frac{q\pi}{2})}+Q_1\tau_0\omega_0^{2q}\sin{(q\pi)}\\ &+Q_2\tau_0\omega_0^q\sin{(\frac{q\pi}{2})}]\cos{(\omega_0\tau_0)}-[2qQ_1\omega_0^{2q-1}\sin{(q\pi)}+qQ_2\omega_0^{q-1}\sin{(\frac{q\pi}{2})}\\ &-Q_1\tau_0\omega_0^{2q}\cos{(q\pi)}-Q_2\tau_0\omega_0^q\cos{(\frac{q\pi}{2})}]\sin{(\omega_0\tau_0)}\\ \end{aligned} \right. \end{equation} $

(5.16)

In order to prove the occurrence of Hopf bifurcation in Eq.(3.2), we give the following lemma.

Lemma 5.2 [8] Assume the roots $ s(\tau)=\psi(\tau)+i\omega(\tau) $ of Eq.(5.4), it satisfies $ \psi(\tau_0)=0, \omega(\tau_0)=\omega_0 $, it means that when $ \tau=\tau_0 $, Eq.(5.4) has a pair of purely imaginary roots. Then the following transversality condition

$ \begin{equation} Re\left[\frac{ds}{d\tau}\right]|_{\omega=\omega_0, \tau=\tau_0}=\frac{N_1N_3+N_2N_4}{N_3^2+N_4^2}\neq0 \end{equation} $

(5.17)

ensures the occurrence of Hopf bifurcation near the bifurcation.

By Lemma 5.1 and Lemma 5.2, we can get the following theorem.

Theorem 5.1 When the feedback gain $ k<0 $, if there is a critical time-delay $ \tau_0 $, Eq.(5.11) has at least one positive real root $ \omega_0 $, and the transversality condition Eq.(5.17) holds, the equilibrium of Eq.(3.2) is asymptotically stable as $ \tau<\tau_0 $ and unstable as $ \tau>\tau_0 $, Hopf bifurcation occurs at the equilibrium of Eq.(3.2).

Proof When Eq.(5.11) has at least one positive real root $ \omega_0 $, it means Eq.(5.4) has at least a pair of pure imaginary roots, and the transversality condition Eq.(5.17) holds, then the Hopf bifurcation will occur at the equilibrium of Eq.(3.2).

In this part, we will test the valid of the above theoretical results.

Example 1 We take $ q=0.83>q_0=0.7612 $, $ k=-5 $. From Eq.(5.11), we can obtain three positive real roots, they are $ \omega_1=0.0226 $, $ \omega_2=0.1169 $ and $ \omega_3=0.6182 $, respectively. Then from Eq.(5.13), the critical delays of the Eq.(3.2) producing Hopf bifurcation at the equilibrium $ (-1.2147, -6.3772, -1.5413) $ are $ \tau_1=1.0128 $, $ \tau_2=3.2398 $, $ \tau_3=2.2122 $. After verification, the transversality condition (5.17) is satisfied. Given $ \tau=1.0 $, $ \tau=1.5 $, $ \tau=1.6 $, $ \tau=2.3 $, the time evolution of the membrane potential is shown in Fig. 2, and we can find that when $ \tau=1.0<\tau_1=1.0128 $, the equilibrium is asymptotically stable. When $ \tau=1.5>\tau_1=1.0128 $, the equilibrium is unstable, and the neuron is chaotic. When $ \tau=1.6<\tau_3=2.2122 $, the equilibrium is asymptotically stable. When $ \tau=2.3>\tau_3=2.2122 $, the neuron is in periodic firing with period-1. According to Theorem 5.1, this confirms the occurrence of Hopf bifurcation near the critical time-delays.

Example 2 We take $ q=0.86>q_0=0.7612 $, $ k=-8 $. From Eq.(5.11), we can obtain three positive real roots, they are $ \omega_1=0.0204 $, $ \omega_2=0.1224 $ and $ \omega_3=5.0037 $, respectively. Then from Eq.(5.13), the critical delays of the Eq.(3.2) producing Hopf bifurcation at the equilibrium $ (-1.2147, -6.3772, -1.5413) $ are $ \tau_1=61.8556 $, $ \tau_2=1.8982 $, $ \tau_3=0.1779 $. After verification, the transversality condition (5.17) is satisfied. Given $ \tau=0.1 $, $ \tau=0.3 $, the time evolution of the membrane potential is shown in Fig. 3. When $ \tau=0.1<\tau_3=0.1779 $, the equilibrium is asymptotically stable. When $ \tau=0.3>\tau_3=0.1779 $, the equilibrium is unstable, and the neuron is in periodic firing with period-1. According to Theorem5.1, this confirms the occurrence of Hopf bifurcation near the critical time-delays.

Remarks When the feedback gain $ k>0 $, the phenomenon at the equilibrium of the system is different from that when the feedback gain $ k<0 $. Here is an example.

We take $ q=0.98>q_0=0.7612 $, $ k=12 $. From Eq.(5.11), we can obtain two positive real roots, they are $ \omega_1=0.0159 $ and $ \omega_2=0.1439 $, respectively. Then from Eq.(5.13), the critical delays of the Eq.(3.2) producing Hopf bifurcation at the equilibrium $ (-1.2147, -6.3772, -1.5413) $ are $ \tau_1=79.8026 $, $ \tau_2=8.6746 $. After verification, the transversality condition (5.17) is satisfied. Given $ \tau=7.0 $, $ \tau=8.75 $, the time evolution of the membrane potential is shown in Fig. 4. When $ \tau=7.0<\tau_2=8.6746 $, the equilibrium is unstable. When $ \tau=8.75>\tau_2=8.6746 $, the equilibrium is asymptotically stable. From Fig. 4 b, we can see that the equilibrium is in a resting state, but the lines in the time evolution diagram are intermittent and not smooth. Obviously, the effect of Fig. 4 (b) is not good. Through the numerical simulation, it can be shown that Hopf bifurcation can occur at the equilibrium of the system when feedback gain $ k $ is positive value or negative value.

In this paper, the stability and Hopf bifurcation of the fractional-order HR neuronal model with time-delay feedback control are studied. In the uncontrolled model, we get the critical order. When the order of the system is smaller than the critical order, the equilibrium is stable; otherwise, the equilibrium is unstable. In the controlled model, we prove the existence of the critical time-delay which generates Hopf bifurcation at the equilibrium of the system, and prove the reliability of the theory through numerical simulations. From Fig. 2 and Fig. 3, the feedback gain is negative, when the time-delay $ \tau $ is less than the critical delay, the equilibrium of the system is in a rest state; when the time-delay $ \tau $ is greater than the critical delay, the equilibrium of the system is in a firing state. From Fig. 4, the feedback gain is positive, when the time-delay $ \tau $ is less than the critical delay, the equilibrium of the system is in a firing state; when the time-delay $ \tau $ is greater than the critical delay, the equilibrium of the system is in a rest state. In summary, it can prove that Hopf bifurcation will occur at the equilibrium of the system.

In the process of numerical simulations, we find an interesting phenomenon, in [17], to realize the occurrence of Hopf bifurcation of the integer-order HR neuronal model with time-delay feedback control, the feedback gain $ k $ is positive, while, in this paper, for the fractional-order HR neuronal model with time-delay feedback control, the feedback gain $ k $ is positive or negative. When the feedback gain $ k $ is negative, the periodic firing occurs while the feedback gain $ k $ is positive, the periodic firing is suppressed. After verification, in the integer-order HR neuronal model with time-delay feedback control, the value of critical time-delay for generating Hopf bifurcation obtained in this paper is consistent with [17].