INTRODUCTION

TOP ABSTRACT INTRODUCTION THE HODGKIN-HUXLEY MODEL SOLITARY PULSE CONCLUSIONS REFERENCES

Understanding the mechanisms of the propagation of nerve activity is one of the fundamental problems in biophysics. The simplest example of such a propagation is a single solitary traveling pulse of action potential in an axon (Katz, 1966). Today it is well established that the changes of the membrane potential in nerve tissue are the result of the complex dynamics of the ionic currents through voltage-sensitive channels (Katz, 1966). The first detailed quantitative measurements of the ionic currents were performed by Hodgkin and Huxley in the early 1950s (Hodgkin and Huxley, 1952). By using the voltage clamp technique they were able to measure the kinetics of Na⁺ and K⁺ currents in the squid giant axon. This led them to a set of differential equations that describe the dynamics of the action potential. Furthermore, by combining these equations with the cable equation for spreading of the current in the axon they were able to calculate the shape and velocity of the propagating action potentials (Hodgkin and Huxley, 1952; Huxley, 1959). The predictions of their model turned out to be in remarkably good agreement with the experimental observations.

The reason that the model introduced by Hodgkin and Huxley (the HH model) admits quantitative comparisons with the experiments is that it contains detailed information about the voltage-dependent kinetics of the Na⁺ and K⁺ channels. Naturally, this makes the models quite complex and intractable analytically. So far, the basic tools for studying the HH model have been numerical simulations. Note that because of the HH model's complexity it was not until recently, with the advent of very fast computers, that simulations could be done routinely. Even then, one is still required to do simulations for each set of the parameters of the model. Therefore, analytical studies that give functional dependencies of the main parameters of the action potentials on the parameters of the model are still highly desirable.

The early analytical works on the HH model relied on the strong separation of the time scales of the (fast) activation and (slow) inactivation processes. These studies made an assumption that the Na⁺ activation is the fastest process and can be eliminated adiabatically, which amounts to assuming that the sodium activation variable m = m(V), where m(V) is the resting value of m at a given membrane voltage V (FitzHugh, 1961; Casten et al., 1975; Carpenter, 1977, 1979). This leads to a cubic-like nonlinearity in the equation for the membrane potential. By further assuming that the Na⁺ inactivation and K⁺ dynamics are much slower than the Na⁺ activation, the problem of the action potential propagation reduces to a single reaction-diffusion equation for the front of the action potential (Casten et al., 1975). A number of simpler models (FitzHugh-Nagumo type) with similar properties had been introduced to mimic the behavior of the membrane (FitzHugh, 1961; Nagumo et al., 1964; Rinzel and Keller, 1973; Casten et al., 1975; Jones et al., 1991). The latter, in fact, became quite popular for explaining traveling wave phenomena in a variety of excitable systems in physics, chemistry, and biology (Vasiliev et al., 1987; Murray, 1989; Mikhailov, 1990; Kerner and Osipov, 1994).

Although this kind of analysis leads to a qualitative explanation of the excitability of the nerve membrane, it fails to give any quantitative predictions for the speed of the propagating action potentials. It also predicts that the traveling wave should have the form of a broad excitation region with the sharp front and back. This is in contrast to the observations in which the pulse is a narrow localized region of excitation (a spike). The reason for this is that in reality it is typically the membrane potential rather than the Na⁺ activation that is the fastest process. For example, in the squid giant axon the time constant of the membrane potential change is _V ~ 0.01 ms, whereas the time constant of the Na⁺ activation is roughly _Na ~ 0.2 ms. Therefore, the FitzHugh-Nagumo-type models are in fact not adequate for any quantitative predictions of the characteristics of the action potential. Also, they only qualitatively reveal the mechanism of the wave propagation.

In this paper we introduce an approximation scheme that does take into account this relationship between the time scales. We will construct an approximate solution for a single traveling pulse in the HH model that is in quantitative agreement with the solutions of the full HH model. We will investigate the structure of the front of the traveling pulse and show that it is substantially different from the conventional case of the FitzHugh-Nagumo-type models. We will also obtain an explicit analytical expression for the speed of the pulses that agrees with the results of the simulations of the HH model within 20% accuracy in a wide parameter range. Using the obtained solutions, we will construct an approximate solution for the entire pulse that is also in quantitative agreement with the solutions of the full HH model.

	THE HODGKIN-HUXLEY MODEL

