Hopf Bifurcation Analysis of a Two-Dimensional Simplified Hodgkin–Huxley Model

Abstract

This paper presents a two-dimensional simplified Hodgkin–Huxley model under exposure to electric fields. The Hopf bifurcations of the simplified Hodgkin–Huxley model are investigated through qualitative analysis and numerical simulations. A necessary and sufficient condition for the existence of Hopf bifurcations is derived, and the conditions for supercritical and subcritical Hopf bifurcations are obtained. Finally, bifurcation diagrams are given for two parameters, and numerical examples are presented to illustrate the effectiveness of the theoretical results.

1. Introduction

The outstanding work of Hodgkin and Huxley on nerve conduction contains many significant results [1]. The Hodgkin–Huxley (H-H) model can be used to determine the action potential in terms of the membrane current, membrane potential, membrane conductance, and the activation and inactivation of certain ions. So far, fruitful results were obtained related to the H-H model [2,3,4,5,6,7,8,9]. The H-H mathematical model describes how action potentials are initiated and propagated in neurons. This set of nonlinear differential equations approximates the electrical characteristics of excitable cells such as neurons and cardiac myocytes. The probability of action potential phase changes can also be calculated. However, as a set of four coupled, highly nonlinear ordinary differential equations, it is difficult to apply detailed mathematical analysis to the H-H model, so it is necessary to impose certain simplifications.

Based on the H-H model, many scholars have used simplified versions to characterize the electrical activity of neurons. FitzHugh et al. simplified the H-H model of spike generation in squid giant axons to a two-dimensional model, called the FitzHugh–Nagumo model [10,11]. Although the H-H model is realistic and biophysically sound, only projections of its four-dimensional phase trajectories can be observed. The FitzHugh–Nagumo model permits the entire solution to be viewed at once. This allows a geometrical explanation of important biological phenomena related to neuronal excitability and spike-generating mechanisms. Morris and Lecar developed a two-dimensional system to reproduce a variety of oscillatory behavior in relation to $C{a}^{2+}$ and $K^{+}$ conductance in the muscle fiber of giant barnacles [12]. The Morris–Lecar model describes the complex relationship between membrane potential and the activation of ion channels within the membrane, showing that the potential depends on the activity of the ion channels, and the activity of the ion channels depends on the voltage. Chay considered three main ion channels across the cell membrane, and identified three variables that describe electrical activity in a nerve model [13]. The Chay model shows that the plasma membrane currents are associated with burst and spike oscillations. Hindmarsh and Rose proposed a phenomenological neuron model as a simplified H-H model. This is capable of mimicking almost all behaviors exhibited by real biological neurons, such as spiking, bursting, and irregular behaviors [14]. Using an algebraic relationship, Moehlis constructed a simplified H-H model as a set of two coupled ordinary differential equations. This allowed the so-called canards (rapid transitions from a small-amplitude limit cycle to a large-amplitude relaxation cycle) that exist in the simplified and original H-H model to be identified [15].

Many simplified H-H models can be divided into two categories. One category includes two or three polynomial differential systems. These simplified H-H models are good at reflecting the dynamic properties of a system, but lose the original H-H structure and biological significance of the neuron model [10,11,14]. The second category of simplified H-H models retains the biological significance of the H-H neuron model, but are still very complicated, making it hard to study their dynamic properties [12,13,15]. A bifurcation is a qualitative change in the dynamical behavior of a nonlinear system as its parameters pass through some critical values [16]. The study of bifurcations in neural models is key to understanding the dynamical origin of many single-neuron and circuit phenomena involved in neural information processing and the organization of behavior [17]. There are many results on the bifurcation analysis of the original H-H model [18,19,20,21,22,23].

The interaction between electric fields and biological tissues has long been of interest. This is primarily because the electric fields can alter neuronal excitability in normal physiological conditions, changing the neuronal ensemble activity of epileptiform behavior [24,25,26,27,28,29,30]. Electric field stimulation has become a powerful form of physiotherapy, so it is important to study the H-H model under exposure to external electric fields. Given the importance of the H-H model, a detailed bifurcation analysis has been applied in the case of external electric fields. Wang and Che studied the Hopf bifurcation of a two-parameter H-H model exposed to external electric fields [31,32,33], but did not give the conditions of supercritical and subcritical bifurcation.

The theory of Hopf bifurcation control may have potential applications in the diagnosis and therapy of dynamical diseases in the H-H model [34,35,36]. The analysis of Hopf bifurcation is significant for the H-H model. An interesting situation that the inhibitory autapse with time delay can induce the resting state change to stable spiking pattern is identified near the subcritical Hopf bifurcation of the Hodgkin–Huxley model in [37]. The generation of double coherence resonances due to white noise is discussed in a three-dimensional reduced H-H model with multiple-timescale feature. In addition, the additive white noise can induce double coherence resonances from the resting state near a subcritical Hopf bifurcation in [38]. Using “the Edge of Chaos” calculation technique, a saddle-node subcritical Hopf bifurcation is a solution of Galvani’s allor-none enigma presented in [5]. So, the conditions for supercritical and subcritical Hopf bifurcations are important and meaningful for the H-H model. As we all know, there are few results about the conditions of supercritical and subcritical Hopf bifurcations for the H-H model. In addition, most existing relevant results about the conditions of supercritical and subcritical Hopf bifurcations are given based on the numerical method, not the analytic method. It is really challenging and important to make a theoretical analysis for the conditions of supercritical and subcritical Hopf bifurcations on the H-H model. Based on the above discussion, this paper presents a simplified H-H model that retains the essential original structure and the biological significance of the neuron model, and reproduces good dynamical properties. The simplified H-H model provides a good approximation of the original H-H model, producing broadly consistent dynamical behaviors. We discuss the Hopf bifurcation in the simplified H-H model under exposure to external electric fields. In addition, we derive a necessary and sufficient condition for the existence of Hopf bifurcations, as well as the conditions for supercritical and subcritical Hopf bifurcations. Bifurcation diagrams are obtained for two parameters, and numerical examples are presented to illustrate the effectiveness of the theoretical results. Our results will be of great assistance in the further investigations of Hopf bifurcations that may occur in the original H-H model.

The remainder of this paper considers the Hopf bifurcations created by two parameters in a simplified H-H model. In Section 2, the simplified two-parameter H-H model is introduced. Section 3 derives the necessary and sufficient condition for the existence of Hopf bifurcations in this model. The conditions for supercritical and subcritical Hopf bifurcations are then obtained. The bifurcation diagrams for the two parameters are presented in Section 4, and numerical examples are conducted to illustrate the effectiveness of the theoretical results. Finally, the conclusions to this study are summarized in Section 5.

2. Simplified H-H Model

In this section, the original H-H model exposed to an external electric field is introduced. The H-H model comprises four coupled nonlinear differential equations, which are difficult to examine by detailed mathematical analysis. So, the original H-H model is simplified as a set of two coupled ordinary differential equations.

2.1. Modified H-H Model Exposed to External Electric Field

The original H-H model [1,3] comprises four coupled nonlinear differential equations describing the temporal evolution of the membrane voltage V under an externally injected current I in a spatially localized axonal compartment. At the cellular level, an exposure of a cell to an external electric field may result in that the induced membrane potential $V_E$ is superimposed on the ionic membrane potential V. Then, the overall membrane potential becomes $V+V_E$. The original H-H model exposed to external electric field is described as follows:

(1)$\begin{array}{c}\frac{dV}{dt}=\frac{1}{C_M}\left[I-{\overline{g}}_{Na}{m}^{3}h(V+V_E-V_{N_a})-{\overline{g}}_K{n}^4(V+V_E-V_K)-{\overline{g}}_l(V+V_E-V_l)\right],\hfill \\ \frac{dm}{dt}={\alpha }_m(1-m)-{\beta }_{m}m,\hfill \\ \frac{dh}{dt}={\alpha }_h(1-h)-{\beta }_{h}h,\hfill \\ \frac{dn}{dt}={\alpha }_n(1-n)-{\beta }_{n}n,\hfill \end{array}$

where $V\left(mV\right)$ is the trans-membrane voltage in muscle cells, m, h, $n\in [0,1]$ are the sodium activation, the sodium inactivation, and the potassium activation variables, respectively, and $t\left(ms\right)$ is the time variable. $C_M$ is the membrane capacitance; ${\overline{g}}_{Na}$, ${\overline{g}}_K$ and ${\overline{g}}_l$ are the maximal conductance of the corresponding ionic currents. In addition, $V_{N_a}$, $V_K$, and $V_l$ are the equilibrium potentials of sodium, potassium, and leak current, respectively. ${\alpha }_x,{\beta }_x\left(m{s}^{-1}\right),$$x\in [m,h,n]$ are functions of V and satisfy the following conditions:

