Effects of Maximal Sodium and Potassium Conductance on the Stability of Hodgkin-Huxley Model

Abstract

Hodgkin-Huxley (HH) equation is the first cell computing model in the world and pioneered the use of model to study electrophysiological problems. The model consists of four differential equations which are based on the experimental data of ion channels. Maximal conductance is an important characteristic of different channels. In this study, mathematical method is used to investigate the importance of maximal sodium conductance 𝑔_Na and maximal potassium conductance 𝑔_K. Applying stability theory, and taking 𝑔_Na and 𝑔_K as variables, we analyze the stability and bifurcations of the model. Bifurcations are found when the variables change, and bifurcation points and boundary are also calculated. There is only one bifurcation point when 𝑔_Na is the variable, while there are two points when 𝑔_K is variable. The (𝑔_Na,  𝑔_K) plane is partitioned into two regions and the upper bifurcation boundary is similar to a line when both 𝑔_Na and 𝑔_K are variables. Numerical simulations illustrate the validity of the analysis. The results obtained could be helpful in studying relevant diseases caused by maximal conductance anomaly.

1. Introduction

Hodgkin-Huxley (HH) equation is created on the foundation of huge experimental data of sodium and potassium channels by Hodgkin and Huxley who are both excellent biology scientists and had long engaged in nerve conduction research. In about 1952, they took squid giant axon as experiment subject and continuously published four papers describing the electrical excitation of this kind of cell [1–4]. In their experiment, all the ion channels were divided into three types, sodium channel, potassium channel, and the others. Now we know there are many ion channels on the cell membrane, such as 𝐼_Na, 𝐼_Kr, 𝐼_Ks, 𝐼_NaCa, 𝐼_K1, 𝐼_CaL, 𝐼_Ca, 𝐼_to, 𝐼_NaK, 𝐼_NaL, and 𝐼_KATP [5–7]. However the discovery of sodium and potassium channels was marvelous at that time. Experimental data was obtained by voltage-clamp technology, while the patch-clamp technology is widely used at present. On this basis, a four-dimensional ordinary differential equation set, called HH model, was proposed, which was autonomous and contained intricate transcendental equations.

The work of Hodgkin and Huxley was recognized as excellent achievement and with significant contribution to the development of electrophysiology. It is the basis of the subsequent models of ion channels. Not only was the HH model consistent with the obtained experiment data accurately, but also it could precisely simulate the change of action potential. The model discovered the relationship of transmembrane potential and current and maximum conductance of ions. This made it possible to research the character of ion channel with mathematical methods. In 1960, Professor Nobel who pioneered the cardiac electrophysiology simulations applied HH model to myocardial cell and got the famous Purkinje fiber cell model [8], which was the first computing myocardial cell model. From then on, HH model was broadly applied to almost all kinds of cardiac cells such as atrial muscle cell model [9] and sinoatrial node cell model [10]. HH model laid the cornerstone of computing electrophysiology. Even today, a large part of electrophysiological models are created on the foundation of HH model. Verkerk’s sinoatrial model [11], Butters’s atrial model [12], O’Hara’s ventricular model [13], and Li’s Purkinje cell model [14], and so forth, all belong to HH model type.

Because of the importance of HH model, the stability has long attracted the researcher’s attention. Hassard et al. were the earlier researchers caring about the bifurcation phenomenon of HH model. And they indicated that bifurcation would occur at the equilibrium points when the external current 𝑖_ext changed which was injected into the neuron from microelectrode [15]. Stable and unstable solutions of the model with regard to 𝑖_ext were analyzed by Rinzel and Miller, and the influence of temperature was also discussed [16]. Two stable steady states were found by Aihara and Matsumoto [17]; when the two states existed, the bifurcation structure was complex, which included a stable limit cycle, two unstable equilibrium points, and one asymptotically stable equilibrium point. Guckenheimer and Labouriau investigated the influence of 𝑖_ext and potassium ion potential 𝑉_K on the bifurcations of the model [18]. Bedrov et al. gave the relationship between the numbers of negative slope regions and presented some results about the possible bifurcation giving rise to maximal sodium conductance 𝑔_Na and maximal potassium conductance 𝑔_K [19, 20]. Fukai and his fellows examined the structure of the model’s bifurcations produced by 𝑖_ext and one of the other parameters [21]. Taking leakage conductance 𝑔_𝑙 and sodium channel effective bias voltage 𝑉_𝑚 as parameters, Terada et al. analyzed bifurcation in Hodgkin-Huxley model for muscles of frogs [22].

