- Email: [email protected]
A note on the asymptotic reduction of the Hodgkin–Huxley equations for nerve impulses
Bulletin of Mathematical Biology 67 (2005) 947–955 www.elsevier.com/locate/ybulm
A note on the asymptotic reduction of the Hodgkin–Huxley equations for nerve impulses Robert Hinch∗ Centre for Mathematical Biology, Mathematical Institute, 24-29 St Giles’, Oxford, OX1 3LB, UK University Laboratory of Physiology, Parks Road, Oxford, OX1 3PT, UK Received 28 April 2004; accepted 11 November 2004
Abstract The (standard) FitzHugh reduction of the Hodgkin–Huxley equations for the propagation of nerve impulses ignores the dynamics of the activation gates. This assumption is invalid and leads to an overestimation of the wave speed by a factor of 5 and the wrong dependence of wave speed on sodium channel conductance. The error occurs because a non-dimensional parameter, which is assumed to be small in the FitzHugh reduction, is in fact large (≈18). We analyse the Hodgkin–Huxley equations for propagating nerve impulses in the limit that this non-dimensional parameter is large, and show that the analytical results are consistent with numerical simulations of the Hodgkin–Huxley equations. © 2005 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
1. Introduction Neurons propagate electrical signals, called action potentials, by changing the conductance of their cell membrane to ions such as sodium and potassium. A standard model has been developed which combines models of the local membrane conductivity to different ions and the response to spatial gradients in the membrane potential (Keener and Sneyd, 1998). The Hodgkin–Huxley equations (Hodgkin and Huxley, 1952) are a model of the membrane potential of axons, and describe how voltage-gated ion ∗ Corresponding address: Centre for Mathematical Biology, Mathematical Institute, 24-29 St Giles’, Oxford, OX1 3LB, UK. Tel.: +44 1865 280614; fax: +44 1865 270515. E-mail address: [email protected].
0092-8240/$30 © 2005 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.11.007
948
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
channels dynamically alter their conductances in response to changes in membrane potential. The spatial variation of the membrane potential along the length of a neuron is modelled using the cable equation. The combined model equations are (Hodgkin and Huxley, 1952) ∂v 2 ∂v = σ 2 − gNa m N h(v − vNa ) − gK n 4 (v − vK ) − gL (v − vL ), ∂t ∂x m ∞,N (v) − m ∂m = , ∂t τm (v) ∂h h ∞ (v) − h = , ∂t τh (v) ∂n n ∞ (v) − n = , ∂t τn (v)
Cm
(1)
where v is the membrane potential; m, h and n are the gating variables; m ∞,N , h ∞ and n ∞ are the steady-state gate variables; τm , τh and τn are the gate time constants; Cm is the membrane capacitance; vNa , vK and vL are the reversal potentials of the sodium, potassium and leak currents respectively; gNa , gK and gL are the conductances of the sodium, potassium and leak currents respectively; σ is the conductance of the cable; and N is the number of activation gates. The potential v has been shifted so that the resting potential is v = 0 and the functional forms for m ∞,N , h ∞ , n ∞ , τm , τh and τn are given in Hodgkin and Huxley (1952) (reprinted in Keener and Sneyd (1998)). While it is simple to accurately measure the current–membrane potential (I –V ) relation for the sodium current, it harder to estimate N and τm accurately. Hodgkin and Huxley estimated N = 3 by fitting the model to voltage-clamp data. These equations (1) have a travelling wave solution with constant profile (Fig. 1A), which represents a propagating action potential. The rapid upstroke of the action potential is due to a large inward sodium current. The sodium channel is controlled by both activation (m) and inactivation (h) gates. Since the upstroke of the action potential is narrow compared with the width of the pulse, the speed of a single pulse is determined by the solution in the region of the upstroke, and the recovery process only affects the velocity by a small amount (see (22) of Rinzel and Keller (1973)). Note that for a periodic wave train of pulses, it is necessary to consider the recovery process, even though the upstroke is narrow (Rinzel and Keller, 1973). While it is easy to solve (1) numerically, many attempts have been made to derive simplified models which can be analysed mathematically. The most famous of these is the FitzHugh reduction (FitzHugh, 1960; Nagumo et al., 1962) on the spatially uniform version of (1). The FitzHugh reduction is based on a separation of time-scales. First it is assumed that the n and h gates vary slowly, so are constant in the wavefront. The second assumptions is that the time constant of the activation gates (τm ) is the fastest time-scale in the model, so that m ≈ m ∞,N (v) and the sodium current is N INa ≈ I¯Na = gNa m ∞,N (v)h.
(2)
A more formal derivation of the FitzHugh reduction is given by Krinsky and Kokov (1973). The FitzHugh reduction has been used to reduce the Hodgkin–Huxley model for nerve conduction (1) to a third-order system, where it was assumed that the time constant of the
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
949
Fig. 1. A. The membrane potential across the wavefront of a propagated action potential in the Hodgkin–Huxley model. Note that the width of the wavefront is approximately 0.3 ms and is much smaller than the width of the of the pulse. B. The activation gate m in the wavefront, the steady-state value m ∞,3 (v) and their difference. In the FitzHugh reduction it is assumed that m takes its steady-state value.
activation gates (τm ) is fast and (2) is used to approximate the sodium channel current (Rinzel and Keller, 1973; Casten et al., 1975). However, Fig. 1B demonstrates that this is not a valid assumption in the wavefront of action potentials, where the time-scale of the depolarisation is the same as that as the activation gates. As mentioned earlier, it is simple to measure the sodium current’s I –V relationship, but much harder to measure the value of N accurately. If we choose m ∞,N (V ) such that the I –V curve is independent of N, then we can measure the effect of N on the propagation velocity. This is easily achieved by setting m ∞,N = (m ∞,3 )3/N , where m ∞,3 is the steady-state function for the m-gates given by Hodgkin and Huxley (1952). Under the FitzHugh approximation with this definition of m ∞,N , the wave speed is independent of N because INa will be independent of N (see (2)). Fig. 2 shows the speed of the wave as a function of gNa for different values of N. Note that the FitzHugh approximation over-estimates the wave speed by a factor of 5 compared with the standard Hodgkin–Huxley parameters (N = 3 and gNa = 120 mS cm−2 ). Additionally, the relationship between wave speed and sodium channel conductance is much weaker the higher the number of activation gates. The FitzHugh reduction is asymptotic in the limit the non-dimensional parameter g = gNa τm (0)h ∞ (0)/Cm 1; however, using physiological parameters g ≈ 17.9. The physiological interpretation of g is the ratio of the time-scale of activation of the m-gate to the time-scale of depolarisation when all sodium channels are activated. The effect of activation gates on the propagation velocity in simplified polynomial (Hunter et al., 1975) and threshold (Hinch, 2002) models has been analysed previously; however, the physiologically realistic Hodgkin–Huxley equations have not been analysed. An alternative approach to the FitzHugh reduction is the method of equivalent potentials, where the gating variables are replaced by equivalent potential variables (e.g. vh (t)), which are defined implicitly via h(t) = h ∞ (vh (t)) (Kepler et al., 1992). The resulting equations can then be solved using a regular perturbation series; however, this technique has not been extended to the problem of propagating action potentials. It should be noted the ordering of time-scales for nerve cells is different from that of cardiac cells.
950
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
Fig. 2. The speed of the action potential as a function of the sodium channel conductance for different numbers of activation gates. The steady-state function for the m-gates was chosen such that the I –V relationship of the sodium current was independent of N (i.e. m ∞,N = (m ∞,3 )3/N ).
In cardiac cells, the time-scale of depolarisation is much slower due to the lower density of sodium channels. The sodium channel inactivation gates then become important, which leads to multiple wavefront solutions with different speeds (Hinch, 2002; Biktashev, 2002) of which only one is stable (Hinch, 2004). We calculate the travelling wave solution of the Hodgkin–Huxley model in the limit g 1 and show that the wave speed is proportional to g 1/(2N+2) . 2. Transition layer in travelling wave solutions We will demonstrate that in the limit g → ∞, a transition layer forms in the Hodgkin–Huxley equations. The approximation is based on a separation of time-scales, which are the same as those in the FitzHugh reduction with the exception that we do not assume that τm is rapid. The first stage of the analysis is to non-dimensionalise the equations. The non-dimensional functions Ti (v) are defined by Ti (v) = τi (v)/τi (0), where i = m, h, n. Define the non-dimensional variables τm (0)Cm t , xˆ = x , (3) tˆ = τm (0) σ the shifted m gate variable m¯ = m − m ∞,N (0),
m¯ ∞ (v) = m ∞,N (v) − m ∞,N (0),
(4)
and the non-dimensional parameters τm (0) τm (0) ,
n = , τh (0) τn (0) (5) gKτm (0)n 4∞ (0) gL τm (0) gNa τm (0)h ∞ (0)
K = ,
L = , g= . Cm Cm Cm Using the Hodgkin–Huxley parameters (Hodgkin and Huxley, 1952), we find that
m , h , n , K , L 1, g 1 and g mN K , L . The physiological interpretations of
m = m ∞,N (0),
h =
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
951
these small parameters are: the time constants of the h- and n-gates are slow compared with that of the m-gate ( h , n 1); at the resting potential the m-gate is almost completely closed ( m 1); the time-scale of the rate of change of the membrane potential due to the potassium and leakage currents at the resting potential is small compared with the time constant of the m-gate ( K , L 1); if the sodium channels were fully open the rate of change of the membrane potential would be quicker than the time constant of the m-gate (g 1); and the resting potential is maintained by balancing the potassium current with the leak current, and the sodium current is not important in determining the resting potential (g mN K , L ). Therefore, to leading order, the resting potential is independent of g. The non-dimensional equations are ∂v h n4 ∂v 2 (v − vNa ) − K 4 ¯ N (v − vK ) − L (v − vL ), = 2 − g( m + m) h ∞ (0) ∂ xˆ n ∞ (0) ∂ tˆ m¯ ∞ (v) − m¯ ∂ m¯ , = Tm (v) ∂ tˆ (6) ∂h h ∞ (v) − h , = h Th (v) ∂ tˆ ∂n n ∞ (v) − n . = n ˆ Tn (v) ∂t Next, we look for travelling wave solutions with constant profile by introducing the travelling wave co-ordinate ξ = tˆ − x/c, ˆ where c is the wave speed. We consider the asymptotic limit g → ∞ and m g → 0. We introduce the re-scaled variables ξ , c and m¯ c = Cg β ,
ξ = Xg −α ,
m¯ = Mg −γ
(7)
where α, β, γ > 0 and are determined later. The physiological meanings of these rescalings are: the wave speed is large, the depolarisation takes place in a narrow transition layer and the depolarisation occurs when the m-gates are only partially open. Inserting (7) into (6) yields dv 1 dv 2 h = g α−2β 2 (v − vNa ) − g 1−α−γ N M N dX C dX 2 h ∞ (0) + O( m g 1−α−γ (N−1) , K g −α , L g −α ), dM m¯ ∞ (v) = g γ −α + O(g −α ), dX Tm (v) dh = O( h g −α ), dX dn = O( n g −α ). dX
(8)
The final two equations tell us that h and n are constant throughout the wavefront, which is one of FitzHugh’s approximations. The leading order balance for the first two equations occurs when α = γ = 1/(1 + N) and β = 1/(2 + 2N). The first two equations in the
952
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
transition layer are 1 d2 v dv = 2 − M N (v − vNa ) + O( m g α , K g −α , L g −α ), dX C dX 2 (9) dM −α = f (v) + O(g ), dX where f (v) = m¯ ∞ (v)/Tm (v) and has the properties f (0) = 0 and f > 0 (note that to leading order v = 0 is the resting potential because g mN K , L ). The equation for v has fixed-points when M = 0 (which occurs when the tissue is at rest so v = 0) and when v = vNa . The travelling wave solution is the hetero-clinic orbit connecting these two points, so the boundary conditions for a travelling wave solution are M(−∞) = v(−∞) = 0 and v(∞) = vNa . It should be noted that at large times when v → vNa , M will grow linearly with time. At first this may seem a problem because m < 1; however, since m = M/g α then m is asymptotically small even though M is large. For long times (X = O(g α )), it is necessary to return to the un-scaled variables and to linearise around the fixed points. This solution can then be linked to the transition layer solution using matched asymptotic expansions; however, it is omitted since it adds nothing of interest because the physiological condition that v is bounded determines the boundary condition v → vNa as X → ∞ within the transition layer. The transition layer problem is an eigenvalue problem for C; however, the leading-order terms are independent of g, therefore c ∝ g 1/(2N+2) ,
(10)
which agrees with a previous analysis of a threshold model (Hinch, 2002). The simple form of (9) suggests that it may be possible to prove existence, uniqueness and stability properties of the transition layer solution. Note that one of the correction terms is of order
m g α (≈0.1 using the Hodgkin–Huxley parameters), which grows at large g. Therefore, to check that the transition layer solution (7) and (9) is asymptotic as g → ∞, it is necessary to ensure that m g α 1. This can be achieved by a minor modification of the I –V relationship of the sodium current (see Fig. 3B). The physiological justification for this modification is that it is known that the resting potential is maintained by the balancing of the potassium and leak currents and not the sodium current. When we take the asymptotic limit g → ∞, the sodium current would become important at the resting potential which would be physiologically incorrect. The modification therefore allows us to examine the asymptotic limit g → ∞ without changing the underlying physiological mechanisms. Fig. 3A shows a comparison of the wave speed of the numerical solution of (1) with the transition layer solution ((7) and (9)) on a log–log plot for different values of N. Details of the numerical schemes used are given in the Appendix. The straight lines in Fig. 3A are the transition layer solutions and the constant C was calculated numerically. Note that, at large g, the numerical solutions asymptote to the transition layer solution. Fig. 4 shows the profile of the wavefront for different values of g compared with the transition layer solution, when N = 1 (A) and N = 3 (B). The time variable is rescaled by a factor of g 1/(N+1) , so that a direct comparison of different values of g can be made. As g → ∞, the N = 1 solution asymptotes quicker to the transition layer solution than the N = 3 solution. This is because the main correction term in (9) is proportional to g −1/(N+1) , which decays
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
953
Fig. 3. A. The wave speed as a function of sodium channel conductance for different values of N . The straight lines are the transition layer solutions and the curves are the full solution of (1). Note that the full solutions asymptote to the transition layer solutions at large g. B. The modification to the I –V relation of the sodium current so that it is possible to show the transition layer solutions are asymptotic.
Fig. 4. The profile of the wavefronts when N = 1 (A) and N = 3 (B). In both graphs, the lines from left to right are: the transition layer solution calculated from (9); then gNa = 10 000 mS cm−2 (g = 1483), gNa = 1000 mS cm−2 (g = 148) and gNa = 120 mS cm−2 (g = 17.8) calculated numerically from (1). Time was rescaled by a factor of g 1/(N +1) in all cases so a direct comparison between different values of g could be made (i.e. in the limit of large g the profile will tend to the same shape). The size of the error in the asymptotic calculation is g −1/(N +1) , which takes the values 0.026, 0.082 and 0.237 when N = 1 (A), and 0.161, 0.287 and 0.487 when N = 3 (B). In both cases, the profile of the solution full equations asymptotes the transition layer solution as g is increased.
more slowly with g when N is increased. Using the Hodgkin–Huxley parameter values, we find that g −1/(N+1) ≈ 17.9−1/4 ≈ 0.49. Unfortunately the power of one quarter dramatically increases the size of the small parameter thus reducing the accuracy of the transition layer solution (Fig. 4B). In this note we have demonstrated that the FitzHugh reduction of the Hodgkin–Huxley equations for propagating nerve impulses is not valid when physiological parameters are used, and this led to an over-estimation of the wave speed by a factor of 5 (Fig. 2). An alternative asymptotic analysis was carried out which included the dynamics of the
954
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
activation gates. The results of this alternative analysis agree with numerical solutions of the Hodgkin–Huxley equations, and demonstrates the importance of including the dynamics of the activation gates when calculating the speed of propagation of nerve impulses. Acknowledgements This work was supported by The Wellcome Trust. RH would like to thank Santiago Schnell for useful comments on this manuscript. Appendix The Hodgkin–Huxley equations were solved using a standard shooting method (FitzHugh and Antosiewicz, 1959). The travelling wave co-ordinate ξ = t − x/c was introduced to (1) which produces a set of five coupled first-order ordinary differential equations (ODEs), creating an eigenvalue problem for c. The resting potential (i.e. fixedpoint) was found using the Newton–Raphson method and the Jacobian matrix was shown to have a single unstable eigenvalue. An initial guess for the eigenvalue c was made. The membrane potential was increased by 0.01 mV above the resting potential and the ODEs were integrated using the fourth-order Runge–Kutta method. The potential either diverged to +∞ implying that the trial value for c was too large, or diverged to −∞ implying that the trial value for c was too small. The eigenvalue was then found iteratively using interval bisection. This method enables the efficient calculation of the wave velocity; however, it is unable to resolve the full pulse profile (FitzHugh and Antosiewicz, 1959; Miller and Rinzel, 1981), which is a standard problem when using shooting in problems with narrow transition layers. Selective points on the graphs were checked by directly solving the partial differential equations using an explicit Euler finite difference scheme, and both methods agreed with each other. However, direct computation of the PDEs is both slow and inaccurate compared with the shooting method. References Biktashev, V.N., 2002. Dissipation of the excitation wavefronts. Phys. Rev. Lett. 89, 168102. Casten, R.G., Cohen, H., Lagerstom, P.A., 1975. Perturbation analysis of an approximation to the Hodgkin–Huxley theory. Quart. Appl. Math. 4, 365–402. FitzHugh, R.A., 1960. Thresholds and plateaus in the Hodgkin–Huxley nerve equations. J. Gen. Physiol. 43, 867–896. FitzHugh, R.A., Antosiewicz, H.A., 1959. Automated computation of nerve excitation-detailed corrections and additions. J. SIAM 7, 447–457. Hinch, R., 2002. An analytical study of the physiology and pathology of the propagation of cardiac action potentials. Prog. Biophys. Mol. Biol. 78, 45–81. Hinch, R., 2004. Stability of cardiac waves. Bull. Math. Biol. 66, 1887–1908. Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerves. J. Physiol. 117, 507–544. Hunter, P.J., McNaughton, P.A., Noble, D., 1975. Analytical models of propagation in excitable cells. Prog. Biophys. Mol. Biol. 30, 99–144. Keener, J., Sneyd, J., 1998. Mathematical Physiology. Springer.
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
955
Kepler, T.B., Abbott, L.F., Marder, E., 1992. Reduction of conductance-based neuron model. Biol. Cybern. 66, 381–387. Krinsky, V.I., Kokov, Y.M., 1973. Analysis of the equations of excitable membranes reducing the Hodgkin–Huxley equations to a second order system. Biofizika 18, 506–511. Miller, R.N., Rinzel, J., 1981. The dependence of impulse propagation speed on firing, frequency, dispersion, for the Hodgkin–Huxley model. Biophys. J. 7, 227–259. Nagumo, J., Arimoto, S., Yoshizawa, S., 1962. An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061. Rinzel, J., Keller, J.B., 1973. Traveling wave solutions of a nerve conduction equation. Biophys. J. 13, 1313–1337.
A note on the asymptotic reduction of the Hodgkin–Huxley equations for nerve impulses Robert Hinch∗ Centre for Mathematical Biology, Mathematical Institute, 24-29 St Giles’, Oxford, OX1 3LB, UK University Laboratory of Physiology, Parks Road, Oxford, OX1 3PT, UK Received 28 April 2004; accepted 11 November 2004
Abstract The (standard) FitzHugh reduction of the Hodgkin–Huxley equations for the propagation of nerve impulses ignores the dynamics of the activation gates. This assumption is invalid and leads to an overestimation of the wave speed by a factor of 5 and the wrong dependence of wave speed on sodium channel conductance. The error occurs because a non-dimensional parameter, which is assumed to be small in the FitzHugh reduction, is in fact large (≈18). We analyse the Hodgkin–Huxley equations for propagating nerve impulses in the limit that this non-dimensional parameter is large, and show that the analytical results are consistent with numerical simulations of the Hodgkin–Huxley equations. © 2005 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
1. Introduction Neurons propagate electrical signals, called action potentials, by changing the conductance of their cell membrane to ions such as sodium and potassium. A standard model has been developed which combines models of the local membrane conductivity to different ions and the response to spatial gradients in the membrane potential (Keener and Sneyd, 1998). The Hodgkin–Huxley equations (Hodgkin and Huxley, 1952) are a model of the membrane potential of axons, and describe how voltage-gated ion ∗ Corresponding address: Centre for Mathematical Biology, Mathematical Institute, 24-29 St Giles’, Oxford, OX1 3LB, UK. Tel.: +44 1865 280614; fax: +44 1865 270515. E-mail address: [email protected].
0092-8240/$30 © 2005 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.11.007
948
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
channels dynamically alter their conductances in response to changes in membrane potential. The spatial variation of the membrane potential along the length of a neuron is modelled using the cable equation. The combined model equations are (Hodgkin and Huxley, 1952) ∂v 2 ∂v = σ 2 − gNa m N h(v − vNa ) − gK n 4 (v − vK ) − gL (v − vL ), ∂t ∂x m ∞,N (v) − m ∂m = , ∂t τm (v) ∂h h ∞ (v) − h = , ∂t τh (v) ∂n n ∞ (v) − n = , ∂t τn (v)
Cm
(1)
where v is the membrane potential; m, h and n are the gating variables; m ∞,N , h ∞ and n ∞ are the steady-state gate variables; τm , τh and τn are the gate time constants; Cm is the membrane capacitance; vNa , vK and vL are the reversal potentials of the sodium, potassium and leak currents respectively; gNa , gK and gL are the conductances of the sodium, potassium and leak currents respectively; σ is the conductance of the cable; and N is the number of activation gates. The potential v has been shifted so that the resting potential is v = 0 and the functional forms for m ∞,N , h ∞ , n ∞ , τm , τh and τn are given in Hodgkin and Huxley (1952) (reprinted in Keener and Sneyd (1998)). While it is simple to accurately measure the current–membrane potential (I –V ) relation for the sodium current, it harder to estimate N and τm accurately. Hodgkin and Huxley estimated N = 3 by fitting the model to voltage-clamp data. These equations (1) have a travelling wave solution with constant profile (Fig. 1A), which represents a propagating action potential. The rapid upstroke of the action potential is due to a large inward sodium current. The sodium channel is controlled by both activation (m) and inactivation (h) gates. Since the upstroke of the action potential is narrow compared with the width of the pulse, the speed of a single pulse is determined by the solution in the region of the upstroke, and the recovery process only affects the velocity by a small amount (see (22) of Rinzel and Keller (1973)). Note that for a periodic wave train of pulses, it is necessary to consider the recovery process, even though the upstroke is narrow (Rinzel and Keller, 1973). While it is easy to solve (1) numerically, many attempts have been made to derive simplified models which can be analysed mathematically. The most famous of these is the FitzHugh reduction (FitzHugh, 1960; Nagumo et al., 1962) on the spatially uniform version of (1). The FitzHugh reduction is based on a separation of time-scales. First it is assumed that the n and h gates vary slowly, so are constant in the wavefront. The second assumptions is that the time constant of the activation gates (τm ) is the fastest time-scale in the model, so that m ≈ m ∞,N (v) and the sodium current is N INa ≈ I¯Na = gNa m ∞,N (v)h.
(2)
A more formal derivation of the FitzHugh reduction is given by Krinsky and Kokov (1973). The FitzHugh reduction has been used to reduce the Hodgkin–Huxley model for nerve conduction (1) to a third-order system, where it was assumed that the time constant of the
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
949
Fig. 1. A. The membrane potential across the wavefront of a propagated action potential in the Hodgkin–Huxley model. Note that the width of the wavefront is approximately 0.3 ms and is much smaller than the width of the of the pulse. B. The activation gate m in the wavefront, the steady-state value m ∞,3 (v) and their difference. In the FitzHugh reduction it is assumed that m takes its steady-state value.
activation gates (τm ) is fast and (2) is used to approximate the sodium channel current (Rinzel and Keller, 1973; Casten et al., 1975). However, Fig. 1B demonstrates that this is not a valid assumption in the wavefront of action potentials, where the time-scale of the depolarisation is the same as that as the activation gates. As mentioned earlier, it is simple to measure the sodium current’s I –V relationship, but much harder to measure the value of N accurately. If we choose m ∞,N (V ) such that the I –V curve is independent of N, then we can measure the effect of N on the propagation velocity. This is easily achieved by setting m ∞,N = (m ∞,3 )3/N , where m ∞,3 is the steady-state function for the m-gates given by Hodgkin and Huxley (1952). Under the FitzHugh approximation with this definition of m ∞,N , the wave speed is independent of N because INa will be independent of N (see (2)). Fig. 2 shows the speed of the wave as a function of gNa for different values of N. Note that the FitzHugh approximation over-estimates the wave speed by a factor of 5 compared with the standard Hodgkin–Huxley parameters (N = 3 and gNa = 120 mS cm−2 ). Additionally, the relationship between wave speed and sodium channel conductance is much weaker the higher the number of activation gates. The FitzHugh reduction is asymptotic in the limit the non-dimensional parameter g = gNa τm (0)h ∞ (0)/Cm 1; however, using physiological parameters g ≈ 17.9. The physiological interpretation of g is the ratio of the time-scale of activation of the m-gate to the time-scale of depolarisation when all sodium channels are activated. The effect of activation gates on the propagation velocity in simplified polynomial (Hunter et al., 1975) and threshold (Hinch, 2002) models has been analysed previously; however, the physiologically realistic Hodgkin–Huxley equations have not been analysed. An alternative approach to the FitzHugh reduction is the method of equivalent potentials, where the gating variables are replaced by equivalent potential variables (e.g. vh (t)), which are defined implicitly via h(t) = h ∞ (vh (t)) (Kepler et al., 1992). The resulting equations can then be solved using a regular perturbation series; however, this technique has not been extended to the problem of propagating action potentials. It should be noted the ordering of time-scales for nerve cells is different from that of cardiac cells.
950
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
Fig. 2. The speed of the action potential as a function of the sodium channel conductance for different numbers of activation gates. The steady-state function for the m-gates was chosen such that the I –V relationship of the sodium current was independent of N (i.e. m ∞,N = (m ∞,3 )3/N ).
In cardiac cells, the time-scale of depolarisation is much slower due to the lower density of sodium channels. The sodium channel inactivation gates then become important, which leads to multiple wavefront solutions with different speeds (Hinch, 2002; Biktashev, 2002) of which only one is stable (Hinch, 2004). We calculate the travelling wave solution of the Hodgkin–Huxley model in the limit g 1 and show that the wave speed is proportional to g 1/(2N+2) . 2. Transition layer in travelling wave solutions We will demonstrate that in the limit g → ∞, a transition layer forms in the Hodgkin–Huxley equations. The approximation is based on a separation of time-scales, which are the same as those in the FitzHugh reduction with the exception that we do not assume that τm is rapid. The first stage of the analysis is to non-dimensionalise the equations. The non-dimensional functions Ti (v) are defined by Ti (v) = τi (v)/τi (0), where i = m, h, n. Define the non-dimensional variables τm (0)Cm t , xˆ = x , (3) tˆ = τm (0) σ the shifted m gate variable m¯ = m − m ∞,N (0),
m¯ ∞ (v) = m ∞,N (v) − m ∞,N (0),
(4)
and the non-dimensional parameters τm (0) τm (0) ,
n = , τh (0) τn (0) (5) gKτm (0)n 4∞ (0) gL τm (0) gNa τm (0)h ∞ (0)
K = ,
L = , g= . Cm Cm Cm Using the Hodgkin–Huxley parameters (Hodgkin and Huxley, 1952), we find that
m , h , n , K , L 1, g 1 and g mN K , L . The physiological interpretations of
m = m ∞,N (0),
h =
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
951
these small parameters are: the time constants of the h- and n-gates are slow compared with that of the m-gate ( h , n 1); at the resting potential the m-gate is almost completely closed ( m 1); the time-scale of the rate of change of the membrane potential due to the potassium and leakage currents at the resting potential is small compared with the time constant of the m-gate ( K , L 1); if the sodium channels were fully open the rate of change of the membrane potential would be quicker than the time constant of the m-gate (g 1); and the resting potential is maintained by balancing the potassium current with the leak current, and the sodium current is not important in determining the resting potential (g mN K , L ). Therefore, to leading order, the resting potential is independent of g. The non-dimensional equations are ∂v h n4 ∂v 2 (v − vNa ) − K 4 ¯ N (v − vK ) − L (v − vL ), = 2 − g( m + m) h ∞ (0) ∂ xˆ n ∞ (0) ∂ tˆ m¯ ∞ (v) − m¯ ∂ m¯ , = Tm (v) ∂ tˆ (6) ∂h h ∞ (v) − h , = h Th (v) ∂ tˆ ∂n n ∞ (v) − n . = n ˆ Tn (v) ∂t Next, we look for travelling wave solutions with constant profile by introducing the travelling wave co-ordinate ξ = tˆ − x/c, ˆ where c is the wave speed. We consider the asymptotic limit g → ∞ and m g → 0. We introduce the re-scaled variables ξ , c and m¯ c = Cg β ,
ξ = Xg −α ,
m¯ = Mg −γ
(7)
where α, β, γ > 0 and are determined later. The physiological meanings of these rescalings are: the wave speed is large, the depolarisation takes place in a narrow transition layer and the depolarisation occurs when the m-gates are only partially open. Inserting (7) into (6) yields dv 1 dv 2 h = g α−2β 2 (v − vNa ) − g 1−α−γ N M N dX C dX 2 h ∞ (0) + O( m g 1−α−γ (N−1) , K g −α , L g −α ), dM m¯ ∞ (v) = g γ −α + O(g −α ), dX Tm (v) dh = O( h g −α ), dX dn = O( n g −α ). dX
(8)
The final two equations tell us that h and n are constant throughout the wavefront, which is one of FitzHugh’s approximations. The leading order balance for the first two equations occurs when α = γ = 1/(1 + N) and β = 1/(2 + 2N). The first two equations in the
952
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
transition layer are 1 d2 v dv = 2 − M N (v − vNa ) + O( m g α , K g −α , L g −α ), dX C dX 2 (9) dM −α = f (v) + O(g ), dX where f (v) = m¯ ∞ (v)/Tm (v) and has the properties f (0) = 0 and f > 0 (note that to leading order v = 0 is the resting potential because g mN K , L ). The equation for v has fixed-points when M = 0 (which occurs when the tissue is at rest so v = 0) and when v = vNa . The travelling wave solution is the hetero-clinic orbit connecting these two points, so the boundary conditions for a travelling wave solution are M(−∞) = v(−∞) = 0 and v(∞) = vNa . It should be noted that at large times when v → vNa , M will grow linearly with time. At first this may seem a problem because m < 1; however, since m = M/g α then m is asymptotically small even though M is large. For long times (X = O(g α )), it is necessary to return to the un-scaled variables and to linearise around the fixed points. This solution can then be linked to the transition layer solution using matched asymptotic expansions; however, it is omitted since it adds nothing of interest because the physiological condition that v is bounded determines the boundary condition v → vNa as X → ∞ within the transition layer. The transition layer problem is an eigenvalue problem for C; however, the leading-order terms are independent of g, therefore c ∝ g 1/(2N+2) ,
(10)
which agrees with a previous analysis of a threshold model (Hinch, 2002). The simple form of (9) suggests that it may be possible to prove existence, uniqueness and stability properties of the transition layer solution. Note that one of the correction terms is of order
m g α (≈0.1 using the Hodgkin–Huxley parameters), which grows at large g. Therefore, to check that the transition layer solution (7) and (9) is asymptotic as g → ∞, it is necessary to ensure that m g α 1. This can be achieved by a minor modification of the I –V relationship of the sodium current (see Fig. 3B). The physiological justification for this modification is that it is known that the resting potential is maintained by the balancing of the potassium and leak currents and not the sodium current. When we take the asymptotic limit g → ∞, the sodium current would become important at the resting potential which would be physiologically incorrect. The modification therefore allows us to examine the asymptotic limit g → ∞ without changing the underlying physiological mechanisms. Fig. 3A shows a comparison of the wave speed of the numerical solution of (1) with the transition layer solution ((7) and (9)) on a log–log plot for different values of N. Details of the numerical schemes used are given in the Appendix. The straight lines in Fig. 3A are the transition layer solutions and the constant C was calculated numerically. Note that, at large g, the numerical solutions asymptote to the transition layer solution. Fig. 4 shows the profile of the wavefront for different values of g compared with the transition layer solution, when N = 1 (A) and N = 3 (B). The time variable is rescaled by a factor of g 1/(N+1) , so that a direct comparison of different values of g can be made. As g → ∞, the N = 1 solution asymptotes quicker to the transition layer solution than the N = 3 solution. This is because the main correction term in (9) is proportional to g −1/(N+1) , which decays
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
953
Fig. 3. A. The wave speed as a function of sodium channel conductance for different values of N . The straight lines are the transition layer solutions and the curves are the full solution of (1). Note that the full solutions asymptote to the transition layer solutions at large g. B. The modification to the I –V relation of the sodium current so that it is possible to show the transition layer solutions are asymptotic.
Fig. 4. The profile of the wavefronts when N = 1 (A) and N = 3 (B). In both graphs, the lines from left to right are: the transition layer solution calculated from (9); then gNa = 10 000 mS cm−2 (g = 1483), gNa = 1000 mS cm−2 (g = 148) and gNa = 120 mS cm−2 (g = 17.8) calculated numerically from (1). Time was rescaled by a factor of g 1/(N +1) in all cases so a direct comparison between different values of g could be made (i.e. in the limit of large g the profile will tend to the same shape). The size of the error in the asymptotic calculation is g −1/(N +1) , which takes the values 0.026, 0.082 and 0.237 when N = 1 (A), and 0.161, 0.287 and 0.487 when N = 3 (B). In both cases, the profile of the solution full equations asymptotes the transition layer solution as g is increased.
more slowly with g when N is increased. Using the Hodgkin–Huxley parameter values, we find that g −1/(N+1) ≈ 17.9−1/4 ≈ 0.49. Unfortunately the power of one quarter dramatically increases the size of the small parameter thus reducing the accuracy of the transition layer solution (Fig. 4B). In this note we have demonstrated that the FitzHugh reduction of the Hodgkin–Huxley equations for propagating nerve impulses is not valid when physiological parameters are used, and this led to an over-estimation of the wave speed by a factor of 5 (Fig. 2). An alternative asymptotic analysis was carried out which included the dynamics of the
954
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
activation gates. The results of this alternative analysis agree with numerical solutions of the Hodgkin–Huxley equations, and demonstrates the importance of including the dynamics of the activation gates when calculating the speed of propagation of nerve impulses. Acknowledgements This work was supported by The Wellcome Trust. RH would like to thank Santiago Schnell for useful comments on this manuscript. Appendix The Hodgkin–Huxley equations were solved using a standard shooting method (FitzHugh and Antosiewicz, 1959). The travelling wave co-ordinate ξ = t − x/c was introduced to (1) which produces a set of five coupled first-order ordinary differential equations (ODEs), creating an eigenvalue problem for c. The resting potential (i.e. fixedpoint) was found using the Newton–Raphson method and the Jacobian matrix was shown to have a single unstable eigenvalue. An initial guess for the eigenvalue c was made. The membrane potential was increased by 0.01 mV above the resting potential and the ODEs were integrated using the fourth-order Runge–Kutta method. The potential either diverged to +∞ implying that the trial value for c was too large, or diverged to −∞ implying that the trial value for c was too small. The eigenvalue was then found iteratively using interval bisection. This method enables the efficient calculation of the wave velocity; however, it is unable to resolve the full pulse profile (FitzHugh and Antosiewicz, 1959; Miller and Rinzel, 1981), which is a standard problem when using shooting in problems with narrow transition layers. Selective points on the graphs were checked by directly solving the partial differential equations using an explicit Euler finite difference scheme, and both methods agreed with each other. However, direct computation of the PDEs is both slow and inaccurate compared with the shooting method. References Biktashev, V.N., 2002. Dissipation of the excitation wavefronts. Phys. Rev. Lett. 89, 168102. Casten, R.G., Cohen, H., Lagerstom, P.A., 1975. Perturbation analysis of an approximation to the Hodgkin–Huxley theory. Quart. Appl. Math. 4, 365–402. FitzHugh, R.A., 1960. Thresholds and plateaus in the Hodgkin–Huxley nerve equations. J. Gen. Physiol. 43, 867–896. FitzHugh, R.A., Antosiewicz, H.A., 1959. Automated computation of nerve excitation-detailed corrections and additions. J. SIAM 7, 447–457. Hinch, R., 2002. An analytical study of the physiology and pathology of the propagation of cardiac action potentials. Prog. Biophys. Mol. Biol. 78, 45–81. Hinch, R., 2004. Stability of cardiac waves. Bull. Math. Biol. 66, 1887–1908. Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerves. J. Physiol. 117, 507–544. Hunter, P.J., McNaughton, P.A., Noble, D., 1975. Analytical models of propagation in excitable cells. Prog. Biophys. Mol. Biol. 30, 99–144. Keener, J., Sneyd, J., 1998. Mathematical Physiology. Springer.
R. Hinch / Bulletin of Mathematical Biology 67 (2005) 947–955
955
Kepler, T.B., Abbott, L.F., Marder, E., 1992. Reduction of conductance-based neuron model. Biol. Cybern. 66, 381–387. Krinsky, V.I., Kokov, Y.M., 1973. Analysis of the equations of excitable membranes reducing the Hodgkin–Huxley equations to a second order system. Biofizika 18, 506–511. Miller, R.N., Rinzel, J., 1981. The dependence of impulse propagation speed on firing, frequency, dispersion, for the Hodgkin–Huxley model. Biophys. J. 7, 227–259. Nagumo, J., Arimoto, S., Yoshizawa, S., 1962. An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061. Rinzel, J., Keller, J.B., 1973. Traveling wave solutions of a nerve conduction equation. Biophys. J. 13, 1313–1337.