(2)$\begin{array}{c}{\alpha }_m=0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right],\hfill \\ {\beta }_m=4.0exp(-V/18.0),\hfill \\ {\alpha }_h=0.07exp(-V/20),\hfill \\ {\beta }_h=1.0/exp((V+30)/010+1.0),\hfill \\ {\alpha }_n=0.01(10.0-V)/\left[exp\left(\right(10.0-V)/10)-1.0\right],\hfill \\ {\beta }_n=0.125exp(-V/80.0).\hfill \end{array}$

The circuit description of the original H-H model (1) is shown in Figure 1. More information on mathematical issues associated with the original H-H model (1) is given in [3].

2.2. Reduction of the H-H Model

In this paper, we reduce the original H-H model (1) to a set of two coupled ordinary differential equations. Note that Equation (2) for the gating variables m, h, and n can be written in the form:

(3)${\tau }_j\left(V\right)\frac{dj}{dt}=j_{\infty }\left(V\right)-j,$

where

$j_{\infty }\left(V\right)\equiv \frac{{\alpha }_j\left(V\right)}{{\alpha }_j\left(V\right)+{\beta }_j\left(V\right)},{\tau }_j\left(V\right)\equiv \frac{1}{{\alpha }_j\left(V\right)+{\beta }_j\left(V\right)},$

with $j\in [m,h,n]$. The time constant ${\tau }_m$ is much smaller than ${\tau }_n$ and ${\tau }_h$ over the entire relevant range of V. Thus, m evolves faster than n and h. This implies that we could replace $m\left(V\right)$ in Equation (1) by $m_{\infty }\left(V\right)$. Next, when the solution of Equation (1) corresponds to a periodic action potential, the following equation approximately holds [15]:

$n\left(t\right)+h\left(t\right)\simeq 0.8.$

According to [15], we choose the gating variable h by taking $h\left(t\right)=0.8-n\left(t\right)$. Considering a linear form for ${\alpha }_n$ and ${\beta }_n$, Equation (1) can be simplified to

(4)$\begin{array}{c}\frac{dV}{dt}=\frac{1}{C_M}\left[I-{\overline{g}}_{Na}{m}_{\infty }{\left(V\right)}^3(0.8-n)(V+V_E-V_{N_a})-{\overline{g}}_K{n}^4(V+V_E-V_K)-{\overline{g}}_l(V+V_E-V_l)\right],\hfill \\ \frac{dn}{dt}=({\alpha }_1-{\alpha }_{2}V)(1-n)-({\beta }_1-{\beta }_{2}V)n,\hfill \end{array}$

where ${\alpha }_1=0.057$, ${\alpha }_2=-0.0037$, ${\beta }_1=0.125$, ${\beta }_2=0.0015$, and

$m_{\infty }\left(V\right)=\frac{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]}{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]+4.0exp(-V/18.0)}.$

Remark 1.

The reason for the replacement of $m\left(V\right)$ by $m_{\infty }\left(V\right)$ is explained in [3]. The gating variable h by taking $h\left(t\right)=0.8-n\left(t\right)$ is chosen in the simplified H-H model (4), and the rationality of this relationship is given in [15]. Note that the different linear forms of $a_n$ and $b_n$ lead to the different dynamical behaviors of model (4). The detailed affects of the linear forms of $a_n$ and $b_n$ are discussed in [39]. Moreover, the linear forms of $a_n={\alpha }_1-{\alpha }_{2}V$ and $b_n={\beta }_1-{\beta }_{2}V$ (${\alpha }_1=0.057$, ${\alpha }_2=-0.0037$, ${\beta }_1=0.125$ and ${\beta }_2=0.0015$) are chosen, and the simplified model (4) approximates the original H-H model (1) reasonably well. Specifically, the dynamical behaviors of the simplified model (4) are consistent with that of the original H-H model (1) [39,40]. Obviously, model (4) is a reasonable and effective simplification for the original H-H model (1).

3. Hopf Bifurcation Analysis

The Hopf bifurcation analysis of the simplified H-H model (4) is discussed in this section. A necessary and sufficient condition for the existence of Hopf bifurcations about the parameters $I,V_E$ is derived in Theorem 1. In addition, the conditions for supercritical and subcritical Hopf bifurcations of the simplified H-H model (4) are obtained in Theorem 2.

Theorem 1.

Model (4) has Hopf bifurcation points if and only if the equilibrium point $(V,n)$ and parameters $I,V_E$ satisfy the quadratic inequality

$V^2-KV+M<0,$

and the quadratic equation

$T=(\overline{\beta }+D{\lambda }_{m_{\infty }})V^2+BV+{\lambda }_{m_{\infty }}H+A=0$

where

• $\alpha ={\alpha }_1+{\beta }_1$, $\beta ={\alpha }_2+{\beta }_2$, ${\lambda }_{m_{\infty }}=3m_{\infty }^{\prime }\left(V\right)/m_{\infty }\left(V\right)C_M$,

• $K=(2\alpha \beta -br-b{\overline{g}}_l)/{\beta }^2$, $M=({\alpha }^2-bI+br{e}_3+b{\overline{g}}_l{e}_2)/{\beta }^2$,

• $e_1=V_E-V_{N_a}$, $e_2=V_E-V_l$, $e_3=V_E-V_K$, $D={\overline{g}}_K{n}^4$,

• $b=\left({\alpha }_2-\beta n\right)/C_M(0.8-n)$, $r={\overline{g}}_K{n}^3(3.2-3n)$, $\overline{\beta }=\beta +{\lambda }_{m_{\infty }}{\overline{g}}_l$, $H=D{e}_1{e}_3+{\overline{g}}_l{e}_1{e}_2-e_{1}I$,

• $A=(D(V_{N_a}-V_K)+{\overline{g}}_l(V_{N_a}-V_l)-I)/C_M-e_1\alpha$, $B=\beta e_1-\alpha +{\lambda }_{m_{\infty }}E$, $E=D(e_1+e_3)+{\overline{g}}_l(e_1+e_2)-I$

Proof. 

Details of the proof can be found in [39,40], therefore the proof is omitted here. □

According to Theorem 1, model (4) has a Hopf bifurcation if and only if $(V,n)$ and parameters $I,V_E$ satisfy the following conditions:

$\left\{\begin{array}{c}T=(\overline{\beta }+D\left({\lambda }_{m_{\infty }}\right))V^2+(B+D(e_1+e_2){\lambda }_{m_{\infty }})V+{\lambda }_{m_{\infty }}H+A=0,\hfill \\ V^2-KV+M<0.\hfill \end{array}\right.$

Thus, we can use a cubic polynomial of V to discuss the bifurcations in the simplified H-H model.

Remark 2.

According to Theorem 1, we use a cubic polynomial of V to discuss the bifurcations in the simplified H-H model (4). Note that the FitzHugh–Nagumo model [10,11] also uses a cubic polynomial of V to reduce the H-H model. However, the proposed model (4) retains the essential original structure and the biological significance of the H-H model, and reproduces good dynamical properties.

According to the method in [41,42], we give conditions for the supercritical and subcritical Hopf bifurcations of model (4) as follows.

Theorem 2.

Model (4) has a Hopf bifurcation when ω($\omega \in (I,V_E)$) passes through ${\omega }_0$, and

• (1) 

  if ${\mu }_2>0(<0)$, the direction of the bifurcation is supercritical (subcritical),

• (2) 

  if ${\gamma }_2>0(<0)$, the bifurcating periodic solutions are orbitally stable (unstable),

• (3) 

  if ${\tau }_2>0(<0)$, the periods of the bifurcating periodic solutions increase (decrease).

Proof. 

Model (4) has the an equilibrium point at $O(V_0,n_0)$, under the bifurcation parameters $I,V_E$. Based on the following linear transformation of $O(V_0,n_0)$,

(5)$\left\{\begin{array}{c}{V}_1=V-V_0,\hfill \\ n_1=n-n_0,\hfill \end{array}\right.$

Model (4) becomes

(6)$\begin{array}{c}\frac{d{V}_1}{dt}=\frac{1}{C_M}[I-{\overline{g}}_{Na}{m}_{\infty }{(V_1+V_0)}^3(0.8-n_1-n_0)(V_1+V_0+V_E-V_{N_a})\hfill \\ -{\overline{g}}_K{(n_1+n_0)}^4(V_1+V_0+V_E-V_K)-{\overline{g}}_l(V_1+V_0+V_E-V_l)],\hfill \\ \\ \frac{d{n}_1}{dt}=({\alpha }_1-{\alpha }_2(V_1+V_0))(1-n_1-n_0)-({\beta }_1-{\beta }_2(V_1+V_0))(n_1+n_0).\hfill \end{array}$

The second equation of model (6) can then be expressed as

(7)$\frac{d{n}_1}{dt}=(-{\alpha }_1+{\alpha }_2{V}_0-{\beta }_1+{\beta }_2{V}_0)n_1+({\alpha }_2{n}_0-{\alpha }_2+{\beta }_2{n}_0)V_1+({\alpha }_2+{\beta }_2)n_1{V}_1+c_0,$

where $c_0={\alpha }_1-{\alpha }_1{n}_0-{\alpha }_2{V}_0+{\alpha }_2{n}_0{V}_0-{\beta }_1{n}_0+{\beta }_2{n}_0{V}_0$.

The function $m_{\infty }\left(V\right)$ in the first equation of model (6) can be described as

$m_{\infty }\left(V\right)=\frac{{\alpha }_m\left(V\right)}{{\alpha }_m\left(V\right)+{\beta }_m\left(V\right)}=\frac{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]}{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]+4.0exp(-V/18.0)}.$