TOP ABSTRACT INTRODUCTION THE HODGKIN-HUXLEY MODEL SOLITARY PULSE CONCLUSIONS REFERENCES

In the following we will use the version of the HH model that was originally introduced by Hodgkin and Huxley to study the behavior of the squid giant axon (Hodgkin and Huxley, 1952) and later adopted by many researchers as a benchmark of the models of nerve activity. Namely, we will consider the following equations

(1)

(2)

(3)

(4)

Here, as usual, Eq. 1 is the cable equation for the membrane potential V, with C = 1 µF/cm² the membrane capacitance per unit area, a = 238 µm is the radius of the axon, and  = 35.4  × cm the resistivity of the intracellular space; g_Na = 120 m¹/cm², g_K = 36 m¹/cm² are the conductances of the open Na⁺ and K⁺ channels per unit area; V_Na = 115 mV and V_K = 12 mV are the equilibrium potentials of Na⁺ and K⁺, and g_l = 0.3 m¹/cm² and V_l = 10.5989 mV are the leakage conductance per unit area and the leakage voltage, respectively. With these definitions the resting potential V_r = 0. Similarly, m and h are the activation and the inactivation variables for the Na⁺ channels, respectively; n is the K⁺ inactivation variable; the rates _m,h,n and _m,h,n as functions of V at temperature T = 6.3°C can be found in Hodgkin and Huxley (1952) (note that Hodgkin and Huxley (1952) have an opposite sign convention for V). The temperature changes are accounted for by a factor  = 3^(T6.3)/10 multiplying all  and  values; T is in degrees Celsius. The lengths are measured in centimeters and the times in milliseconds. The voltage V is measured in millivolts.

Equations 1-4 constitute a closed system of partial differential equations that quantitatively describes the changes in the membrane as functions of time and space. Let us emphasize that their ingredients are obtained by measurements and fitting of the parameters to the actual experiments, so it is important to understand the relationships between the characteristic parameters, namely the time scales, in this system. From the functional form of  and  values (Hodgkin and Huxley, 1952) we can make the following estimates for the time constants _m,h,n for m, h, and n, respectively, at T = 6.3°C:

(5)

(6)

(7)

Also, from Eq. 1 one gets the following estimate for the time scale _V of the variation of the voltage, assuming that all the Na⁺ channels are open:

(8)

One can see that the following hierarchy of time scales holds in the system:

(9)

The first inequality holds better for sufficiently low temperatures and remains qualitatively correct up to the temperatures T ~ 30°C, at which the pulses fail to propagate in the HH model (Huxley, 1959). As we pointed out in the Introduction, this is an important property of the system which is not taken into account in most of the analytical studies of the HH model. In the following, we will use this hierarchy of time scales to introduce the approximation scheme for the traveling pulses in this model.

Another important point about the HH model is the fact that the very nonlinearities in Eq. 1, namely the powers with which the variables m, h, and n enter the equation, are determined experimentally (Hodgkin and Huxley, 1952). Furthermore, these powers correspond to the number of particles involved in the operation of the respective channels and therefore represent significant physical quantities. As will be seen below, these powers play crucial roles in our studies.

Before going to the analysis of the traveling pulses, let us discuss the basic physics of the excitability in the HH model. In the rest state, the Na⁺ channels are basically closed; at V = 0 the equilibrium values for m and h are m₀  0.05 and h₀  0.60, respectively, while the K⁺ channels are partially open, with n₀  0.32. If, by applying an external stimulus, the membrane voltage V is increased to ~10 mV, the Na⁺ channels will start opening on the time scale of order _m. The influx of the Na⁺ ions will in turn lead to the increase of the membrane potential V on the time scale intermediate between _V and _m (see below), resulting in positive feedback. The membrane potential V will come close to the resting potential V_Na, while the Na⁺ channels will become mostly open with m  1. During this time, the slow inactivation variables h and n will remain almost unchanged. After that, the slow inactivation variables h and n will start to react, closing the Na⁺ and opening the K⁺ channels, which will drive the potential V back to equilibrium. In the spatially extended system the diffusive spreading of the current in front of the excitation region in the axon will provide the sustaining force for the propagation of the pulse along the axon. In that sense, from the physical point of view the traveling pulse in the nerve axon is a classical example of an autosolitonself-sustained solitary inhomogeneous state in an active dissipative system whose existence is determined only by the nonlinearities of the system and not the initial conditions (Kerner and Osipov, 1994).

	SOLITARY PULSE