Wang et al. [23–26] did a lot of research on the stability of HH model. They studied the bifurcations caused by leakage conductance 𝑔_𝑙 and sodium ions antielectromotive force when ELF external electric field was considered. The stability and bifurcation control were analyzed and controllers were designed. Bifurcation in HH model exposed to DC electric fields was investigated in detail.

Bifurcation means qualitative changes in the solution structure of a dynamic system when the parameters vary. From analyzing the bifurcation, we can get the effects of the parameters. Further, changing the corresponding parameters, we could make the solution into an ideal condition. Bifurcation is an important branch in mathematics and applied to much different field [27–29]. In recent years, it is also widely studied in electrophysiology. Indeed, there are many diseases having close relations with bifurcations, such as Parkinson’s, epilepsy, and pathological heart rhythms [30].

In the past, for HH model, external current 𝑖_ext and leakage conductance 𝑔_𝑙 have been most investigated, because they were easily measured. The sodium current is the contributor which leads to depolarization of the neuron while it is potassium current that plays the major role of repolarization. However, 𝑔_Na and 𝑔_K are seldom taken into consideration to analyze the stability of model, as the relevant data is not abundant. In this study, the effects of  𝑔_Na and 𝑔_K on the stability and bifurcations of the model will be discussed, respectively, and collectively. And we will give the critical points of 𝑔_Na and 𝑔_K when they play the role separately, and the critical boundaries in 𝑔_Na-𝑔_K plane will be provided when together. Simulation results demonstrate the validity of the theoretic analysis.

The rest of the paper is organized as follows. The HH equations are introduced in detail in Section 2. In Section 3 we analyze the effects of 𝑔_Na and 𝑔_K on the model and calculate the bifurcation points (line). Finally, discussion and conclusion are presented in Section 4.

2. Hodgkin-Huxley Equations

HH model was proposed on the foundation of ion channels. The electrophysiological activities of a cell are shown in Figure 1(a). The gray circle is membrane, which ensures orderly biochemical reaction. 𝐼_Na, 𝐼_K, and 𝐼_L are the ion currents corresponding to respective channels on the membrane. When an electrical stimulation makes the sodium channels open, a large number of Na⁺ flow inward, forming current 𝐼_Na, resulting in the rise of transmembrane potential. The open of potassium channels makes a large outflux of K⁺, creating the current 𝐼_K and the reduction of potential. The model is comprised of four autonomous ordinary differential equations to describe the electrophysiological activities of cell shown in Figure 1(a). In the model, membrane is taken as a constant capacitance and the ion channels are seen as variable resistances. Figure 1(b) shows the equivalent circuit in detail, in which 𝑅_Na=1/𝑔_Na, 𝑅_K=1/𝑔_K, and 𝑅_𝑙=1/𝑔_𝑙. 𝑅_Na and 𝑅_K vary with time.

(a) Electrophysiological activities of cell

(b) Electrical equivalent circuit of HH model

Figure 1 

The electrophysiological process and equivalent circuit of neuron.