Then, $m_{\infty }\left(V\right)$ can be written in the following form

$m_{\infty }\left(\overline{V}\right)=\frac{0.1\overline{V}/\left[exp(\overline{V}/10)-1.0\right]}{0.1\overline{V}/\left[exp(\overline{V}/10)-1.0\right]+4.0exp(-\overline{V}/18.0+25/18)},$

where $\overline{V}=25.0-V$. When $\underset{\overline{V}\to 0}{lim}{m}_{\infty }\left(\overline{V}\right)=\frac{1}{1+4exp\left(\frac{25}{18}\right)}$, $m_{\infty }\left(V\right)$ has infinite-order continuous partial derivatives. The Taylor expansion of $m_{\infty }{(V_1+V_0)}^3$ at point $V_0$ is

(8)$m_{\infty }{(V_1+V_0)}^3=\sum _{i=1}^{\infty }{a}_i{V}_1^i.$

The expression in the first equation of model (6)

(9)${\overline{g}}_{Na}{m}_{\infty }{(V_1+V_0)}^3(0.8-n_1-n_0)(V_1+V_0+V_E-V_{N_a}),$

can be written as

(10)${\overline{g}}_{Na}(a_0+a_1{V}_1+a_2{V}_1^2+a_3{V}_1^3+\sum _{i=4}^{\infty }{a}_i{V}_1^i)(0.8-n_1-n_0)(V_1+V_0+V_E-V_{N_a}).$

With $k_1={\overline{g}}_{Na}(0.8-n_0)$, $k_2={\overline{g}}_{Na}(0.8-n_0)(V_0+V_E-V_{N_a})$ and $k_3={\overline{g}}_{Na}(V_0+V_E-V_{N_a})$, we obtain

${\overline{g}}_{Na}(a_0+a_1{V}_1+a_2{V}_1^2+a_3{V}_1^3+\sum _{i=4}^{\infty }{a}_i{V}_1^i)(0.8-n_1-n_0)(V_1+V_0+V_E-V_{N_a})$

$=k_1{a}_0{V}_1+k_1{a}_1{V}_1^2+k_1{a}_2{V}_1^3+k_2{a}_0+k_2{a}_1{V}_1+k_1{a}_2{V}_1^2+k_1{a}_3{V}_1^3-{\overline{g}}_{Na}{a}_0{n}_1{V}_1$

$-{\overline{g}}_{Na}{a}_1{n}_1{V}_1^2-{\overline{g}}_{Na}{a}_2{n}_1{V}_1^3-k_3{a}_0{n}_1-k_3{a}_1{n}_1{V}_1-k_3{a}_2{n}_1{V}_1^2+H_1(V_1,n_1),$

where

$H_1(V_1,n_1)=k_1\sum _{i=3}^{\infty }{a}_i{V}_1^{i+1}+k_2\sum _{i=4}^{\infty }{a}_i{V}_1^i-{\overline{g}}_{Na}\sum _{i=2}^{\infty }{a}_i{n}_1{V}_1^{i+1}-k_3\sum _{i=3}^{\infty }{a}_i{n}_1{V}_1^i.$

Let $k_4={\overline{g}}_K(V_0+V_E-V_K)$, we can derive

$-{\overline{g}}_K{(n_1+n_0)}^4(V_1+V_0+V_E-V_K)-{\overline{g}}_l(V_1+V_0+V_E-V_l)$

$=-6{\overline{g}}_K{n}_1^2{n}_0^2{V}_1-4{\overline{g}}_K{n}_1{n}_0^3{V}_1-{\overline{g}}_K{n}_0^4{V}_1-4k_4{n}_0{n}_1^3-6k_4{n}_1^2{n}_0^2-4k_4{n}_1{n}_0^3-k_4{n}_0^4-{\overline{g}}_l{V}_1-{\overline{g}}_l(V_0+V_E-V_l)-H_2(V_1,n_1),$

where

$H_2(V_1,n_1)={\overline{g}}_K{n}_1^4{V}_1+4{\overline{g}}_K{n}_1^3{V}_1+k_4{n}_1^4.$

Hence, the changed first equation of model (6) becomes

$\begin{array}{ll}\frac{d{V}_1}{dt}& =\frac{1}{C_M}(I-k_1{a}_0{V}_1-k_1{a}_1{V}_1^2-k_1{a}_2{V}_1^3-k_2{a}_0-k_2{a}_1{V}_1-k_1{a}_2{V}_1^2-k_1{a}_3{V}_1^3+{\overline{g}}_{Na}{a}_0{n}_1{V}_1+{\overline{g}}_{Na}{a}_1{n}_1{V}_1^2\\ & \hspace{1em}+{\overline{g}}_{Na}{a}_2{n}_1{V}_1^3+k_3{a}_0{n}_1+k_3{a}_1{n}_1{V}_1+k_3{a}_2{n}_1{V}_1^2-6{\overline{g}}_K{n}_1^2{n}_0^2{V}_1-4{\overline{g}}_K{n}_1{n}_0^3{V}_1-{\overline{g}}_K{n}_0^4{V}_1\end{array}$

(11)$-4k_4{n}_0{n}_1^3-6k_4{n}_1^2{n}_0^2-4k_4{n}_1{n}_0^3-k_4{n}_0^4-{\overline{g}}_l{V}_1-{\overline{g}}_l(V_0+V_E-V_l)-H_1(V_1,n_1)-H_2(V_1,n_1)).$

Presuming that

(12)$\begin{array}{c}{d}_1=-\frac{1}{C_M}(k_1{a}_0+k_2{a}_1+{\overline{g}}_l+{\overline{g}}_K{n}_0^4),\hfill \\ d_2=\frac{1}{C_M}(k_3{a}_0-4k_4{n}_0^3),\hfill \\ d_3={\alpha }_2{n}_0-{\alpha }_2+{\beta }_2{n}_0,\hfill \\ d_4=-{\alpha }_1+{\alpha }_2{V}_0-{\beta }_1+{\beta }_2{V}_0,\hfill \end{array}$

Equation (6) can be simplified as

(13)$\dot{X_1}=A{X}_1+B\left(X_1\right),$

where $X_1=\left(\begin{array}{c}{V}_1\\ n_1\end{array}\right)$, $A=\left(\begin{array}{cc}{d}_1& d_2\\ d_3& d_4\end{array}\right)$, $B\left(X_1\right)=\left(\begin{array}{c}Q(V_1,n_1)\\ P(V_1,n_1)\end{array}\right)$.

The Jacobian matrix of model (6) at point $O(0,0)$ is

(14)$J\left(0\right)=A=\left(\begin{array}{cc}{d}_1& d_2\\ d_3& d_4\end{array}\right),$

and its characteristic equation is

(15)$f\left(\lambda \right)={\lambda }^2-(d_1+d_4)\lambda +d_1{d}_4-d_2{d}_3=0.$

When $d_1+d_4=0$, $d_1{d}_4-d_2{d}_3>0$, the solutions of Equation (15) are ${\lambda }_{1,2}=\pm \sqrt{d_1{d}_4-d_2{d}_3}i$.

The characteristic value is $\lambda =\lambda \left(\omega \right)$, $\omega \in (I,V_E)$. From Equation (15), we obtain