TOP ABSTRACT INTRODUCTION THE HODGKIN-HUXLEY MODEL SOLITARY PULSE CONCLUSIONS REFERENCES

We are now going to construct an approximate traveling wave solution in the form of a solitary pulse, using the ideas introduced in the preceding section. Let us introduce a self-similar variable z = x  ct, where c is the propagation speed of the pulse. Then, Eqs. 1-4 for a traveling wave with speed c will become

(10)

(11)

(12)

(13)

The boundary conditions for these equations are

(14)

where m₀, h₀, and n₀ are the values of m, n, and h in the rest state, respectively. For the chosen functions  and  the rest state V = 0 is unique and stable.

The solution of Eqs. 10-13 in the form of a traveling solitary pulse obtained numerically at the "standard" temperature T = 6.3°C is shown in Fig. 1. From this figure one can see several features of the solution we will use in the approximation scheme that we are going to construct. First, observe that the length scale of the rise of the potential is substantially smaller than that of the fall of the potential. Second, during the rise of the potential the variables h and n remain almost unchanged at their resting values h₀ and n₀. Third, in front of the spike the value of m (that is, m₀) is practically zero.

View larger version (18K): [in this window] [in a new window]   FIGURE 1   The numerical solution of Eqs. 10-13 at T = 6.3°C.

Let us use the above facts to simplify Eqs. 10-13. Because the values of h and n change little in the front of the spike, we may replace them by their values h₀ and n₀ at rest and disregard Eqs. 12 and 13. Furthermore, because the value of m₀ is very small, with very good accuracy, we may assume it to be zero. Therefore, in the rest state we will have g_Kn₀⁴V_K + g_lV_l = 0 with very good accuracy, so these terms drop out of Eq. 10. Also, the coefficient multiplying V in the last two terms of Eq. 10 is of order 0.7 and is much smaller than the contribution from the term g_Nam³h during the rise of the potential when m is not close to zero, so these terms can be dropped as well. What we are then left with is an equation for V coupled only to the equation for m with a number of terms dropped. Observe that because V is much faster than m when m ~ 1, the value of m has to be sufficiently small for the nontrivial collective dynamics of V and m to be possible.

This allows further simplification of Eq. 11 by neglecting the terms proportional to m. After making all these approximations, we are left with the following set of equations

(15)

(16)

instead of Eqs. 10 and 11. Note that we added a term _m(0) to Eq. 16 for this equation to be consistent with the approximate boundary conditions ahead of the spike front

(17)

where V_z = dV/dz. We can do this in our approximation scheme because the value of _m(0) is in practice very small compared to _m(V_Na).

Let us assume that the characteristic value of m in the front is 1 and the characteristic width of the front is l. Then, as all the terms in Eqs. 15 and 16 should be of the same order of magnitude, we obtain the following estimates:

(18)

From these one can also estimate the characteristic time scale for the rise of the potential in the pulse as  ~ l/c ~ (C/_m³g_Nah₀)^1/4, where

(19)

so

(20)

One can see from this equation that the dynamics in the front of the traveling pulse will indeed occur on the time scale intermediate between _V and _m.

From the estimates above we immediately conclude that in the traveling spike

(21)

where is a constant of order 1. Substituting the parameters of the HH model at T = 6.3°C, we see that   0.6, what corresponds to the relevant quantity ³  0.2, which is indeed rather small.

Let us introduce the following new variables:

(22)

where

(23)

The latter is plotted in Fig. 2. Note that this way the time scale in Eq. 16 has been absorbed into the constant _m.

View larger version (9K): [in this window] [in a new window]   FIGURE 2   The dependence (u).

In terms of the variables introduced in Eq. 22 and after a few manipulations one can rewrite Eqs. 15 and 16 in the following form:

(24)

(25)

(26)

where now s is an independent variable. These transformations can be done for 0 < u < 1 because (u) is always positive for u  0 (see Fig. 2). Note that these equations do not have any  dependence in their right-hand side, so it suffices to solve only Eqs. 24 and 25. The solution for (s) can then be obtained by a simple integration.

The problem now became substantially simpler because instead of solving the nonlinear boundary value problem for Eqs. 10-13, one now needs to solve the initial value problem for Eqs. 24 and 25. Indeed, according to Eq. 17, when z  + we have m  0, so s  0 as   +. This means that u = 0 and v = 0 at s = 0, because du/d = v  0 as   + (see Eqs. 17, 24, and 26). One should be careful to specify what exactly happens near s = 0, because (0) = 0. To do this, let us divide Eq. 25 by Eq. 24. We get