The equations were obtained according to electrical formulas and experimental data, which are shown as follows: 𝑑𝑉=1𝑑𝑡𝐶_𝑀𝑖_ext−𝑔_Na𝑚³ℎ𝑉−𝑉_Na−𝑔_K𝑛⁴𝑉−𝑉_K−𝑔_𝑙𝑉−𝑉_𝑙,𝑑𝑚𝑑𝑡=𝛼_𝑚(𝑉)(1−𝑚)−𝛽_𝑚(𝑉)𝑚,𝑑ℎ𝑑𝑡=𝛼_ℎ(𝑉)(1−ℎ)−𝛽_ℎ(𝑉)ℎ,𝑑𝑛𝑑𝑡=𝛼_𝑛(𝑉)(1−𝑛)−𝛽_𝑛(𝑉)𝑛, (1) where 𝛼_𝑚(𝑉)=0.1(𝑉−25.0)[−()],𝛽1−exp𝑉−25.0/10_𝑚(𝑉)−𝑉=4.0exp,𝛼18.0_ℎ(𝑉)−𝑉=0.07exp,𝛽20.0_ℎ(𝑉)=1[−()],𝛼1+exp𝑉−30.0/10_𝑛(𝑉)=0.01(𝑉−10.0)[−()],𝛽1−exp𝑉−10.0/10_𝑛(𝑉)−𝑉=0.125exp.80.0 (2)

In these equations, 𝑉 is the transmembrane potential. 0≤𝑚≤1 and 0≤ℎ≤1 are the gating variables indicating activation and inactivation of sodium ion current, respectively. 0≤𝑛≤1 is the gating variable showing activation of potassium ion current. 𝑔_Na, 𝑔_K, and 𝑔_𝑙 represent the maximal conductance of corresponding currents. 𝐶_𝑚=1.0 μF/cm² is membrane capacitance. 𝑖_ext is the current injected into the neuron. In our paper, we suppose 𝑖_ext=0 and 𝑔_Na=120 mS/cm², 𝑔_K=36 mS/cm², and 𝑔_𝑙=0.3 mS/cm², which are the ideal experimental data.

3. Stability Analysis of HH Model

Stability is one of a model’s important properties. If the model is stable, it will reach a rest state at last. Otherwise, periodic phenomenon or chaos may appear. To analyze an ordinary differential system, equilibrium points are one of its most important aspects, which may be the final state of the system. Suppose (𝑉_∗, 𝑚_∗, ℎ_∗, 𝑛_∗) is the equilibrium points of the model. So it should make the right side of (1) equal to zero. That is, 𝑖_ext−𝑔_Na𝑚³_∗ℎ_∗𝑉_∗−𝑉_Na−𝑔_K𝑛⁴_∗𝑉_∗−𝑉_K−𝑔_𝑙𝑉_∗−𝑉_𝑙𝛼=0,_𝑚𝑉_∗1−𝑚_∗−𝛽_𝑚𝑉_∗𝑚_∗𝛼=0,_ℎ𝑉_∗1−ℎ_∗−𝛽_ℎ𝑉_∗ℎ_∗𝛼=0,_𝑛𝑉_∗1−𝑛_∗−𝛽_𝑛𝑉_∗𝑛_∗=0. (3) Then the linearization of (1) around the equilibrium could be obtained as follows: 𝑑𝑉𝑑𝑡=𝐽₁₁𝑉+𝐽₁₂𝑚+𝐽₁₃ℎ+𝐽₁₄𝑛,𝑑𝑚𝑑𝑡=𝐽₂₁𝑉+𝐽₂₂𝑚,𝑑ℎ𝑑𝑡=𝐽₃₁𝑉+𝐽₃₃ℎ,𝑑𝑛𝑑𝑡=𝐽₄₁𝑉+𝐽₄₄𝑛, (4) where 𝐽₁₁=−𝑔_Na𝑚³_∗ℎ_∗+𝑔_K𝑛⁴_∗+𝑔_𝑙𝐶_𝑀,𝐽₁₂3=−𝑔_Na𝑚²_∗ℎ_∗𝑉_∗−𝑉_Na𝐶_𝑀,𝐽₁₃=−𝑔_Na𝑚³_∗𝑉_∗−𝑉_Na𝐶_𝑀,𝐽₁₄4=−𝑔_K𝑛³_∗𝑉_∗−𝑉_K𝐶_𝑀,𝐽₂₁=2𝑚_∗9exp𝑉_∗−𝑉/18.00.1_∗−25.0exp−𝑉_∗−25.0/101−exp−𝑉_∗−25/10²−0.11−exp−𝑉_∗−25/101−𝑚_∗,𝐽₂₂=0.1𝑉_∗−25.01−exp−𝑉_∗−𝑉−25.0/10−4.0exp^∗,𝐽18.0₃₁𝑉=−7exp_∗/20.020001−ℎ_∗−exp−𝑉_∗−30.0/10ℎ_∗101+exp−𝑉_∗−30.0/10²,𝐽₃₃−𝑉=0.07exp^∗−120.01+exp−𝑉_∗,𝐽−30.0/10₄₁=𝑛_∗exp𝑉_∗/80.0+16401001−exp−𝑉_∗−−𝑉−10.0/100.01exp_∗−10.0/10𝑉_∗−10.0101−exp−𝑉_∗−10.0/10²1−𝑛_∗,𝐽₄₄=0.01𝑉_∗−10.01−exp−𝑉_∗−𝑉−10.0/10−0.125exp^∗.80.0 (5)