(16)$\frac{d\lambda }{d\omega }=\frac{{(d_1+d_4)}^{\prime }\lambda -{(d_1{d}_4-d_2{d}_3)}^{\prime }}{2\lambda -(d_1+d_4)}.$

Therefore,

(17)$\begin{array}{c}{\alpha }^{\prime }\left(0\right)=Re\left({\lambda }^{\prime }\left(\omega \right)\right)=\frac{{(d_1+d_4)}^{\prime }}{2}\hfill \\ {\omega }^{\prime }\left(0\right)=Im\left({\lambda }^{\prime }\left(\omega \right)\right)=\frac{{(d_1{d}_4-d_2{d}_3)}^{\prime }}{2k},\hfill \end{array}$

where $k=\sqrt{d_1{d}_4-d_2{d}_3}$. From Equation (17), we have that ${\alpha }^{\prime }\left(0\right)\ne 0$, ${\omega }^{\prime }\left(0\right)\ne 0$. Therefore, model (6) has a Hopf bifurcation.

The eigenvector $\nu$ of ${\lambda }_1=ki$ is

$\nu =\left(\begin{array}{c}1\\ \frac{ki-d_1}{d_2}\end{array}\right).$

Define

$T=(Re\nu ,-Im\nu )=\left(\begin{array}{cc}1& 0\\ -\frac{d_1}{d_2}& -\frac{k}{d_2}\end{array}\right)$

and

$\left(\begin{array}{c}{V}_1\\ n_1\end{array}\right)$ = T$\left(\begin{array}{c}{V}_2\\ n_2\end{array}\right)$.

Then, we have

$\left(\begin{array}{c}{V}_1\\ n_1\end{array}\right)$ = $\left(\begin{array}{c}{V}_2\\ -\frac{d_1{V}_2+k{n}_2}{d_2}\end{array}\right),$ from which we obtain