(27)

When s  0, we have dv/du  1, provided that v(s) has a non-zero derivative at s = 0 (the latter follows from the physical considerations). Therefore, according to Eqs. 24 and 25, as s  0, the solution will behave like

(28)

where the prime means differentiation. Here we expanded (u) in the neighborhood of zero and took into account that '(0)  0.

It is not difficult to see from Eq. 27 that for 0 < u < 1 and v > 0 the projection of the phase trajectory on the u-v plane will lie below the line u = v. Because du/ds > 0 and there are no fixed points in this region of u and v, the solution u(s), v(s) will cross either the line u = 1 or v = 0 in the u-v plane. By direct inspection of Eqs. 24 and 25 one can see that this intersection is transversal. Observe that the intersection of the lines u = 1 and v = 0 is a fixed point in this plane.

According to Eqs. 24 and 25, once the solution leaves the region bounded by the lines u = v, u = 1, and v = 0, it can never come back to this region. Indeed, if the solution crosses the line u = 1 at some v > 0 in the u-v plane, it will then have dv/ds > 0 for all s, so v will stay positive. If, however, the solution crosses the line v = 0 at some u < 1, we will have both dv/ds < 0 and du/ds < 0 afterward. In fact, it is easy to see that the only trajectory on which u and v will remain bounded for all s is the one that connects the point u = 0, v = 0 with u = 1, v = 0. Therefore, this trajectory is precisely the one that corresponds to the front of the traveling pulse.

Of course, not all the speeds will produce this kind of trajectory. It is clear that if is very large, the trajectory will cross the line u = 1 in the u-v plane at v close to 1. However, if is very small, the trajectory will cross the line v = 0 at very small u. Fig. 3 shows the results of the numerical solution of Eqs. 24 and 25 at different values of . From this numerical solution we found that the trajectory that connects u = 0, v = 0 and u = 1, v = 1 exists only for a unique value of  = *.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   The numerical solution of Eqs. 24 and 25 in the u-v plane v(u) at different .

In fact, it is possible to prove that such a trajectory indeed exists and is unique at a unique value of c. To do this, let us see what happens with the trajectory as the value of is changed. For convenience, we will use u instead of s as an independent variable. Let v₀(u) and s₀(u) be a trajectory in the region bounded by u = v, u = 1, and v = 0 with the initial conditions v₀(0) = 0, s₀(0) = 0 for some  = ₀. To calculate the changes in the trajectory v(u), s(u) as is changed by , we write the equations in variations for v and s obtained from Eqs. 24 and 25

(29)

(30)

Because the change in does not affect the initial conditions, we should have

(31)

Note that according to Eqs. 28 we have v₀ ~ u and s₀ ~ u, so s ~ u³ and v ~ u³ in the neighborhood of u = 0.

According to Eq. 30, when u  0 we have (d/du) v > 0 for  > 0, so v > 0. In turn, according to Eq. 29, (d/du) s < 0 and therefore s < 0. This means that the derivatives (d/du) v and (d/du) s do not change signs, so v > 0 everywhere for  > 0 and vice versa. Therefore, when the value of changes from 0 to , the point at which the trajectory crosses either the line u = 1 or the line v = 0 in the u-v plane will go monotonically from u = 0, v = 0 to u = 1, v = 1. Because it depends continuously on , there is a unique value of  = *, at which this point coincides with u = 1, v = 0. Numerically, the value of * =  up to the fourth digit. Thus, the dynamics in the pulse front uniquely determines its propagation speed.

Thus, we have obtained an approximate analytical expression for the speed of the traveling pulses in the HH model:

(32)

Equation 32 predicts the speed of the traveling pulse as a function of the parameters. To compare this predicted speed with the results obtained from the numerical solution of the HH model, we plot the speed c as a function of temperature obtained from Eq. 32 and from the numerical simulations of the HH model (see also Huxley, 1959) in Fig. 4 (recall that the temperature dependence is contained in the value of _m).

View larger version (11K): [in this window] [in a new window]   FIGURE 4   The speed c of the traveling pulse as a function of T. The solid line shows the results of the numerical solution of the full HH model. The dashed line is the prediction of Eq. 32. The dotted line is the result of the solution of the HH model without the h and n dynamics.