We can get the eigenmatrix of (4): 𝐽𝐽=¹¹𝐽₁₂𝐽₁₃𝐽₁₄𝐽₂₁𝐽₂₂𝐽00₃₁0𝐽₃₃0𝐽₄₁00𝐽₄₄ (6) and then the characteristic equation can be obtained: 𝜆⁴+𝑎𝜆³+𝑏𝜆²+𝑐𝜆+𝑑=0, (7) where 𝑎=−𝐽₁₁+𝐽₂₂+𝐽₃₃+𝐽₄₄,𝑏=𝐽₁₁𝐽₂₂+𝐽₃₃+𝐽₄₄+𝐽₂₂𝐽₃₃+𝐽₄₄+𝐽₃₃𝐽₄₄−𝐽₁₂𝐽₂₁−𝐽₁₃𝐽₃₁−𝐽₁₄𝐽₄₁,𝑐=𝐽₁₂𝐽₂₁𝐽₃₃+𝐽₄₄+𝐽₁₃𝐽₃₁𝐽₂₂+𝐽₄₄+𝐽₁₄𝐽₄₁𝐽₂₂+𝐽₃₃−𝐽₁₁𝐽₂₂𝐽₃₃+𝐽₄₄−𝐽₁₁+𝐽₂₂𝐽₃₃𝐽₄₄,𝑑=𝐽₁₁𝐽₂₂𝐽₃₃𝐽₄₄−𝐽₁₂𝐽₂₁𝐽₃₃𝐽₄₄−𝐽₁₃𝐽₂₂𝐽₃₁𝐽₄₄−𝐽₁₄𝐽₂₂𝐽₃₃𝐽₄₁. (8)

According to Routh-Hurwitz criterion, if 𝑎>0, 𝑎𝑏>𝑐, 𝑑>0,  𝑎𝑏𝑐>𝑐²+𝑎²𝑑, the real parts of all the roots are minus. Otherwise, all the real parts are not negative. In the following, we will analyze the stability of the model according to this criterion.

3.1. Effect of 𝑔_Na on the Stability

In this section, we will investigate the influence of 𝑔_Na on the equilibrium, stability, and bifurcations of the model. 𝑔_Na is taken as variable, and the other parameters are all kept with desired values. Because the desired value of 𝑔_Na is around 120 mS/cm², 𝑔_Na∈[0,500] is taken into consideration. When 𝑔_Na changes, making the right side of (1) equal to zero, corresponding 𝑉^∗ can be acquired. Taking Matlab as a tool, we could obtain the relationship between 𝑔_Na and the equilibrium point 𝑉^∗ shown in Figure 2.

Figure 2 

The relationship between 𝑔Na and 𝑉∗.

From Figure 2, we can see that 𝑉^∗ changes slowly when 𝑔_Na∈[0,300] and increases rapidly when 𝑔_Na∈[350,500]. This means that equilibrium points are sensitive to 𝑔_Na when 𝑔_Na∈[350,500]; a slight change of 𝑔_Na may lead the model to a totally different state even though model is still stable.

Applying bifurcation theory and using the method of bisection, we can get one bifurcation point 𝑔^∗_Na=212.648720656 when 𝑔_Na changes.