(18)$\left\{\begin{array}{c}\dot{V_2}=-k{n}_2+Q(V_2,n_2),\hfill \\ \\ \dot{n_2}=k{V}_2+R(V_2,n_2),\hfill \end{array}\right.$

where

(19)$\left\{\begin{array}{c}Q(V_2,n_2)=l_1{V}_2^2+l_2{V}_2^3+l_3{V}_2{n}_2+l_4{V}_2^2{n}_2+l_5{V}_2{n}_2^2+l_6{n}_2^2+l_7{n}_2^3+c_1+H(V_2,n_2),\hfill \\ \\ R(V_2,n_2)=-\frac{1}{k}(d_{1}Q(V_2,n_2)+d_{2}P(V_2,n_2)),\hfill \end{array}\right.$

$\left\{\begin{array}{c}{l}_1=\frac{1}{C_M}(-k_1{a}_1-k_2{a}_2-\frac{{\overline{g}}_{Na}{a}_0{d}_1}{d_2}-\frac{k_3{a}_1{d}_1}{d_2}+\frac{4{\overline{g}}_K{n}_0^3{d}_1}{d_2}-\frac{6k_4{n}_0^2{d}_1^2}{d_2^2}),\hfill \\ l_2=\frac{1}{C_M}(-k_1{a}_2-k_2{a}_3-\frac{{\overline{g}}_{Na}{a}_1{d}_1}{d_2}-\frac{k_3{a}_2{d}_1}{d_2}-\frac{6{\overline{g}}_K{n}_0^2{d}_1^2}{d_2^2}+\frac{4k_4{n}_0{d}_1^3}{d_2^3}),\hfill \\ l_3=\frac{1}{C_M}(-\frac{{\overline{g}}_{Na}{a}_{0}k}{d_2}-\frac{k_3{a}_{1}k}{d_2}+\frac{4{\overline{g}}_K{n}_0^{3}k}{d_2}+\frac{12k_4{n}_0^2{d}_{1}k}{d_2^2}),\hfill \\ l_4=\frac{1}{C_M}(-\frac{{\overline{g}}_{Na}{a}_{1}k}{d_2}-\frac{k_3{a}_{2}k}{d_2}+\frac{12{\overline{g}}_K{n}_0^2{d}_{1}k}{d_2^2}+\frac{12k_4{n}_0{d}_1^{2}k}{d_2^3}),\hfill \\ l_5=\frac{1}{C_M}(-\frac{6{\overline{g}}_K{n}_0^2{k}^2}{d_2^2}+\frac{12k_4{n}_0{d}_1{k}^2}{d_2^3}),\hfill \\ l_6=\frac{1}{C_M}(-\frac{6k_4{n}_0^2{k}^2}{d_2^2}),\hfill \\ l_7=\frac{1}{C_M}\left(\frac{4k_4{n}_0{d}_1{k}^3}{d_2^3}\right),\hfill \end{array}\right.$

and $c_1=\frac{1}{C_M}(I-k_2{a}_0-k_3{n}_0^4-{\overline{g}}_l(V_0+V_E-V_l))$, $H(V_2,n_2)=\frac{1}{C_M}(-H_1(V_2,n_2)-H_2(V_2,n_2))$, $P(V_2,n_2)=-\frac{\beta }{d_2}(d_1{V}_2^2+k{V}_2{n}_2)+c_0$, $\beta ={\alpha }_2+{\beta }_2$.

From Equation (19), we can derive

$\begin{array}{c}R(V_2,n_2)=-\frac{1}{k}(d_{1}Q(V_2,n_2)+d_{2}P(V_2,n_2))\hfill \\ =-\frac{1}{k}(d_1{l}_1{V}_2^2+d_1{l}_2{V}_2^3+d_1{l}_3{V}_2{n}_2+d_1{l}_4{V}_2^2{n}_2\hfill \\ +d_1{l}_5{V}_2{n}_2^2+d_1{l}_6{n}_2^2+d_1{l}_7{n}_2^3-d_1\beta V_2^2-\beta k{V}_2{n}_2+d_1{c}_1+d_2{c}_0+d_{1}H(V_2,n_2)).\hfill \end{array}$

Taking $c_2=-\frac{1}{k}(d_1{c}_1+d_2{c}_0)$, $\overline{H(V_2,n_2)}=-\frac{d_{1}H(V_2,n_2)}{k}$, then (19) can become

(20)$\left\{\begin{array}{c}Q(V_2,n_2)=l_1{V}_2^2+l_2{V}_2^3+l_3{V}_2{n}_2+l_4{V}_2^2{n}_2+l_5{V}_2{n}_2^2+l_6{n}_2^2+l_7{n}_2^3+c_1+\overline{H(V_2,n_2)},\hfill \\ \\ R(V_2,n_2)=-\frac{(d_1{l}_1-d_1\beta )}{k}{V}_2^2-\frac{d_1{l}_2}{k}{V}_2^3-\frac{(d_1{l}_3-k\beta )}{k}{V}_2{n}_2-\frac{d_1{l}_4}{k}{V}_2^2{n}_2-\frac{d_1{l}_5}{k}{V}_2{n}_2^2-\frac{d_1{l}_6}{k}{n}_2^2-\frac{d_1{l}_7}{k}{n}_2^3.\hfill \end{array}\right.$

According to the method in [41,42] and Equation (20), we have

(21)$\begin{array}{c}{g}_{11}=\frac{1}{4}(\frac{\partial Q^2}{\partial {V_2}^2}+\frac{\partial Q^2}{\partial {n_2}^2}+i(\frac{\partial R^2}{\partial {V_2}^2}+\frac{\partial R^2}{\partial {n_2}^2}))\hfill \\ \\ =\frac{1}{4}(l_1+l_6-\frac{(d_1{l}_1-d_1\beta +d_1{l}_6)}{k}i),\hfill \end{array}$

(22)$\begin{array}{c}{g}_{02}=\frac{1}{4}(\frac{\partial Q^2}{\partial {V_2}^2}-\frac{\partial Q^2}{\partial {n_2}^2}-2\frac{\partial R^2}{\partial V_2\partial n_2}+i(\frac{\partial R^2}{\partial {V_2}^2}-\frac{\partial R^2}{\partial {n_2}^2}+2\frac{\partial Q^2}{\partial V_2\partial n_2}))\hfill \\ \\ =\frac{1}{4}(l_1-l_6+2\frac{(d_1{l}_3-k\beta )}{k}-\frac{(d_1{l}_1-d_1\beta -d_1{l}_6+2k{l}_3)}{k}i),\hfill \end{array}$

(23)$\begin{array}{c}{g}_{20}=\frac{1}{4}(\frac{\partial Q^2}{\partial {V_2}^2}-\frac{\partial Q^2}{\partial {n_2}^2}+2\frac{\partial R^2}{\partial V_2\partial n_2}+i(\frac{\partial R^2}{\partial {V_2}^2}-\frac{\partial R^2}{\partial {n_2}^2}-2\frac{\partial Q^2}{\partial V_2\partial n_2}))\hfill \\ \\ =\frac{1}{4}(l_1-l_6-2\frac{(d_1{l}_3-k\beta )}{k}-\frac{(d_1{l}_1-d_1\beta -d_1{l}_6-2k{l}_3)}{k}i),\hfill \end{array}$

(24)$\begin{array}{c}{G}_{21}=\frac{1}{8}(\frac{\partial Q^3}{\partial {V_2}^3}+\frac{\partial Q^3}{\partial V_2\partial {n_2}^2}+\frac{\partial R^3}{\partial V_2^2\partial n_2}+\frac{\partial R^3}{\partial {n_2}^3}+i(\frac{\partial R^3}{\partial {V_2}^3}+\frac{\partial R^3}{\partial V_2\partial {n_2}^2}-\frac{\partial Q^3}{\partial V_2^2\partial n_2}-\frac{\partial Q^3}{\partial {n_2}^3}))\hfill \\ \\ =\frac{1}{8}(l_2+l_5-\frac{d_1{l}_4}{k}-\frac{d_1{l}_7}{k}-\frac{(d_1{l}_2+d_1{l}_5-k{l}_4-k{l}_7)}{k}i).\hfill \end{array}$

Setting $g_{21}=G_{21}$, we obtain

$M_1\left(0\right)=\frac{i}{2k}(g_{20}{g}_{11}-2\mid g_{11}{\mid }^2-\frac{1}{3}\mid g_{02}{\mid }^2)+\frac{1}{2}{g}_{21}.$

Thus, we can compute the following quantities:

(25)$\begin{array}{c}{\mu }_2=-\frac{Re\left(M_1\left(0\right)\right)}{{\alpha }^{\prime }\left(0\right)}\hfill \\ =-\frac{1}{\left({(d_1+d_4)}^{\prime }\right)8k^2}[\frac{1}{2}(3l_1+2\beta )(l_1-\beta )-(2k+\frac{d_1^2}{k})(l_1{l}_3+l_3{l}_6)+\frac{d_1^2\beta l_3}{k}\hfill \\ -\frac{3}{2}{d}_1{l}_6^2+l_1{k}^2+l_5{k}^2-(d_1{l}_4+d_1{l}_7)],\hfill \end{array}$

(26)$\begin{array}{c}{\gamma }_2=2Re\left(M_1\left(0\right)\right)=-{(d_1+d_4)}^{\prime }{\mu }_2\hfill \\ \\ =\frac{1}{8k^2}[\frac{1}{2}(3l_1+2\beta )(l_1-\beta )-(2k+\frac{d_1^2}{k})(l_1{l}_3+l_3{l}_6)+\frac{d_1^2\beta l_3}{k}-\frac{3}{2}{d}_1{l}_6^2+l_1{k}^2+l_5{k}^2-(d_1{l}_4+d_1{l}_7)],\hfill \end{array}$

(27)$\begin{array}{c}{\tau }_2=-\frac{Im\left(M_1\left(0\right)\right)+{\mu }_2{\omega }^{\prime }\left(0\right)}{{\omega }_0}\hfill \\ \\ =\frac{1}{16k^2}[(\frac{2}{3}{l}_1^2+\frac{5}{3}{l}_6^2+\frac{2}{3}{\beta }^2+4l_3^2+\frac{5}{3}{l}_1{l}_6-\frac{1}{3}\beta l_1-\frac{5}{3}\beta l_6+d_1{l}_2+d_1{l}_5-k{l}_4-k{l}_7)\hfill \\ \\ +\frac{1}{k}(\frac{5}{3}{d}_1{l}_3{l}_6-\frac{4}{3}{d}_1{l}_3\beta -\frac{17}{3}{d}_1{l}_1{l}_3+6l_1{l}_3\beta +2d_1{l}_3{l}_6)\hfill \\ \\ +\frac{1}{k^2}{(\frac{2}{3}{d}_1^2{l}_3^2+\frac{13}{6}{d}_1^2{l}_1^2-\frac{13}{3}{d}_1^2{l}_1\beta -\frac{5}{3}{d}_1{l}_6\beta +d_1^2{\beta }^2+d_1^2{l}_6)}^2\hfill \\ \\ +\frac{2k^{\prime }}{\left({(d_1+d_4)}^{\prime }\right)k}(\frac{1}{2}(3l_1+2\beta )(l_1-\beta )-(2k+\frac{d_1^2}{k})(l_1{l}_3+l_3{l}_6)+\frac{d_1^2\beta l_3}{k}-\frac{3}{2}{d}_1{l}_6^2+l_1{k}^2\hfill \\ \\ +l_5{k}^2-(d_1{l}_4+d_1{l}_7)\left)\right].\hfill \end{array}$

If ${\mu }_2>0(<0)$, the direction of the bifurcation of model (4) is supercritical (subcritical); if ${\gamma }_2>0(<0)$, the bifurcating periodic solutions of model (4) are orbitally stable (unstable); if ${\tau }_2>0(<0)$, the periods of the bifurcating periodic solutions increase (decrease). We complete the proof. □

Remark 3.

It is a challenging and meaningful problem to obtain the conditions for the supercritical and subcritical Hopf bifurcations of the original H-H model (1) as a set of four coupled, highly nonlinear ordinary differential equations. The dynamical behaviors of the simplified H-H model (4) are consistent with that of the original H-H model (1). In addition, Theorem 2 presents conditions for the supercritical and subcritical Hopf bifurcations of model (4). Thus, we can gain a better insight into the direction of Hopf bifurcations in the original H-H model (1) by Theorem 2.

Remark 4.

Generally, the biological significance of the coefficients $C_M$, $\overline{g_{Na}}$, $\overline{g_K}$, $\overline{g_l}$, $V_{N_a}$, $V_K$, and $V_l$ need to be considered in the numerical value choices. If we choose different values for these coefficients, the equilibrium points of model (4) may change, and the Hopf bifurcations about the bifurcation parameters $I,V_E$ may also alter. As long as the coefficients and bifurcation parameters satisfy the conditions of Theorems 1 and 2, we can get the corresponding conclusions on Hopf bifurcation. Theorems 1 and 2 are general results for the Hopf bifurcations of model (4) with the parameters $I,V_E$.

4. Numerical Analysis and Discussion

In this section, some numerical results for the conditions of supercritical and subcritical Hopf bifurcations of model (4) are given to illustrate the above theoretical results.

According to [3,15,33,39,40], we choose the coefficients $C_M=1.0$ $\mathsf{\mu }$F/cm${}^2$, ${\overline{g}}_{Na}=120$ ms/cm${}^2$, ${\overline{g}}_K=36.0$ ms/cm${}^2$, ${\overline{g}}_l=0.3$ ms/cm${}^2$, $V_{N_a}=115.0$ mV, $V_K=-12$ mV, $V_l=10.599$ mV, ${\alpha }_1=0.057$, ${\alpha }_2=-0.0037$, ${\beta }_1=0.125$ and ${\beta }_2=0.0015$. Submitting these coefficients into model (4), we can obtain a more specific model as follows:

(28)$\begin{array}{c}\frac{dV}{dt}=I-120{\left[m_{\infty }\left(V\right)\right]}^3(0.8-n)(V+V_E-115)-36n^4(V+V_E+12)-0.3(V+V_E-10.599),\hfill \\ \frac{dn}{dt}=(0.057+0.0037V)(1-n)-(0.125-0.0015V)n,\hfill \end{array}$

where $m_{\infty }\left(V\right)=\frac{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]}{0.1(25.0-V)/\left[exp\left(\right(25.0-V)/10)-1.0\right]+4.0exp(-V/18.0)}.$

The solutions of model (28) are solved by the fourth-order Runge–Kutta algorithm. In addition, the bifurcation diagrams for the parameters I and $V_E$ are given by MATCONT which is a MATLAB continuation package for research of parameterized ODE systems [43]. Then, the dynamical behaviors of model (28) are analyzed with respect to the parameters I and $V_E$.

4.1. Numerical Analysis

The eigenvalues of the Jacobian of model (28), quadratic inequality, and T can be computed from the root of the equilibrium point equation which is obtained with particular values of the background stimulus I and $V_E$. All the values are listed in Table 1. In addition, the last column in Table 1 describes the stability of the fixed points with S (denoting stable) and U (denoting unstable).

Examining Table 1, model (28) has a Hopf bifurcation as the parameter I ($V_E=0$) increases through $I=7.5036$. The equilibrium point is $(V_0,n_0)=(4.6441,0.3859)$. We compute $a_0=0.00073046$, $a_1=0.00024228$, $a_2=0.00003813$, $a_3=0.00000372$, $k_1=49.6920$, $k_2=-5483.8$, $k_3=-13243$, $k_4=599.1876$, $d_1=0.1921$, $d_2=-147.4094$, $d_3=0.0029$, $d_4=-0.1922$, ${\alpha }^{\prime }\left(0\right)=\frac{{(d_1+d_4)}^{\prime }}{2}=0.5374$, ${\omega }^{\prime }\left(0\right)=\frac{{(d_1{d}_4-d_2{d}_3)}^{\prime }}{2k}=271.0423$, $l_1=0.1811$, $l_2=0.0178$, $l_3=-0.0634$, $l_4=-0.0017$, $l_5=-0.00064$, $l_6=-0.0096$, $l_7=-0.000013$, ${\mu }_2=-0.0401$, ${\gamma }_2=0.0431$, ${\tau }_2=17.3983$. In light of Theorem 2, the fact that ${\mu }_2<0$ implies that the Hopf bifurcation is subcritical. This means that the equilibrium point of model (28) is stable when $I<7.5036$, and the equilibrium loses its stability. A Hopf bifurcation occurs when I increases past $I=7.5036$.

For $I=6$, the equilibrium point of model (28) is stable as shown in Figure 2a,b. When $I=7.41$, there is a center in model (28) as shown in Figure 2c,d. For $I=7.51$, the equilibrium points of model (28) are unstable, and the phase trajectory converges to the stable limit cycle, as shown in Figure 2d,e. Overall, Figure 2 shows that the stable equilibrium solutions of model (28) come from the unstable equilibrium solutions and the stable limit cycles after a Hopf bifurcation about the parameter I. According to [44], the direction of the Hopf bifurcation is subcritical near $I=7.5036$.

Similarly, examining Table 1, model (28) has a Hopf bifurcation as the parameter $V_E$ ($I=0$) increases through $V_E=-6.2825$. We obtain the equilibrium point $(V_0,n_0)=(4.3190,0.3811)$, and compute $k_1=49.6920$, $k_2=-5879.5$, $k_3=-14036$, $k_4=361.3140$, $d_1=0.1914$, $d_2=-89.1968$, $d_3=0.0029$, $d_4=-0.1915$, ${\alpha }^{\prime }\left(0\right)=\frac{{(d_1+d_4)}^{\prime }}{2}=-0.058$, ${\omega }^{\prime }\left(0\right)=\frac{{(d_1{d}_4-d_2{d}_3)}^{\prime }}{2k}=0.0266$, $l_1=0.1655$, $l_2=0.0172$, $l_3=-0.0979$, $l_4=-0.0017$, $l_5=-0.00097$, $l_6=-0.0088$, $l_7=-0.000016$. We have ${\mu }_2=5.1$, ${\gamma }_2=0.0561$, ${\tau }_2=3.2505$. Because ${\mu }_2>0$, the Hopf bifurcation is supercritical, which means that the equilibrium point of model (28) is stable when $V_E>-6.2825$ and loses its stability through a Hopf bifurcation as $V_E$ decreases past $V_E=-6.2825$.

For $V_E=-6.284$, the $V-n$ phase plane trajectories are shown in Figure 3a,b. From Figure 3a,b, the equilibrium point is unstable, and the limit cycle is stable. When $V_E=-6.265$, model (28) has a center as shown in Figure 3c,d. When $V_E=-5.6$, the equilibrium point of model (28) is stable, as shown in Figure 3e,f. Overall, Figure 3 shows that the unstable equilibrium solution and stable limit cycles of the model (28) come from the stable equilibrium following a Hopf bifurcation near $V_E=-6.2825$. According to [44], the direction of this Hopf bifurcation is supercritical near $V_E=-6.2825$.

4.2. Discussion

In this subsection, we mainly discuss the bifurcation diagram for the parameters I and $V_E$, respectively. Moreover, we analyze the dynamical behaviors of model (28) about the parameters I and $V_E$. Especially, we focus on the Hopf bifurcation.

The bifurcation diagram for the parameter I is shown in Figure 4. The $I-V_{max}$, $I-V_{min}$ curves and equilibrium points of model (28) are shown in Figure 4a. A locally enlarged view of Figure 4a is shown in Figure 4b. The three points $A,B,E$ in Figure 4a are explained now. Point $A(I=5.26475)$ is a canard point (for details of the canards in model (28), see [39,40]). Points $B(I=7.5036)$ and $E(I=355)$ are Hopf bifurcation points. As these two points exhibit similar behavior, we only consider the point B. In Figure 4b, there are stable equilibrium points, unstable limit cycles, and stable limit cycles between points A and B.

To facilitate analysis of the Hopf bifurcation at point B, we take $I_A\approx 5.26475$ and $I_B\approx 7.5036$. When $0<I<I_A$, the equilibrium points of model (28) are stable. When $I_A<I<I_B$, the equilibrium points of model (28) are stable, and the limit cycles near the equilibrium points are unstable. However, model (28) has stable limit cycles that are totally different from the limit cycles near the equilibrium points. The reason is that model (28) has canard points. When $I_B<I<I_E$, model (28) has unstable equilibrium points and stable limit cycles. The stable equilibrium points become unstable equilibrium points, and the unstable limit vanishes when model (28) passes a Hopf bifurcation near $I_B\approx 7.5036$. This dynamic behavior is totally consistent with the original H-H model.

Take $I=5.5$, $I=6$, $I=6.5$ and $I=7$, respectively, then they lie in $I_A<I<I_B$. Choose the initial condition near the equilibrium points, the solutions of model (28) have converged to the equilibrium points that are shown in Figure 5a. When we choose the initial condition near the stable limit cycles, the solutions of model (28) have converged to the equilibrium points. The stable equilibrium points and stable limit cycles are shown in Figure 5b. Because the dynamic behaviors of model (28) about the parameters $I=5.5$, $I=6$, $I=6.5$ and $I=7$ are uniform, we choose $I=6$ to show the unstable limit cycles in model (28). The stable equilibrium points, unstable limit cycles and stable limit cycles are shown in Figure 5c for $I=6$. A locally enlarged view of Figure 5c is shown in Figure 5d. According to Figure 5c,d, the stable equilibrium points, unstable limit cycles, and stable limit cycles coexist in model (28). So, the model (28) has stable equilibrium points, unstable limit cycles, and stable limit cycles for $I_A<I<I_B$.

Similarly, the bifurcation diagram for parameter $V_E$ is shown in Figure 6. The $V_E-V_{max}$, $V_E-V_{min}$ curves, and equilibrium points of model (28) are shown in Figure 6a. A locally enlarged view of Figure 6a is shown in Figure 6b. The point $A(V_E=-5.4456)$ in Figure 6a is a canard point. Points $B(V_E=-6.2825)$ and $E(V_E=-37.2)$ are Hopf bifurcation points. Again, we only consider the point B here. In Figure 6b, there are stable equilibrium points, unstable limit cycles, and stable limit cycles in the region between points A and B.

Consider $V_E\approx -5.4456$ and $V_E\approx -6.2825$, respectively. When $-5.4456<V_E<0$, the equilibrium points of model (28) are stable. When $-6.2825<V_E<-5.4456$, the equilibrium points of model (28) are stable, and the limit cycles near the equilibrium points are unstable. However, model (28) has stable limit cycles that are totally different from the limit cycles near the equilibrium points. The reason is that model (28) has canard points. When $-37.2<V_E<-6.2825$, model (28) has unstable equilibrium points and stable limit cycles. The stable equilibrium points become unstable equilibrium points, and the unstable limit vanish when model (28) passes a Hopf bifurcation near $V_E\approx -6.2825$. The dynamic behaviors of model (28) about the parameters I and $V_E$ are summarized in Table 2.

5. Conclusions

This paper has examined the Hopf bifurcation of a two-dimensional simplified H-H model with two parameters. Based on theoretical analysis and numerical simulations, a sufficient and necessary condition for the existence of Hopf bifurcations was derived. The Hopf bifurcations of the two-dimensional simplified H-H model were discussed using normal form theory. The direction of the Hopf bifurcations and the stability of bifurcating periodic solutions were also analyzed. It is significant and meaningful to obtain the conditions for supercritical and subcritical Hopf bifurcations of the original H-H model as a set of four coupled and highly nonlinear ordinary differential equations. The dynamical behavior of the considered two-dimensional simplified H-H model is consistent with that of the original H-H model, so the direction of the Hopf bifurcations in the original H-H model is better interpreted by the obtained results. In future, the conditions of supercritical and subcritical Hopf bifurcations in the H-H model with multiple parameters will be discussed, and codimension-two bifurcations in the H-H model will be investigated.

Author Contributions

Conceptualization, methodology, software, H.W.; writing original draft preparation, validation, S.W.; writing review and editing, Y.G.; methodology, supervision, Y.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the respected reviewers for their kind comments and the editorial office for their advice. The authors also thank Stuart Jenkinson, for editing the English text of a draft of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Circuit description of the H-H model.

Figure 2. The convergent behaviors of model (28) about parameter I. (a) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (b) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (c) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (d) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (e) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (f) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.].

Figure 3. The convergent behaviors of model (28) about parameter [Forumla omitted. See PDF.]. (a) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (b) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (c) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (d) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (e) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (f) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.].

Figure 3. The convergent behaviors of model (28) about parameter [Forumla omitted. See PDF.]. (a) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (b) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (c) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (d) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (e) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.]. (f) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.].

