Integrated Approach in Systems BiologyView this Special Issue
Research Article | Open Access
Effects of Maximal Sodium and Potassium Conductance on the Stability of Hodgkin-Huxley Model
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 and maximal potassium conductance . Applying stability theory, and taking and 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 is the variable, while there are two points when is variable. The (, ) plane is partitioned into two regions and the upper bifurcation boundary is similar to a line when both and 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.
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 , , , , , , , , , , and [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 , 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  and sinoatrial node cell model . 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 , Butters’s atrial model , O’Hara’s ventricular model , and Li’s Purkinje cell model , 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 changed which was injected into the neuron from microelectrode . Stable and unstable solutions of the model with regard to were analyzed by Rinzel and Miller, and the influence of temperature was also discussed . Two stable steady states were found by Aihara and Matsumoto ; 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 and potassium ion potential on the bifurcations of the model . 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 and maximal potassium conductance [19, 20]. Fukai and his fellows examined the structure of the model’s bifurcations produced by and one of the other parameters . Taking leakage conductance and sodium channel effective bias voltage as parameters, Terada et al. analyzed bifurcation in Hodgkin-Huxley model for muscles of frogs .
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 .
In the past, for HH model, external current 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, and are seldom taken into consideration to analyze the stability of model, as the relevant data is not abundant. In this study, the effects of and on the stability and bifurcations of the model will be discussed, respectively, and collectively. And we will give the critical points of and when they play the role separately, and the critical boundaries in 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 and 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. , , and 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 , resulting in the rise of transmembrane potential. The open of potassium channels makes a large outflux of K+, creating the current 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 , , and . and vary with time.
(a) Electrophysiological activities of cell
(b) Electrical equivalent circuit of HH model
The equations were obtained according to electrical formulas and experimental data, which are shown as follows: where
In these equations, is the transmembrane potential. and are the gating variables indicating activation and inactivation of sodium ion current, respectively. is the gating variable showing activation of potassium ion current. , , and represent the maximal conductance of corresponding currents. μF/cm2 is membrane capacitance. is the current injected into the neuron. In our paper, we suppose and mS/cm2, mS/cm2, and mS/cm2, 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, Then the linearization of (1) around the equilibrium could be obtained as follows: where
We can get the eigenmatrix of (4): and then the characteristic equation can be obtained: where
According to Routh-Hurwitz criterion, if , , , , 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 on the Stability
In this section, we will investigate the influence of on the equilibrium, stability, and bifurcations of the model. is taken as variable, and the other parameters are all kept with desired values. Because the desired value of is around 120 mS/cm2, is taken into consideration. When 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 and the equilibrium point shown in Figure 2.
From Figure 2, we can see that changes slowly when and increases rapidly when . This means that equilibrium points are sensitive to when ; a slight change of 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 when changes.
Substituting into the original equation, we get the equilibrium , and then substituting both and into eigenmatrix of (4), we can gain the eigenvalues as follows:
Here, we regard as 0. With the help of computer, we can get that all the real parts of () are negative when . And there exist positive real parts when . According to the stability theory, the system is stable around equilibrium when and it is unstable when . The model undergoes Hopf bifurcation at .
Figure 3 shows the response of and , , and to different . As the analysis above, when , 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 ( )
(b) The trajectory of , , and ( )
(c) The curve ( )
(d) The trajectory of , , and ( )
Figures 3(c) and 3(d) demonstrate that the system is unstable when . 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 on the Model
In this part, we choose as variable and keep the other parameters with ideal values. The same method with analysis of is taken to analyze the effect of on the equilibrium, stability, and bifurcation of HH model. is taken into consideration because the desired value of is 36.0. First, the relationship between and the equilibrium is obtained in Figure 4.
From Figure 4, we can see that varies rapidly when and decreases slowly when . This means that equilibrium points are sensitive to when . A slight change of may make the final state of model change greatly.
Using the method of bisection to calculate the eigenvalues, we can find two bifurcation points and when varies. Substitute () into (1), and obtain the corresponding equilibrium points . Both and are substituted into (7), and then corresponding eigenvalues could be obtained as follows:
Here, and can be approximately regarded as 0. From computing, all the real parts of eigenvalues are negative when , and all of them are not negative when . According to the stability theory, the system is stable around equilibrium when and it is unstable when . The model undergoes Hopf bifurcations at (). The system is from locally stable state () to unstable state () and becomes stable () again. Responses of and , , and to different are shown in Figure 5.
(a) The curve
(b) The trajectory of , , and ( )
(c) The curve ( )
(d) The trajectory of , , and ( )
(e) The curve ( )
(f) The trajectory of , , and ( )
Figure 5 shows the response of and , , and to different . When , 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 . 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 . 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 and on the Model
Both and are taken as variables in this part to study the stability and bifurcation of the model when and . Keeping the other parameters with desired values, using Matlab as a tool, we get the equilibrium points first when and 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 and belong to the pink region and positive real parts appear if and are in white area.
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 . However, the lower boundary is not regular. According to stability theory, we can easily know that the model is stable when and are in pink region and unstable when and are in white. This means that the electrophysiological activity can reach a steady state when and are in pink region and it is periodic when and are in white. The system undergoes bifurcations when and are on the boundary.
4. Discussion and Conclusion
The effects of and on the stability and bifurcation of HH model are analyzed in the paper. The critical values and boundary are obtained. When increases to the critical value, the model will have bifurcation phenomenon, which means system will reach stable state when is less than the critical value and the cell will have continuous action potential after stimulation when is greater. However, there are two critical values about . The system will be stable when is less than the smaller critical value and there are periodic solutions when is greater than the value and meanwhile is less than the larger one. The model will reach steady state again when is greater than the larger critical value. From analyzing and collectively, we can get a critical line which divides the plane into two parts. The system will be stable when , is in the upper half plane and model will have periodic solutions when , is in the lower half.
In our analysis, when or 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 () and/or potassium () 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 . Ivabradine is a sinus node channel inhibitor, which is selective for the current but does not affect other cardiac ionic currents . 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 Ca2+ current, or even inward-rectifier K+ current . 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 . In other words, and may be the causes of the similar diseases to cardiac arrhythmias. So given appropriate medicine to change or 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.
This work is supported by the National Natural Science Foundation of China (NSFC) under Grant no. 61173086 and no. 61179009.
- A. L. Hodgkin and A. F. Huxley, “The components of membrane conductance in the giant axon of Loligo,” The Journal of Physiology, vol. 116, no. 4, pp. 473–496, 1952.
- A. L. Hodgkin and A. F. Huxley, “The components of membrane conductance in the giant axon of Loligo,” The Journal of Physiology, vol. 116, no. 4, pp. 473–496, 1952.
- A. L. Hodgkin and A. F. Huxley, “The dual effect of membrane potential on sodium conductance in the giant axon of Loligo,” The Journal of Physiology, vol. 116, no. 4, pp. 497–506, 1952.
- A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve.,” The Journal of Physiology, vol. 117, no. 4, pp. 500–544, 1952.
- C. Soeller, I. D. Jayasinghe, P. Li, A. V. Holden, and M. B. Cannell, “Three-dimensional high-resolution imaging of cardiac proteins to construct models of intracellular Ca2+ signalling in rat ventricular myocytes,” Experimental Physiology, vol. 94, no. 5, pp. 496–508, 2009.
- M. Wussling and G. Szymanski, “Simulation by two calcium store models of myocardial dynamic properties: potentiation, staircase, and biphasic tension development,” General Physiology and Biophysics, vol. 5, no. 2, pp. 135–152, 1986.
- A. Mahajan, Y. Shiferaw, D. Sato et al., “A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates,” Biophysical Journal, vol. 94, no. 2, pp. 392–410, 2008.
- D. Noble, “Cardiac action and pacemaker potentials based on the Hodgkin-Huxley equations,” Nature, vol. 188, no. 4749, pp. 495–497, 1960.
- J. Kneller, R. J. Ramirez, D. Chartier, M. Courtemanche, and S. Nattel, “Time-dependent transients in an ionically based mathematical model of the canine atrial action potential,” American Journal of Physiology—Heart and Circulatory Physiology, vol. 282, no. 4, pp. H1437–H1451, 2002.
- H. Zhang, A. V. Holden, I. Kodama et al., “Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node,” American Journal of Physiology: Heart and Circulatory Physiology, vol. 279, no. 1, pp. H397–H421, 2000.
- A. O. Verkerk and R. Wilders, “Hyperpolarization-activated current, If, in mathematical models of rabbit sinoatrial node pacemaker cells,” BioMed Research International, vol. 2013, Article ID 872454, 18 pages, 2013.
- T. D. Butters, O. V. Aslanidi, J. Zhao, B. Smaill, and H. Zhang, “A novel computational sheep atria model for the study of atrial fibrillation,” Interface Focus, vol. 3, no. 2, Article ID 20120067, 2013.
- T. O'Hara, L. Virág, A. Varró, and Y. Rudy, “Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation,” PLoS Computational Biology, vol. 7, no. 5, Article ID e1002061, 2011.
- P. Li and Y. Rudy, “A model of canine purkinje cell electrophysiology and Ca2+ cycling: rate dependence, triggered activity, and comparison to ventricular myocytes,” Circulation Research, vol. 109, no. 1, pp. 71–79, 2011.
- B. Hassard, “Bifurcation of periodic solutions of the Hodgkin-Huxley model for the squid giant axon,” Journal of Theoretical Biology, vol. 71, no. 3, pp. 401–420, 1978.
- J. Rinzel and R. N. Miller, “Numerical calculation of stable and unstable periodic solutions to the Hodgkin-Huxley equations,” Mathematical Biosciences, vol. 49, no. 1, pp. 27–59, 1980.
- K. Aihara and G. Matsumoto, “Two stable steady states in the Hodgkin-Huxley axons,” Biophysical Journal, vol. 41, no. 1, pp. 87–89, 1983.
- J. Guckenheimer and J. S. Labouriau, “Bifurcation of the Hodgkin and Huxley equations: a new twist,” Bulletin of Mathematical Biology, vol. 55, no. 5, pp. 937–952, 1993.
- Y. A. Bedrov, G. N. Akoev, and O. E. Dick, “On the relationship between the number of negative slope regions in the voltage-current curve of the Hodgkin-Huxley model and its parameter values,” Biological Cybernetics, vol. 73, no. 2, pp. 149–154, 1995.
- Y. A. Bedrov, G. N. Akoev, and O. E. Dick, “Partition of the Hodgkin-Huxley type model parameter space into the regions of qualitatively different solutions,” Biological Cybernetics, vol. 66, no. 5, pp. 413–418, 1992.
- H. Fukai, S. Doi, T. Nomura, and S. Sato, “Hopf bifurcations in multiple-parameter space of the Hodgkin-Huxley equations I. Global organization of bistable periodic solutions,” Biological Cybernetics, vol. 82, no. 3, pp. 215–222, 2000.
- K. Terada, H. A. Tanaka, and S. Yoshizawa, “Two-parameter bifurcations in the Hodgkin–Huxley equations for muscle fibers,” Electronics and Communications in Japan, vol. 83, no. 6, pp. 86–94, 2000.
- W. Jiang, C. Yanqiu, F. Xiangyang, and L. Li, “Multi-parameter Hopf-bifurcation in Hodgkin-Huxley model exposed to ELF external electric field,” Chaos, Solitons & Fractals, vol. 26, no. 4, pp. 1221–1229, 2005.
- J. Wang, L. Chen, and X. Fei, “Analysis and control of the bifurcation of Hodgkin-Huxley model,” Chaos, Solitons and Fractals, vol. 31, no. 1, pp. 247–256, 2007.
- J. Wang, L. Q. Chen, and X. Y. Fei, “Bifurcation control of the Hodgkin-Huxley equations,” Chaos, Solitons & Fractals, vol. 33, no. 1, pp. 217–224, 2007.
- Y. Che, J. Wang, B. Deng, X. Wei, and C. Han, “Bifurcations in the Hodgkin-Huxley model exposed to DC electric fields,” Neurocomputing, vol. 81, pp. 41–48, 2012.
- P. M. Hao, D. J. Fan, J. J. Wei, and Q. Liu, “Dynamic behaviors of a delayed HIV model with stage-structure,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 12, pp. 4753–4766, 2012.
- X. Y. Chang and J. J. Wei, “Hopf bifurcation and optimal control in a diffusive predator-prey system with time delay and prey harvesting,” Nonlinear Analysis: Modelling and Control, vol. 17, no. 4, pp. 379–409, 2012.
- X. Y. Chang and J. J. Wei, “Stability and Hopf bifurcation in a diffusive predator-prey system incorporating a prey refuge,” Mathematical Biosciences and Engineering, vol. 10, no. 4, pp. 979–996, 2013.
- D. M. Durand and M. Bikson, “Suppression and control of eplleptiform activity by electrical stimulation: a review,” Proceedings of the IEEE, vol. 89, no. 7, pp. 1065–1082, 2001.
- B. Venkatesh, S. Q. Lu, N. Dandona, S. L. See, S. Brenner, and T. W. Soong, “Genetic basis of tetrodotoxin resistance in pufferfishes,” Current Biology, vol. 15, no. 22, pp. 2069–2072, 2005.
- S. Sulfi and A. D. Timmis, “Ivabradine—the first selective sinus node If channel inhibitor in the treatment of stable angina,” International Journal of Clinical Practice, vol. 60, no. 2, pp. 222–228, 2006.
- G. Li, H. Wang, G. Qin et al., “Acacetin, a natural flavone, selectively inhibits human atrial repolarization potassium currents and prevents atrial fibrillation in dogs,” Circulation, vol. 117, no. 19, pp. 2449–2457, 2008.
- H. O. Wang, D. Chen, and L. G. Bushnell, “Control of bifurcations and chaos in heart rhythms,” in Proceedings of the 36th IEEE Conference on Decision and Control, vol. 1, pp. 395–400, San Diego, Calif, USA, December 1997.
Copyright © 2014 Yue Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.