All published worksheets from http://sagenb.org
Image: ubuntu2004
Choose two out of the following three exercises on the Hodgkin-Huxley (HH) model. You may use the XPP or Matlab codes on our course website. Choose a or for numerical integration so that tightening the criterion gives comparable results. State the error criteria and integration method that you use.
Problem 1
Do some simulations with the standard HH model; set the temperature to C.
- Use an amplitude for a -ms duration that produces a spike, and a weaker stimulus that produces a just-subthreshold response. Show time course plots of and of , , and (these three on one set of axes), and show phase plane projections: vs. and vs. . Describe the features of the response as viewed in the phase plane.
- Use a long duration , say, ms, to see repetitive firing for A/cm. Produce corresponding plots as for (a) and describe the features.
- Develop a three-variable HH model by setting . Redo (a) and (b), and compare your results to the full model. Describe any differences.
- Make a two-variable HH model by setting and . Redo (a) and (b), and compare your results to the full model. Describe any differences.
- If you use XPP, use the command to compute the steady state and stability properties (including the eigenvalues) of the three versions of the model for and A/cm. Compare and discuss the results.
Solution
We implement the core functions for the required HH models below in Python, making use of Cython for speed, following the Matlab codes provided.
This is followed by a general driver, which can be used to run simulations of the full, three-, and two-variable models, as desired.
A note on the numerics: we use the solver , as wrapped in SciPy, for numerical integration. As performs automatic timestep adjustment with regard to numerical stability, we need only specify the timestep with regard to time resolution. To do this, we use the criterion where , i.e., is the vector of relative derivatives; and where and are parameters.
To initialize the models, we use the parameters from the Matlab codes; the initial condition is set by assuming a resting potential of mV and setting , , and to their corresponding equilibrium values at . For the applied current, we use the wave form
The first simulation (a) is of the full model for a spike in response to a short stimulus. To do this, we use A/cm and ms. For resolution, we set and .
As desired, we obtain a spike, with characteristic time courses for both the membrane potential and the gating variables , , and . The spiking behavior is especially evident in the - phase plane, which produces a full cycle of the relaxation dynamics. Futhermore, interestingly, we observe a strong linear dependence between the two variables in the - phase plane, suggesting a possible model reduction as we will see later.
The next simulation is of the full model for a short subthreshold stimulation ( A/cm, ms).
Indeed, we do not produce a spike, with the membrane voltage remaining around mV, and the gating variables more or less constant throughout. In the - phase plane, we see an attempt by the system to leave the rest state, as in (a), but it fails to cross the threshold and quickly restabilizes. The - phase plane is similar to that previously obtained.
We now run a simulation (b) of the full model in response to a long stimulus ( A/cm, ms).
Repetitive firing is observed, though, interestingly, the first spike appears larger than those ensuing. Biochemically, this could be due to, e.g., synaptic depression, though this, of course, does not appear explicitly in the model. The phase planes show, generically, multiple passes over the trajectories in (a), with the system locking onto a limit cycle sustained by the elevated as evident in the - phase plane.
The fourth simulation (c) is of the three-variable model, obtained by setting . For the remainder of this problem, we consider only the short and long spike-generating stimuli ( A/cm with ms and ms, respectively). For the short stimulus, we obtain the following:
The plots are generally similar to those obtained in (a), suggesting a justified model reduction, though opens and hence allows the cell to fire a bit sooner. The time courses of and are, in response, a bit sharper to account for the instantaneous relaxation introduced in . Furthermore, the - phase plane appears rotated relative to that in (a), no doubt a result of the instantaneous dynamics.
Likewise, the long stimulus gives:
which, again, is quite similar to that obtained in (b). Similar comments apply as in the short stimulus case. Moreover, effective synaptic depression appears to be weaker.
Finally, we run simulations (d) of the two-variable model, obtained by further setting , as suggested by essentially all of the above - phase planes.
For the short stimulus, we obtain fairly good agreement with previous results, the primary point of contention (apart from the obvious collapse of the - phase plane) being that the time course appears segmented on the downstroke. Note that the - phase plane produced here is an artifact of the matplotlib graphics system; inspection of the actual data reveals a sustained line much like that to appear below.
The long stimulus produces results very similar to those in (c). Most striking, though, is the loss of variability in the limit cycle in the - phase plane.
We study the steady states (e) of the three models above by numerically solving for , where is the instantaneous current through the membrane, on setting and , and for the full and three-variable models, and for the two-variable model. The stabilities of the steady states are evaluated by symbolically constructing the Jacobian of the system and computing its eigenvalues at the steady states.
The steady-state potential satisfying is obtained using MINPACK, and the eigenvalues of computed using LAPACK, both as wrapped by NumPy/SciPy.
We compute the steady states and their associated eigenvalues for and A/cm. For the full model, these are:
i.e., the steady state at A/cm is stable, while that at A/cm is unstable. There is evidently only one steady state as a later plot will show. The results for the three-variable model are similar:
For both of these models, therefore, the steady applied current induces a bifurcation from a steady rest state to an oscillating limit cycle, producing repetitive firing as has already been seen earlier.
For the two-variable model, however, we get something a little different. In particular, there are now three steady states: first,
second:
and third:
The first steady state, again, shares the same stable-unstable transition (Hopf bifurcation) as the other two models, while the second and third steady states are inherently unstable, generating natural thresholds for spiking behavior. Note that the first and second activated ( A/cm) steady states may be switched, judging based on the value of , though the character of the eigenvalues may suggest otherwise.
For justification of the number of steady states that each model admits, we appeal to the following graphical analysis, as obtained by plotting with , , and set accordingly at their steady-state values.
The monotonic nature of the curves for the full and three-variable models gives clearly that only one steady state exists, whereas the more complicated shape of that of the two-variable model suggests the existence of three steady states.
Problem 2
Consider HH without (i.e., ). Show that with adjustments in (and maybe ), the HH model is still excitable and generates an action potential [do it with ]. Study this two-variable (-) model in the phase plane: nullclines, stability of rest state, trajectories, etc. Then consider a range of to see if you get repetitive firing. Compute the frequency vs. relation; study in the phase plane. Do analysis to see that the rest point must be on the middle branch to get a limit cycle.
Solution
The model setup is exactly that of the three-variable model but with , hence we will make use of the three-variable model in the following. To facilitate analysis, we first implement code for computing nullclines in the - phase plane.
We then run the three-variable model with zero input ( A/cm), using the previous values of and .
As observed, the system spontaneously generates one spike, then stabilizes to a depolarized state. This may clearly be attributed to the lack of potassium outflux to repolarize the membrane. Hence to restore proper spiking behavior, we must intuitively decrease and increase . We do this by setting and , where, for illustration, we choose and , and use and in place of and . For example, using a stimulus with A/cm and ms gives repetitive firing:
The phase plane is particular illuminating, demonstrating the limit cycle about the unstable steady state. We note here that we have plotted only the nullclines as the stability of the rest state is obvious given these and the trajectory. Furthermore, that repetitive firing requires the rest state to be on the middle downward branch of the -nullcline is also immediate given the diagram. To see this, we adopt a geometric view, and suppose for simplicity that all trajectories start from rest, i.e., the upper left of the above phase plane; moreover, we assume that the dynamics in are fast, as supported by the plots. Then if the steady state lies on the rightmost branch of the -nullcline, the trajectory cannot cross and so hugs the -nullcline until it collides with the steady state, muting any potential oscillation. Similarly, if the steady state lies on the leftmost branch, the same argument applies after the trajectory crosses over past the bend of the middle branch. Therefore, a limit cycle is possible only if the steady state lies on the middle branch, in which case, we have, generically, oscillations between the left and right branches of the -nullcline, with each escape made possible by the corresponding bends defining the extent of the middle branch. To compute the frequency for a given repetitive firing trajectory, we make use of the fast Fourier transform (FFT). The following function uniformly resamples the data, then uses the FFT to pick out the dominant .
Initial numerical explorations indicate that repetitive firing is observed from to A/cm, so we use this as our parameter range.
The frequency relation is shown below, from which we see the characteristic upward trend indicative, physiologically, of a rate code (with now identified with ).
Problem 3
Convert the HH model into "phasic" mode. By "phasic" I mean that the point-neuron does not fire repetitively for any values—only one to a few spikes at onset and then it settles to a stable steady state. Many neurons in the auditory system behave phasically. Do this by, say, sliding some channel-gating dynamics along the -axis (probably just for ). [If you slide , then you must also slide .] If it can be done using and , then do the phase plane analysis.
Solution
We begin by implementing modified versions of the HH models in Problem 1, taking as an additional argument the voltage shift , manifested in the dynamics of the potassium activation gating variable (to modify ) as i.e., we slide and .
Numerical simulation then shows that taking mV suffices for phasic conversion.
In anticipation of the corresponding two-variable model, we implement the following code for plotting nullclines in the - phase plane.
There are four roots of to the equation for and . To demonstrate that we have selected the correct root in the code above, we check which of the roots satisfies at standard resting values:
We now first study the two-variable model with mV, necessary for proper interpretation by comparison later, with a sustained applied current of A/cm and ms.
In this case, we obtain repetitive firing as demonstrated in Problem 1. The oscillating steady state indeed lies on the middle branch of the -nullcline as suggested by Problem 2. In contrast, for mV:
we see that the -nullcline has shifted rightward so that only one (stable) steady state now exists, which in particular does not lie on the middle branch. The trajectory hence relaxes to rest along the -manifold, producing phasic firing.
We note also that phasic mode may be induced, for example, by modifying the conductances , , and , as shown in Problem 2.