Substituting 𝑔^∗_Na into the original equation, we get the equilibrium 𝑉^∗, and then substituting both 𝑔^∗_Na and 𝑉^∗ into eigenmatrix of (4), we can gain the eigenvalues as follows: 𝜆₁𝜆=−4.9711711484,₂𝜆=−0.1259717048,₃=1.9×10⁻¹⁶𝜆−0.3798402483𝑖,₄=1.9×10⁻¹⁶+0.3798402483𝑖. (9)

Here, we regard 1.9×10⁻¹⁶ as 0. With the help of computer, we can get that all the real parts of 𝜆_𝑖 (𝑖=1,2,3,4) are negative when 𝑔_Na<𝑔^∗_Na. And there exist positive real parts when 𝑔_Na>𝑔^∗_Na. According to the stability theory, the system is stable around equilibrium when 𝑔_Na∈[0,𝑔^∗_Na) and it is unstable when 𝑔_Na∈(𝑔^∗_Na,500]. The model undergoes Hopf bifurcation at 𝑔_Na=𝑔^∗_Na.

Figure 3 shows the response of 𝑉 and 𝑚, ℎ, and 𝑛 to different 𝑔_Na. As the analysis above, when 𝑔_Na=198<𝑔^∗_Na, the system is stable. Figure 3(a) is the potential-time (𝑉-𝑡) curve, which shows that the action potential 𝑉 becomes steady. Figure 3(b) displays the trajectory of gating variables 𝑚, ℎ, and 𝑛 with time. We can see that the electrophysiological activity of cell reaches an equilibrium state at last.

(a) The 𝑉-𝑡 curve (𝑔Na=198<𝑔∗Na)

(b) The trajectory of 𝑚, ℎ, and 𝑛 (𝑔Na=198<𝑔∗Na)

(c) The 𝑉-𝑡 curve (𝑔Na=250>𝑔∗Na)

(d) The trajectory of 𝑚, ℎ, and 𝑛 (𝑔Na=250>𝑔∗Na)

Figure 3 

The response of 𝑉 and 𝑚, ℎ, and 𝑛 to different 𝑔Na.

Figures 3(c) and 3(d) demonstrate that the system is unstable when 𝑔_Na=250>𝑔^∗_Na. Figure 3(c) is the 𝑉-𝑡 graph, from which we can see the potential changes periodically. Figure 3(d) describes the trend of 𝑚, ℎ, and 𝑛, whose trajectory is a circle finally. Both Figures 3(c) and 3(d) imply that the system is unstable and the electrophysiological activity of cell is in a periodical state at a certain frequency.

3.2. Effect of 𝑔_K on the Model

In this part, we choose 𝑔_K as variable and keep the other parameters with ideal values. The same method with analysis of 𝑔_Na is taken to analyze the effect of 𝑔_K on the equilibrium, stability, and bifurcation of HH model. 𝑔_K∈[0,200] is taken into consideration because the desired value of 𝑔_K is 36.0. First, the relationship between 𝑔_K and the equilibrium 𝑉^∗ is obtained in Figure 4.

Figure 4 

The relationship between 𝑔K and 𝑉∗.

From Figure 4, we can see that 𝑉^∗ varies rapidly when 𝑔_K∈[0,20] and decreases slowly when 𝑔_K∈[30,200]. This means that equilibrium points are sensitive to 𝑔_K when 𝑔_K∈[0,20]. A slight change of 𝑔_K may make the final state of model change greatly.

Using the method of bisection to calculate the eigenvalues, we can find two bifurcation points 𝑔^∗_K1=3.843499029 and 𝑔^∗_K2=19.762260771 when 𝑔_K varies. Substitute 𝑔^∗_K𝑖 (𝑖=1,2) into (1), and obtain the corresponding equilibrium points 𝑉^∗. Both 𝑔^∗_K𝑖 and 𝑉^∗ are substituted into (7), and then corresponding eigenvalues could be obtained as follows: 𝜆¹₁𝜆=−5.3218099843,₁₂𝜆=−0.4223840650,₁₃=3.2×10⁻¹⁶𝜆−1.1305093754𝑖,₁₄=3.2×10⁻¹⁶𝜆+1.1305093754𝑖,₂₁𝜆=−4.5370272278,₂₂𝜆=−0.1319002182,₂₃=4.2×10⁻¹⁶𝜆−0.3436440068𝑖,₂₄=4.2×10⁻¹⁶+0.3436440068𝑖. (10)