As can be seen from this figure, the approximate expression for the speed of the pulse given by Eq. 32 agrees with the results for the HH model within 20% at temperatures below 15°C. We emphasize that this is the kind of accuracy with which the HH model itself predicts the speeds of the traveling pulses as compared to the experiments. At higher temperatures the agreement between Eq. 32 and the results of the simulations of the HH model becomes worse, and at the temperatures of the propagation threshold Eq. 32 fails completely. We have also checked that Eq. 32 predicts the correct dependences on the other parameters with similar accuracy at low enough temperatures. For example, Fig. 5 shows a comparison of the prediction of Eq. 32 with the results of the numerical simulations of the HH model at T = 6.3°C as the value of the membrane capacitance C per unit area is varied on the log-log plot. Note that the slopes of the two graphs in Fig. 5 agree very well with each other. The agreement of the slopes is not as good for the log-log plot of the dependence of c on g_Na with other parameters fixed. This is probably the consequence of the fact that the errors introduced by our approximation depend on g_Na stronger than the prediction of the approximation c ~ g_Na^1/8.

View larger version (10K): [in this window] [in a new window]   FIGURE 5   The speed c of the traveling pulse as a function of the membrane capacitance C. The solid line is the result of the numerical solution of the HH model, the dashed line is the prediction of Eq. 32.

Incidentally, if the factor of 2/3 in Eq. 32 is replaced by 5/9, it will give the speed of the pulse within a few percent of that found in the full HH model for T < 15°C. Note that if one assumes that m is the fastest variable (FitzHugh, 1961; Casten et al., 1975; Carpenter, 1977, 1979) and calculates the speed of the traveling wave, one will get a value an order of magnitude greater than the actual value.

Observe that Fig. 4 also shows the dependence of the speed of the pulse on temperature obtained from the simulations of the HH model without the h and n dynamics (dotted line). Note that the insignificance of these variables is one of the key assumptions in deriving Eq. 32. One can see that this solution gives an even better approximation to the speed of the pulse. This problem, however, cannot be treated analytically in the same manner as that for Eqs. 15 and 16.

Let us emphasize that the existence of the front solution is essentially determined by the complicated interplay of the V and m kinetics, so the problem does not reduce to simple phase plane analysis, in contrast to most studies of the traveling waves in excitable systems (FitzHugh, 1961; Rinzel and Keller, 1973; Casten et al., 1975; Carpenter, 1977, 1979; Vasiliev et al., 1987; Mikhailov, 1990). Note that a similar situation takes place in a class of excitable systems in which the so-called spike autosolitons are realized (Osipov and Muratov, 1995; Muratov and Osipov, 1999). These models also give rise to the complicated front structures that are similar to the one realized in the HH model.

The validity of the approximations made by us is violated in two cases. First, when the temperature becomes sufficiently high, the dynamics of the m variable becomes faster, so the separation of the time scales _m and _V used in our arguments will no longer be justified. One of the implications of the absence of this scale separation is the fact that the characteristic value of m =  in the front can no longer be treated as small. This allows us to derive a criterion for the validity of our approximations

(33)

which is equivalent to   1 (see Eq. 21). In the case of the squid giant axon this criterion shows the applicability of our results up to T  25°C, in good agreement with Fig. 4.

Another problem may occur when the temperature becomes too low and the variable m too slow. In this case the effective time scale  of the dynamics in the front of the pulse slows down (see Eq. 20), so at some point it may become comparable to the leakage time scale _l ~ C/g_l ~ 3 ms. In this case one can no longer discard the leakage and the K⁺ contributions to the membrane current in Eq. 10. Thus, the second validity criterion becomes (see Eq. 18)

(34)

For the squid giant axon, this quantity becomes comparable to 1 only for the unrealistically low temperatures T  30°C.

As can be seen from Eqs. 33 and 34, the quality of the approximations used by us should increase with the increase of g_Na. In fact, our procedure for finding the spike's speed and the front profile is the leading order of the asymptotic limit m₀  0 and g_Na  .

Fig. 6 shows the functions u() and s() for  = * obtained numerically. One can see that u() has the form of a front connecting the rest state u = 0 with the excited state u = 1, which in the original variables corresponds to the saturation value V = V_Na.

View larger version (11K): [in this window] [in a new window]   FIGURE 6   The functions u() and s() obtained numerically from Eqs. 24-26 for  = *.