Figure 4. Bifurcation diagram for parameter I. (a) The I–Vmax, I–Vmin curves and equilibrium points of model (28). (b) A locally enlarged view of (a) near Hopf point B.

Figure 5. The convergent behaviors of model (28) about parameter I. (a) Response curve of trans-membrane voltage V for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. (b) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. (c) Trajectories in [Forumla omitted. See PDF.] phase plane for [Forumla omitted. See PDF.]. (d) A locally enlarged view of (c) for [Forumla omitted. See PDF.].

Figure 6. Bifurcation diagram for parameter [Forumla omitted. See PDF.]. (a) The [Forumla omitted. See PDF.]–Vmax, [Forumla omitted. See PDF.]–Vmin curves and equilibrium points of model (28). (b) A locally enlarged view of (a) near Hopf point B.

Table 1

Eigenvalues of the Jacobian, quadratic inequality, and T.

Parameter I	${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	Inequality	T
5.264	−0.107653 + 0.577980i	−0.107653 − 0.577980i	−0.304819	−0.215306	S
5.265	−0.107612 + 0.578004i	−0.107612 − 0.578004i	−0.304854	−0.215225	S
5.266	−0.107571 + 0.578028i	−0.107571 − 0.578028i	−0.304888	−0.215143	S
⋯	⋯	⋯	⋯	⋯	⋯
7.503	−0.000038 + 0.619185i	−0.000038 − 0.619185i	−0.383375	−7.624340 × 10${}^{-5}$	S
7.504	+0.000017 + 0.619196i	+0.000017 − 0.619196i	−0.383411	+3.542764 × 10${}^{-5}$	U
Parameter ${\mathit{V}}_{\mathit{E}}$	${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	Inequality	$\mathit{T}$
−6.283	0.000208 + 0.467422i	0.000208 − 0.467422i	−0.287757	4.168339 × 10${}^{-4}$	U
−6.282	+0.000097 + 0.467435i	+0.000097 − 0.467435i	−0.287722	+1.958441 × 10${}^{-4}$	U
−6.281	−0.000012 + 0.467449i	−0.000012 − 0.467449i	−0.287687	−2.503515 × 10${}^{-5}$	S
⋯	⋯	⋯	⋯	⋯	⋯
−5.446	−0.076029 + 0.469630i	−0.076029 − 0.469630i	−0.262058	−0.152059	S
−5.445	−0.076104 + 0.469625i	−0.076104 − 0.469625i	−0.262031	−0.152209	S

Table 2

The dynamic behaviors of model (28) about the parameters I and $V_E$.

Parameter I	$0<\mathit{I}<5.26475$	$5.26475<\mathit{I}<7.5306$	$7.5306<\mathit{I}<355$	$\mathit{I}>355$
Dynamical Behaviors	stable equilibrium points	stable equilibrium points	unstable equilibrium points	stable equilibrium points
		stable limit cycles	stable limit cycles
		unstable limit cycles
Parameter ${\mathit{V}}_{\mathit{E}}$	$\mathbf{0}>{\mathit{V}}_{\mathit{E}}>-\mathbf{5}.\mathbf{4456}$	$-\mathbf{5}.\mathbf{4456}>{\mathit{V}}_{\mathit{E}}>-\mathbf{6}.\mathbf{2825}$	$-\mathbf{6}.\mathbf{2825}>{\mathit{V}}_{\mathit{E}}>-\mathbf{37}.\mathbf{2}$	$-\mathbf{37}.\mathbf{2}>{\mathit{V}}_{\mathit{E}}$
Dynamical Behaviors	stable equilibrium points	stable equilibrium points	unstable equilibrium points	stable equilibrium points
		stable limit cycles	stable limit cycles
		unstable limit cycles

References

1. Hodgkin, A.L.; Huxley, A.F. A quantitative description of membrane and its application to conduction and excitation in nerve. J. Physiol.; 1952; 117, pp. 500-544. [DOI: https://dx.doi.org/10.1113/jphysiol.1952.sp004764] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12991237]

2. Izhikevich, E.M. Which model to use for cortical spiking neurons?. IEEE Trans. Neural Netw.; 2004; 15, pp. 1063-1070. [DOI: https://dx.doi.org/10.1109/TNN.2004.832719] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15484883]

3. Izhikevich, E.M. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting; The MIT Press: Cambridge, MA, USA, 2010.

4. Liu, Y.; Iu, H.H.C.; Qian, Y.H. Implementation of Hodgkin-Huxley neuron model with the novel memristive oscillator. IEEE Trans. Circuits Syst. II Express Briefs; 2021; 68, pp. 2982-2986. [DOI: https://dx.doi.org/10.1109/TCSII.2021.3066471]

5. Chua, L. Hodgkin-Huxley equations implies edge of chaos kernel. Jpn. J. Appl. Phys.; 2022; 61, SM0805. [DOI: https://dx.doi.org/10.35848/1347-4065/ac64e1]

6. Valle, J.; Madureira, A.L. Parameter identification problem in the Hodgkin-Huxley model. Neural Comput.; 2022; 34, pp. 939-970. [DOI: https://dx.doi.org/10.1162/neco_a_01487]

7. Baysal, V.; Yilmaz, E. Chaotic signal induced delay decay in Hodgkin-Huxley neuron. Appl. Math. Comput.; 2021; 411, 126540. [DOI: https://dx.doi.org/10.1016/j.amc.2021.126540]

8. Haghiri, S.; Naderi, A.; Ghanbari, B.; Ahmadi, A. High speed and low digital resources implementation of Hodgkin-Huxley neuronal model using base-2 functions. IEEE Trans. Circuits Syst. I Regul. Pap.; 2021; 68, pp. 275-287. [DOI: https://dx.doi.org/10.1109/TCSI.2020.3026076]

9. Chen, Z.D.; Raman, B.; Stern, A. Structure-preserving numerical integrators for Hodgkin-Huxley-Type systems. SIAM J. Sci. Comput.; 2020; 42, pp. 273-298. [DOI: https://dx.doi.org/10.1137/18M123390X]

10. FitzHugh, R. Impulses and physiological states in theoretical models of nerve membrane. J. Biophys.; 1961; 1, pp. 445-466. [DOI: https://dx.doi.org/10.1016/S0006-3495(61)86902-6]

11. Nagumo, J.; Arimoto, S.; Yoshizawa, S. An active pulse transmission line simulating nerve axon. Proc. IRE; 1962; 50, pp. 2061-2070. [DOI: https://dx.doi.org/10.1109/JRPROC.1962.288235]

12. Morris, C.; Lecar, H. Voltage Oscillations in the barnacle giant muscle fiber. J. Biophys.; 1981; 35, pp. 193-213. [DOI: https://dx.doi.org/10.1016/S0006-3495(81)84782-0]

13. Chay, T.R. Electrical bursting and luminal calcium oscillation in excitable cell models. Biol. Cybern.; 1996; 75, pp. 419-431. [DOI: https://dx.doi.org/10.1007/s004220050307]

14. Hindmarsh, J.L.; Rose, R.M. A model of neuronal bursting using three coupled first order differential equations. Proc. R. Soc. Lond. Ser. B; 1984; 221, pp. 87-102.

15. Moehlis, J. Canards for a reduction of the Hodgkin-Huxley equations. J. Math. Biol.; 2006; 52, pp. 141-153. [DOI: https://dx.doi.org/10.1007/s00285-005-0347-1]

16. Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields; Springer: Berlin/Heidelberg, Germany, 1983.

17. Rabinovich, M.I.; Varona, P.; Selverston, A.I.; Abarbanel, H.D.I. Dynamical principles in neuroscience. Rev. Mod. Phys.; 2006; 78, pp. 1213-1265. [DOI: https://dx.doi.org/10.1103/RevModPhys.78.1213]

18. Fukai, H.; Doi, S.; Nomura, T.; Sato, S. Hopf bifurcations in multiple parameter space of the Hodgkin-Huxley equations. I. Global organization of bistable periodic solutions. Biol. Cybern.; 2000; 82, pp. 215-222. [DOI: https://dx.doi.org/10.1007/s004220050021]

19. Fukai, H.; Doi, S.; Nomura, T.; Sato, S. Hopf bifurcations in multiple parameter space of the Hodgkin-Huxley equations. II. Singularity theoretic approach and highly degenerate bifurcations. Biol. Cybern.; 2000; 82, pp. 223-229. [DOI: https://dx.doi.org/10.1007/s004220050022]

20. Guckenheimer, J.; Oliva, A. Chaos in the Hodgkin-Huxley model. SIAM J. Appl. Dyn. Syst.; 2002; 1, pp. 105-114. [DOI: https://dx.doi.org/10.1137/S1111111101394040]

21. Drover, J.; Rubin, J.; Su, J.; Ermentrout, B. Analysis of a canard mechanism by which excitatory synaptic coupling can synchronize neurons at low firing frequencies. SIAM J. Appl. Math.; 2005; 65, pp. 69-92. [DOI: https://dx.doi.org/10.1137/S0036139903431233]

22. Calitoiu, D.; Oommen, B.J.; Nussbaum, D. Spikes annihilation in the Hodgkin-Huxley neuron. Biol. Cybern.; 2008; 98, pp. 239-257. [DOI: https://dx.doi.org/10.1007/s00422-007-0207-8]

23. Kang, Y.M.; Chen, Y.Q.; Fu, Y.X.; Wang, Z.L.; Chen, G.R. Formation of spiral wave in Hodgkin-Huxley neuron networks with Gamma-distributed synaptic input. Commun. Nonlinear Sci. Numer. Simul.; 2020; 83, 105112. [DOI: https://dx.doi.org/10.1016/j.cnsns.2019.105112]

24. Erhardt, A.H. Bifurcation analysis of a certain Hodgkin-Huxley model depending on multiple bifurcation parameters. Mathematics; 2018; 6, 103. [DOI: https://dx.doi.org/10.3390/math6060103]

25. Li, S.L.; Lv, W.P.; Chen, Z.Y.; Xue, M.; Bi, Q.S. Slow-fast dynamics behaviors under the comprehensive effect of rest spike bistability and timescale difference in a Filippov slow-fast modified Chua’s circuit Model. Mathematics; 2022; 10, 4606. [DOI: https://dx.doi.org/10.3390/math10234606]

26. Berzhanskaya, J.; Gorchetchnikov, A.; Schiff, S.J. Switching between gamma and theta: Dynamic network control using subthreshold electric fields. Neurocomputing; 2007; 70, pp. 2091-2095. [DOI: https://dx.doi.org/10.1016/j.neucom.2006.10.124]

27. Park, E.H.; So, P.; Barreto, E.; Gluckman, B.J.; Schiff, S.J. Electric field modulation of synchronization in neuronal networks. Neurocomputing; 2002; 52–54, pp. 169-175. [DOI: https://dx.doi.org/10.1016/S0925-2312(02)00820-2]

28. Bikson, M.; Inoue, M.; Akiyama, H.; Deans, J.K.; Fox, J.E.; Miyakawa, H.; Jefferys, G.R. Effects of uniform extracellular DC electric fields on excitability in rat hippocampal slices in vitro. J. Physiol.; 2004; 557, pp. 175-190. [DOI: https://dx.doi.org/10.1113/jphysiol.2003.055772] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14978199]

29. Bujan, A.F.; Aertsen, A.; Kumar, A. Role of input correlations in shaping the variability and noise correlations of evoked activity in the neocortex. J. Neurosci.; 2015; 35, pp. 8611-8625. [DOI: https://dx.doi.org/10.1523/JNEUROSCI.4536-14.2015] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26041927]

30. Yu, H.; Guo, X.; Wang, J. Stochastic resonance enhancement of small-world neural networks by hybrid synapses and time delay. Commun. Nonlinear Sci. Numer. Simul.; 2017; 42, pp. 532-544. [DOI: https://dx.doi.org/10.1016/j.cnsns.2016.06.021]

31. Wang, J.; Chen, L.Q.; Fei, X.Y. Bifurcation control of the Hodgkin-Huxley equations. Chaos Solitons Fractals; 2007; 33, pp. 217-224. [DOI: https://dx.doi.org/10.1016/j.chaos.2006.01.035]

32. Che, Y.Q.; Wang, J.; Si, W.J.; Fei, X.Y. Phase-locking and chaos in a silent Hodgkin-Huxley neuron exposed to sinusoidal electric field. Chaos Solitons Fractals; 2009; 39, pp. 454-462. [DOI: https://dx.doi.org/10.1016/j.chaos.2007.03.008]

33. Che, Y.Q.; Wang, J.; Deng, B.; Wei, X.L.; Han, C.X. Bifurcations in the Hodgkin-Huxley model exposed to DC electric fields. Neurocomputing; 2012; 81, pp. 41-48. [DOI: https://dx.doi.org/10.1016/j.neucom.2011.11.019]

34. Xie, Y.; Chen, L.; Kang, Y.M.; Aihara, K. Controlling the onset of Hopf bifurcation in the Hodgkin-Huxley model. Phys. Rev. E; 2008; 77, 061921. [DOI: https://dx.doi.org/10.1103/PhysRevE.77.061921] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18643314]

35. Ding, L.; Hou, C. Stabilizing control of Hopf bifurcation in the Hodgkin-Huxley model via washout filter with linear control term. Nonlinear Dynam.; 2010; 60, pp. 131-139. [DOI: https://dx.doi.org/10.1007/s11071-009-9585-x]

36. Doruk, R. Control of Hopf bifurcations in Hodgkin-Huxley neurons by automatic temperature manipulation. Neuroquantology; 2018; 16, pp. 59-74. [DOI: https://dx.doi.org/10.14704/nq.2018.16.2.1140]

37. Zhao, Z.; Li, L.; Gu, H.; Gao, Y. Different dynamics of repetitive neural spiking induced by inhibitory and excitatory autapses near subcritical Hopf bifurcation. Nonlinear Dyn.; 2020; 99, pp. 1129-1154. [DOI: https://dx.doi.org/10.1007/s11071-019-05342-6]

38. Li, L.; Zhao, Z. White-noise-induced double coherence resonances in reduced Hodgkin-Huxley neuron model near subcritical Hopf bifurcation. Phys. Rev. E; 2022; 105, 034408. [DOI: https://dx.doi.org/10.1103/PhysRevE.105.034408] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35428043]

39. Wang, H.; Yu, Y.G.; Wang, S.; Yu, J.Z. Bifurcation analysis of a two-dimensional simplified Hodgkin-Huxley model exposed to external electric fields. Neural Comput. Appl.; 2014; 24, pp. 37-44. [DOI: https://dx.doi.org/10.1007/s00521-013-1462-3]

40. Wang, H.; Yu, Y.G.; Zhao, R.; Wang, S. Two-parameter bifurcation in a two-dimensional simplified Hodgkin-Huxley model. Commun. Nonlinear Sci. Numer. Simul.; 2013; 18, pp. 184-193. [DOI: https://dx.doi.org/10.1016/j.cnsns.2012.06.022]

41. Hassard, B.D.; Kazarinoff, N.D.; Wan, Y.H. Theory and Applications of Hopf Biufrcation; Cambridge University Press: Cambridge, UK, 1981.

42. Yu, Y.G.; Zhang, S.C. Hopf bifurcation in the Lü system. Chaos Solitons Fractals; 2003; 17, pp. 901-906. [DOI: https://dx.doi.org/10.1016/S0960-0779(02)00573-8]

43. Dhooge, A.; Govaerts, W.; Kuznetsov, Y.A. MATCONT: A Matlab package for numerical bifurcation analysis of ODEs. ACM Trans. Math. Softw.; 2003; 29, pp. 141-164. [DOI: https://dx.doi.org/10.1145/779359.779362]

44. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory; 3rd ed. Springer: New York, NY, USA, 2004.