Here, 3.2×10⁻¹⁶ and 4.2×10⁻¹⁶ can be approximately regarded as 0. From computing, all the real parts of eigenvalues are negative when 𝑔_K∈[0,𝑔^∗_K1)∪(𝑔^∗_K2,200], and all of them are not negative when 𝑔_K∈(𝑔^∗_K1,𝑔^∗_K2). According to the stability theory, the system is stable around equilibrium when 𝑔_K∈[0,𝑔^∗_K1)∪(𝑔^∗_K2,200] and it is unstable when 𝑔_K∈(𝑔^∗_K1,𝑔^∗_K2). The model undergoes Hopf bifurcations at 𝑔_K=𝑔^∗_K𝑖 (𝑖=1,2). The system is from locally stable state (𝑔_K∈[0,𝑔^∗_K1)) to unstable state (𝑔_K∈(𝑔^∗_K1,𝑔^∗_K2)) and becomes stable (𝑔_K∈(𝑔^∗_K2,200]) again. Responses of 𝑉 and 𝑚, ℎ, and 𝑛 to different 𝑔_K are shown in Figure 5.

(a) The 𝑉-𝑡 curve (𝑔K=2.8<𝑔∗K1)

(b) The trajectory of 𝑚, ℎ, and 𝑛 (𝑔K=2.8<𝑔∗K1)

(c) The 𝑉-𝑡 curve (𝑔∗K1<𝑔K=15<𝑔∗K2)

(d) The trajectory of 𝑚, ℎ, and 𝑛 (𝑔∗K1<𝑔K=15<𝑔∗K2)

(e) The 𝑉-𝑡 curve (𝑔K=21>𝑔∗K2)

(f) The trajectory of 𝑚, ℎ, and 𝑛 (𝑔K=21>𝑔∗K2)

Figure 5 

The response of 𝑉 and 𝑚, ℎ, and 𝑛 to different 𝑔K.

Figure 5 shows the response of 𝑉 and 𝑚, ℎ, and 𝑛 to different 𝑔_K. When 𝑔_K=2.8<𝑔^∗_K1, the system is stable. Figure 5(a) shows the trend of potential with time, from which we can see that the potential reaches a fixed value. Figure 5(b) is the trajectory of 𝑚, ℎ, and 𝑛 with time. All the gating variables also stay at fixed values (a steady point in Figure 5(b)) at last. These mean that the electrophysiological activity of cell reaches a steady state ultimately.

Figures 5(c) and 5(d) are 𝑉-𝑡 and 𝑚-ℎ-𝑛 graphs, respectively, when 𝑔^∗_K1<𝑔_K=15<𝑔^∗_K2. Figure 5(c) shows that the action potential changes in a certain period. And Figure 5(d) describes the trajectory of 𝑚, ℎ, and 𝑛 with time, from which we can find that the shape of the trajectory is a loop. Figures 5(c) and 5(d) imply all the gating variables and potential change periodically, which means that the electrophysiological activity of cell is in a periodical state.

Figures 5(e) and 5(f) are 𝑉-𝑡 and 𝑚-ℎ-𝑛 curves, respectively, when 𝑔_K=21>𝑔^∗_K2. From Figure 5(e) we can see that the potential reaches the resting state at this occasion. Figure 5(f) describes the trajectory of 𝑚, ℎ, and 𝑛, which shows that the three variables stay at a fixed point at last. Both Figures 5(e) and 5(f) show that all the potential and gating variables no longer change with time, which implies the cell reaches the resting state finally.

3.3. Effect of 𝑔_Na and 𝑔_K on the Model