As can be seen from Fig. 6, the distribution s() behind the front approaches a straight line with the slope *. In terms of the original variables, this should correspond to the unlimited growth of m behind the front. This, however, should not be a problem because this happens only when m ~  1. When m becomes of order 1, the approximations used to derive Eq. 16 ceases to be valid. However, when this happens, V should already be very close to V_Na, so on the time scale _{m h,n} the variable m will simply exponentially relax to m = m (V_Na)  1 behind the front, where, as usual

(35)

This will happen at distances of order c_m l, as _m  (see Eq. 20). If we assume that on the time scale _m the front was located at z = 0, the distribution of m that properly matches with that in Fig. 6 will be

(36)

As was already noted, in the back of the spike and in the refractory tail the voltage V changes substantially slower than in the front. Because this happens when m ~ 1, the voltage is indeed the fastest variable, so we can eliminate it adiabatically from the equations. If we put all the derivatives in Eq. 10 to zero, we find

(37)

This expression uniquely determines the value of V as a function of m, h, and n.

To find the approximate distributions of all the variables in the back and behind the spike we simply need to solve the initial value problem given by Eqs. 11-13 with c given by Eq. 32 and the following initial conditions:

(38)

where we assumed that on the even larger length scale c_h,n the front is located at z = 0. This initial value problem can be straightforwardly solved numerically. The result of this solution is shown in Fig. 7. Note that the changes in temperature will only change the characteristic length of this solution, not its shape.

View larger version (17K): [in this window] [in a new window]   FIGURE 7   The numerical solution of Eqs. 11-13, 37, and 38.

One can simplify the procedure of finding the distributions of m, h, and n by using the fact that _{m h,n} by adiabatically eliminating m. This will amount to replacing m by m(V) from Eq. 35 in Eq. 37 and then solving for V as a function of h and n. The result of the numerical solution of Eqs. 12 and 13 with these approximations is shown in Fig. 8. This figure shows a good agreement of the slow variables h and n in the refractory tail. Also, observe an abrupt jump in the back of the spike. This is due to the fact that now the value of V is not uniquely determined by h and n, so at some point in the solution a jump occurs from one branch of the dependence V(h, n) to the other (see also Carpenter (1977, 1979)). Note that while the adiabatic elimination of m works well for the refractory tail, it fails in the back of the pulse (compare Figs. 7 and 8).

View larger version (17K): [in this window] [in a new window]   FIGURE 8   The numerical solution of Eqs. 12, 13, 37, and 38 with m = m(V) given by Eq. 35.

The results in Figs. 6, 7, and Eq. 36 can be combined to give a quantitative approximation for the whole pulse. This is done in Fig. 9 for T = 6.3°C. One can see a remarkable similarity between the solution of the full HH model shown in Fig. 1 and the one shown in Fig. 9. Thus, our approximation scheme has been able to quantitatively capture the essential features of the traveling pulses in the HH model.

View larger version (18K): [in this window] [in a new window]   FIGURE 9   An approximate solution for the entire pulse.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THE HODGKIN-HUXLEY MODEL SOLITARY PULSE CONCLUSIONS REFERENCES

We have developed a method that allows approximate computation of the shape and parameters of the traveling spikes in the HH model of conductance along an axon. Our method is different from the conventional approach (FitzHugh, 1961; Rinzel and Keller, 1973; Casten et al., 1975; Carpenter, 1977, 1979) in the fact that it treats the membrane potential, rather than the sodium activation variable, as the fast variable. We show that this is in fact the case for the typical set of the parameters of the Hodgkin-Huxley model. This leads to a good quantitative agreement between the predictions of the theory and the results of the numerical simulations of the HH model.

Let us emphasize that the HH model itself gives only an approximate, although quite accurate, description of the dynamics of the action potential in an actual axon. What we find is that in a broad range of the parameters the approximation introduced by us gives an error that is comparable to the error produced by the HH model itself, as opposed to the experiments. For example, at T = 18.5°C the speed of the pulse in the squid giant axon was found to be 21.2 m/s (Hodgkin and Huxley, 1952). The direct numerical simulation of the HH model produces the speed of 18.8 m/s, while our procedure, which for this temperature is already near the limit of its applicability, gives 24.7 m/s. Therefore, we suggest that the ideas of our analysis can be built into simpler and more tractable models of nerve conductance that will yet be able to give quantitative agreement with the experimental observations.

One of the observations from the analysis made by us is the fact that the speed of the traveling spikes in the HH model depends very weakly on the slow-state variables of the membrane. Indeed, according to Eq. 32, the speed of the spike is independent of the value of n in front of the spike and is proportional to h₀^1/8, so a change by a factor of 2 in h₀ will result in only a 10% change in the speed. This makes a perfect biological sense. Thus, propagation of the nerve pulses is indeed a very robust phenomenon.

Another observation one can make from Eq. 32 is that if one assumes that in addition to the membrane conductance C there is an extra capacitance associated with each sodium channel, there exists a density of the channels at which the speed is maximal (Hodgkin, 1975). Indeed, let us assume that C = C₀ + NC^*_Na and g_Na = Ng^*_Na, where C₀ is the capacitance of the membrane without the channels, N is the channel density, C^*_Na is the capacitance associated with a single channel, and g^*_Na is the maximum conductance of a single channel. For the squid giant axon these parameters are estimated to be C₀ = 0.8 µF/cm², C^*_Na = 4 × 10¹⁸ F, g^*_Na = 2.4 × 10^{12 1}, and N = 500 µm² (Hodgkin, 1975). Substituting these expressions into Eq. 32, one gets the speed of the pulse as a function of N. It is not difficult to see that this function has a maximum at N = N_max = C₀/(4C^*_Na). For the numerical values above we find N_max = 500 µm², which suggests that the channel density in the axon is indeed at the optimum level. The fact that we get exactly the same value of N as the one observed may be a bit fortuitous because of the approximate nature of Eq. 32. Note that because of the very slow dependence of the speed on g_Na, the maximum of the dependence c(N) is in fact very flat, so a change of N by a factor of 2 from N_max results only in a 7% difference in c.

So far, we were talking only about the traveling wave solutions in the form of the solitary spikes. It is not difficult to see that our method can be extended to spike trains as well. Indeed, the speed of a spike in a spike train is determined by the value of the slow variable h in front of the spike (see Eq. 32), which, however, is now different from the equilibrium value h₀ and must be determined. Outside of the spike fronts one has to solve the equations of the slow dynamics given by Eqs. 11-13 in which the fast variable V has been eliminated adiabatically via Eq. 37. These equations have to be solved as an initial value problem for L  z  0 with m(0) = m(V_Na), h(0) = h_s, and n(0) = n_s. Here h_s and n_s are the values of h and n in the spike, respectively, L is the spatial period of the spike train, and we assumed that the front of one of the spikes in the spike train is located at z = 0. The spikes are also assumed to travel to the right with the speed given by Eq. 32, in which h₀ is replaced by h_s. Then, the values of h_s and n_s have to be found self-consistently from the condition that h(L) = h_s and n(L) = n_s.

We implemented this procedure numerically. To find the values of h_s and n_s as functions of L, we started with a reasonable initial guess for h_s and n_s and then solved the initial value problem up to z = L. The values of h and n at z = L were then taken as the new values for h_s and n_s, respectively, and the procedure was iterated. We found a fast convergence to the periodic solution. Then, from the value of h_s we obtained the speed of the spike train as a function of L and therefore the dependence of the spike frequency on the period. Our main finding was that in the region of the applicability of our approximation (Eq. 33) the speed of the spike trains is practically independent of the period, so they have almost no dispersion. This also makes good biological sense. In fact, the amount of the dispersion we found is comparable to the error introduced by our approximation scheme. Because we are only interested in quantitative predictions, we do not present these results in detail here. Also, our method fails for periods L  10 cm, for which a substantial amount of dispersion was found in the simulations of the full HH model (Miller and Rinzel, 1981). Nevertheless, let us point out that the results obtained with our method are in good qualitative agreement with those obtained for the full HH model (Miller and Rinzel, 1981). In particular, according to our numerical solution outlined above there exists a period of the spike train for which the speed of the spikes reaches maximum, greater than that of the solitary spike due to a slight overshoot of the h variable behind the spike. However, the magnitude of this overshoot is so small that it only changes the speed of the spike by a fraction of a percent, so for practical purposes the spike trains with period L  10 cm may be considered dispersionless.

In short, we have introduced an approximation scheme that allows making quantitative predictions of the shape and the parameters of the traveling pulses in the HH model of nerve conductance. We hope that our results will provide an easy and convenient tool for analyzing the fascinating complexity of neural activity.