Both 𝑔_Na and 𝑔_K are taken as variables in this part to study the stability and bifurcation of the model when 𝑔_Na∈[0,400] and 𝑔_K∈[0,60]. Keeping the other parameters with desired values, using Matlab as a tool, we get the equilibrium points first when 𝑔_Na and 𝑔_K both vary. Then the points are substituted into eigenmatrix of (4) and the eigenvalues of the model can be calculated. At last, Figure 6 is gained, in which all the real parts of eigenmatrix are negative if 𝑔_Na and 𝑔_K belong to the pink region and positive real parts appear if 𝑔_Na and 𝑔_K are in white area.

Figure 6 

The 𝑔Na-𝑔K plane and the critical boundary.

From Figure 6, we can find that the upper boundary of the regions is similar to a line. With the least square method applied, the expression of the line can be gotten as 𝑔_K=0.175×𝑔_Na−1.675. However, the lower boundary is not regular. According to stability theory, we can easily know that the model is stable when 𝑔_Na and 𝑔_K are in pink region and unstable when 𝑔_Na and 𝑔_K are in white. This means that the electrophysiological activity can reach a steady state when 𝑔_Na and 𝑔_K are in pink region and it is periodic when 𝑔_Na and 𝑔_K are in white. The system undergoes bifurcations when 𝑔_Na and 𝑔_K are on the boundary.

4. Discussion and Conclusion

The effects of 𝑔_Na and 𝑔_K on the stability and bifurcation of HH model are analyzed in the paper. The critical values and boundary are obtained. When 𝑔_Na increases to the critical value, the model will have bifurcation phenomenon, which means system will reach stable state when 𝑔_Na is less than the critical value and the cell will have continuous action potential after stimulation when 𝑔_Na is greater. However, there are two critical values about 𝑔_K. The system will be stable when 𝑔_K is less than the smaller critical value and there are periodic solutions when 𝑔_K is greater than the value and meanwhile is less than the larger one. The model will reach steady state again when 𝑔_K is greater than the larger critical value. From analyzing 𝑔_Na and 𝑔_K collectively, we can get a critical line which divides the 𝑔_Na-𝑔_K plane into two parts. The system will be stable when (𝑔_Na, 𝑔_K) is in the upper half plane and model will have periodic solutions when (𝑔_Na, 𝑔_K) is in the lower half.

In our analysis, when 𝑔_Na or 𝑔_K are taken as the variable(s), all the other parameters are kept with desired values. However, almost all the biological systems are coupled. All the components influence one another and work together forming the overall functionality. Therefore, when the sodium (𝑔_Na) and/or potassium (𝑔_K) channels vary, do the other parameters remain unchanged? Is it reasonable to keep the other parameters still with desired values? We could not ensure that it must be reasonable. Nevertheless, some evidences may explain a certain rationality of the method. For example, tetrodotoxin (TTX) selectively binds to the outer vestibule voltage-gated sodium channels, preventing channels from opening [31]. Ivabradine is a sinus node 𝐼_𝑓 channel inhibitor, which is selective for the 𝐼_𝑓 current but does not affect other cardiac ionic currents [32]. Acacetin could suppress the ultrarapid delayed rectifier K⁺ current and the transient outward K⁺ current and block the acetylcholine-activated K⁺ current; however, it has no effect on Na⁺ current, L-type Ca²⁺ current, or even inward-rectifier K⁺ current [33]. All of these demonstrate that to an extent when one channel changes, the others may not be affected. That is, when the parameter describing a channel varies, it is reasonable to keep parameters describing the other channels unchanged.

Stable states indicate that the electrophysiological activity of cell will get to corresponding resting state at last, while periodic phenomenon looks like response of pathological cell’s action potentials caused by cardiac arrhythmias [34]. In other words, 𝑔_Na and 𝑔_K may be the causes of the similar diseases to cardiac arrhythmias. So given appropriate medicine to change 𝑔_Na or 𝑔_K to reasonable intervals, the corresponding diseases could be abolished or the discomfort can be ameliorated. After all, our research could be a reference to treat relevant diseases. Some diseases led to by abnormal ion channels may be eased by medicine to adjust the conductance into corresponding intervals.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work is supported by the National Natural Science Foundation of China (NSFC) under Grant no. 61173086 and no. 61179009.