CONTENTS Passive Properties
Nernst equation: Equilibrium and reversal potential Goldman equation: Goldman-Hodgkin-Katz constant field equation Pumps: Na/K exchange pump Permeability vs Conductance
Basic RC properties Input resistance Membrane time constant Low-pass filter
Attenuation and delay: Basic cable properties Space (length) constant: : Attenuation with distance Time constant: Time delay with distance Input resistance Non-uniform diameter: A compartmental model Spatial summation
Units: Passive property units Action Potentials (Spikes) Spike Conceptual Summary
Threshold: All-or-none? Sodium Dependence Hyperkalaemia: Too much extracellular potassium Spike Currents and Conductance Molecular Events: Single-channel patch recording Strength-Duration Curve Refractory Period Threshold Accommodation Rebound Excitation: Post-inhibitory rebound Pacemaker: Continuous regular spiking
Basic Conduction Features Conduction Velocity Collision Spatial Distribution: Voltage profile along the axon Passive vs Active Conduction: A qualitative comparison Refractory Conduction Myelination and Saltatory Conduction
Theory Classic HH: How to get a Nobel prize Sodium IV Curve Potassium IV Curve Capacitive Current: Importance in voltage clamp mechanism Patch Recording: Single channel recording under voltage clamp
Endogenous Bursters Plateau Potentials A Current: Fast (transient) potassium current
Gate Voltage Dependency Channel Voltage Dependency Putting it Together: Finding the membrane potential Non-Ohmic Channels: The GHK current equation Reduced Models: Morris-Lecar and others Phase Plane Analysis
Synapses
Excitation and Inhibition Voltage Clamp: Applied to synaptic potentials Summation and Facilitation Depolarizing and Silent IPSPs Non-Spiking Synapses Voltage-Dependent Synapses: NMDA-type Hebbian Synapses: The basis of memory? Quantal Release
Frequency Matters Size Matters Rectification
Proximal and Distal Inhibtion Directional Sensitivity
Networks
Flip-Flop Circuit Reciprocal Inhibition Oscillator Multi Phase Rhythms Photoinactivation and Optogenetics Sensory Feedback Phase Resetting Tests Tadpole Swimming: A case study Synchronization and Entrainment
Stochastic Resonance: Noise matters Lateral Inhibition
Primary Afferent Depolarization (PAD)
Phase Ambiguity
Classical Conditioning: Pavlov's dog Associative Learning: Pattern completion Facilities for Memory Models
Theory: Population-level coding in the Wilson-Cowan model Neurosim Implementation Stable States Oscillations Multi-Unit Models
Single-Channel Kinetics
Sojourn Duration and Transition Rate Constants Sojourn Distribution is Exponential
Open-Time Distribution Closed-Time Distribution is Multi-Exponential
What's the Point?: A reminder
Open-Time Distribution Closed-Time Distribution
Acetylcholine Receptor Model Auto-Correlation Compartmental Models Active Currents Neurosim Implementation Implementation
Passive Properties Noise Spiking Properties Voltage-Dependent Channels Ca-Dependent Channels Sodium Pump
Spiking Chemical Synapses Non-Spiking Chemical Synapses Carrier Ions: Calcium and Sodium Electrical Synapses Drugs: Blockers and Activators
Exponential Euler Integration Integration Components Adaptive Step Size Wilson-Cowan integration
References
$\setCounter{30}$
Contents

Copyright © W. J. Heitler (2019)

Advanced Kinetics

This section describes in detail the kinetic equations that underlie the Hodgkin-Huxley model, and some of the advanced techniques used to analyse them. It can be safely IGNORED if this is beyond what is needed by the user.

As described previously, in the classic HH model each sodium channel has 3 activation (m) gates and one inactivation (h) gate, while the potassium channel has 4 activation (n) gates. In order for a channel to be open, all the gates within it must be open. Furthermore, all the gates are independentThe assumption of gate independence is probably one of the main simplifications of the HH formalism. More sophisticated model systems can be built using Markov chain rules, in which the possible future states depend upon the current state (e.g. the inactivation gate can only close if the activation gates are open). However, most published neural models use the HH formalism described here. of each other. In the HH formalism, the individual gates act like a first-order chemical reaction with 2 states:

\begin{equation} \ce{shut <-->[\alpha][\beta] open} \label{eq:eqAlphaBeta} \end{equation}

The factors \(\alpha\) and \(\beta\) are the transition rate constants. \(\alpha\) is the number of times per second that a gate which is in the shut state opens, while \(\beta\) is the number of times per second that a gate which is in the open state shuts.

The key factor in the HH model which allows action potentials to be generated is that \(\alpha\) and \(\beta\) are voltage dependent. The equations defining the voltage-dependency for \(\alpha\) and \(\beta\) for the m, h and n gates are given in the table below.

Voltage-Dependency of Alpha and Beta
Na Activation (m) Na Inactivation (h) K Activation (n)
\[{\large \alpha_{m}} = \frac{\small{-0.1(V_{m}+45)}}{ e^ {\Large {\frac{V_{m}+45}{-10}}} \small{-1}}\] \[{\large \alpha_{h}} = 0.07 e ^ {\Large {\frac{V_{m}+70}{-20}}} \] \[{\large \alpha_{n}} = \frac{\small{-0.01(V_{m}+60)}}{ e^ {\Large {\frac{V_{m}+60}{-10}}} \small{-1}}\]
\[ {\large\beta_{m}} = 4 e ^ {\Large {\frac{V_{m}+70}{-18}}} \] \[{\large \beta_{h}} = \frac{1}{ 1 + e^ {\Large {\frac{V_{m}+40}{-10}}} }\] \[{\large \beta_{n}} = 0.125 e ^ {\Large {\frac{V_{m}+70}{-80}}} \]

As can be seen, the only variable on the right-hand side of each equation is the membrane potential Vm. This means that \(\alpha\) and \(\beta\) are directly and instantaneously dependent on the membrane potential.

It is worth making three points about the equations:

So how does the open probability of a gate depend upon \(\alpha\) and \(\beta\)? For a whole population of gates, a proportion P are in the open state, where P varies between 0 and 1. This means that a proportion 1 - P will be in the closed state. The fraction of the total population which open in a given time is dependent on the proportion of gates which are shut, and the rate at which shut gates open:

\begin{equation} \textsf{fraction of gates opening} =\alpha (1-P) \label{eq:eqFracGatesOpening} \end{equation}

and similarly

\begin{equation} \textsf{fraction of gates closing} =\beta P \label{eq:eqFracGatesClosing} \end{equation}

From this we can construct a differential equation. The rate of change of P will be the difference between the fraction opening per unit time and the fraction closing:

\begin{equation} \frac{dP}{dt} =\alpha (1-P) - \beta P \label{eq:eqDerivAlphaBetaP} \end{equation}

If \(\alpha\) and \(\beta\) are fixed, which they are if the voltage is stable, then this equation will naturally tend towards equilibrium (no change in P; dP/dt = 0). This is because if dP/dt is positive, then the rate of opening exceeds the rate of closing and P increases, which increases the rate of closing (βP) and decreases the rate of opening (α(1-P) so dP/dt becomes less positive. The same argument (in reverse) applies if dP/dt starts negative. So dP/dt ends up at zero.

When the system reaches equilibrium, the fraction of gates opening must equal the fraction of gates closing in any given period of time:

\begin{equation*} \alpha (1-P) = \beta P \end{equation*}

which rearranges as

\begin{equation} P_{\infty} = \frac{\alpha}{\alpha + \beta} \label{eq:eqGateInfinity} \end{equation}

This means that if \(\alpha\) is high and \(\beta\) is low, the gate has a high probability of being open, and vice versa. (The infinity subscript is used for P because the system only achieves equilibrium if \(\alpha\) and \(\beta\) remain stable for a relatively long period of time.)

If the membrane potential changes, and consequently the values of \(\alpha\) and \(\beta\) change, then the open probability P will also change. However, this does not occur instantly, even though the \(\alpha\) and \(\beta\) changes are instant. The rate at which P achieves its new value from a starting value of P0 is equal to the difference in the rate of shutting and the rate of opening given in the differential equation \eqref{eq:eqDerivAlphaBetaP}. This has the explicit solution in the form of a bounded exponential curve:

\begin{equation} P_{t} = P_{\infty} - (P_{\infty}-P_{0})e^{-t/\tau} \label{eq:eqGateAlphaBetaP} \end{equation}

where

\begin{equation} \tau = \frac{1}{\alpha + \beta} \label{eq:eqGateTimeConst} \end{equation}

Thus, following a change in voltage, the rate of change of P, as well as the direction and size of change, is dependent on the values of \(\alpha\) and \(\beta\). Depending on these values, some classes of gates will respond more rapidly to changes in voltage than others.

These equations can be understood as follows. We start with assuming that the system has been at a fixed constant voltage for a long period of time, and therefore P is at a starting equilibrium value \(P_{\infty}\) defined in equation \eqref{eq:eqGateInfinity}. The voltage is then changed suddenly, and \(\alpha\) and \(\beta\) immediately switch to new values appropriate to the new voltage, as defined by the equations in the table above. P then starts to change, and approaches its new equilibrium value \(P_{\infty}\) (also defined in equation \eqref{eq:eqGateInfinity}, but with the new values for \(\alpha\) and \(\beta\)) with an exponential time course with a time constant of \(\tau\) from equation \eqref{eq:eqGateTimeConst}. If either \(\alpha\) or \(\beta\) are large, then the time constant is short and P arrives at its new value rapidly. If both are small, then the time constant is long and it takes longer for P to reach equilibrium.

By combining equations \eqref{eq:eqGateInfinity} and \eqref{eq:eqGateTimeConst} we can express \(\alpha\) and \(\beta\) in terms of \(P_{\infty}\) and \(\tau\) :

\begin{equation} \alpha = \frac{P_{\infty}}{\tau} \label{eq:eqAlphaAsPInf} \end{equation}

and

\begin{equation} \beta = \frac{1-P_{\infty}}{\tau} \label{eq:eqBetaAsPInf} \end{equation}

and the differential equation \eqref{eq:eqDerivAlphaBetaP} as:

\begin{equation} \frac{dP}{dt} = \frac{P_{\infty} - P}{\tau} \label{eq:eqDerivInfTauP} \end{equation}

There is thus a simple relationship between \(\alpha\) and \(\beta\), and the equilibrium value of P and the time constant \(\tau\) with which P attains this equilibrium value.

In the equations above, P is the open-probability of an unspecified gate type. However, in the HH model, we have the specific m, h and n gates. The same equations apply, just substitute the appropriate letter identifier for P.

Take-home message: To specify the voltage-dependent characteristics of a gate type, we can

The two forms are mathematically completely interchangeable. The \(P_{\infty}\) and \(\tau\) form is usually easier to interpret in terms of qualitative function, but the \(\alpha\) and \(\beta\) form more closely reflects the underlying chemical concepts in the HH model.

Gate Voltage-Dependency

The Results display shows a voltage-clamp experiment, with the clamp potential in the top trace. However, the second trace is not current (as it would be in a real experiment), but rather the values of the open probability of the sodium and potassium activation gates (m, orange trace, and n, blue trace, respectively). The sodium inactivation gate (h) is not shown to avoid confusion with too many lines crossing over. [Remember, these probabilities could never be measured directly in a real experiment - they are just part of the model, but in terms of understanding how the model works it is useful to display their values.]

Each time there is a step change in the voltage, the m and n values change to new steady-state values. However, the change is not instant, but rather occurs with what looks like an exponential time course. This is particularly obvious for the n gate curve (blue trace) in the first clamp step.

So far so good - the model matches the theory given above.

Steady-State Probability

The HH Kinetics dialog box opens with the Probability (\(t_{\infty}\)) display option selected.

The graph shows the steady-state (time infinity) values of m, h and n across a range of physiological membrane potentials. The sodium and potassium activation gates (m, orange and n, blue respectively) both increase their open probability with depolarization, and the sodium inactivation gate (green) decreases its open probability with depolarization, as they should.

Take-home message: m, h and n gates have steady-state open probabilities whose voltage-dependency reflects their roles in sodium channel activation, sodium channel inactivation and potassium channel activation respectively.

Question: The graph shows that the n line (blue potassium activation gate) is above the m line (orange sodium activation gate) at voltages below about -30 mV, but below it at more depolarized voltages. Does this qualitatively match-up with the steady-state m and n values shown in the Results view? The clamp potentials are -70, -60 and -10 mV.

Time Constant of Probability Change

Note that m (orange line) is low down on the graph and thus has a short time constant at all membrane potentials (less than 0.5 ms). In contrast, both h and n have relatively long time constants (6 - 8 ms) around the resting potential (left-hand part of the graph), but this decreases significantly with depolarization (right-hand part of the graph).

Compare the m and n responses in the Results view. m (the orange trace) changes to its new value almost instantlyIn some published HH-type models m is assumed to have a fixed 0 or very short time constant, which can speed up computation by avoiding calculating the kinetic equations for m, without sacrificing too much accuracy., and this is true for both the rising and falling phases of the small depolarization (the first step) and the large depolarization (the second step). In contrast, n (the blue trace) responds quite slowly to both rising and falling phases of the small depolarization. This is because the membrane potential remains within the range in which the time constant is long. But it responds much more quickly to the rising phase of the large depolarization because the membrane potential has stepped into the region where its time constant is short. However, when the membrane potential steps back down to rest (the falling phase of the large step), the rate of change in n is again slow (much slower than it was on the rising phase), because the membrane potential is now once again in the long-time-constant region.

If you put this together, it means the m gates will always lead the other two types of gate. A small but significant depolarization from rest will open m gates rapidly, before the h and n gates respond. This allows the sodium channels to “escape”, and generate the rising phase of the spike. However, once the spike reaches its peak potential the other gates can “catch up”; the potassium n gates will open and sodium h gates will shut, thus terminating the spike and repolarizing the membrane. After repolarization the m gates close (de-activate) almost instantly, and thus are ready to open again if there were another depolarization. However, it takes much longer for the n gates to close (de-activate) and the h gates to open (de-inactivate), so there is a period when the membrane is less excitable than normal – the refractory period.

Take-home message: The long time constants of h and n relative to m at the resting potential allow relatively small depolarization to produce a “window of opportunity” in which there is increased sodium conductance relative to potassium conductance, and it is this that allows the spike to occur. Their long time constants at the resting potential compared to the peak potential also explain why the duration of the refractory period is long compared to that of the spike.

For completeness let us look at the voltage-dependency of \(\alpha\) and \(\beta\):

This is exactly the same file as previously, except that the Results view shows the alpha and beta values for the n gates (the potassium activation gates) in the bottom axis.

First note that the waveform transitions are abrupt (they look like square pulses). This fits with the fact that \(\alpha\) and \(\beta\), unlike n, are instantaneous functions of membrane potential. Second note that \(\alpha\) and \(\beta\) are quite small during the first clamp pulse, but that \(\alpha\) is much bigger during the second pulse. This fits with the faster time constant of n during this pulse (see \eqref{eq:eqGateTimeConst}, which shows that the time constant is short if either \(\alpha\) and \(\beta\) is large).

At the left of the graph in the hyperpolarized region the beta n line is above the alpha n line, showing that open n gates close more rapidly than closed n gates open which is why the n gates tend to be shut in this voltage region. On the right in the depolarized region, the opposite is true. However, the opening rate in the depolarized region is much faster (higher line) than the closing rate in the hyperpolarized region. This is why the time constant is faster at depolarized than hyperpolarized voltages.

As can be seen, the quantitative implications of the alpha/beta curves completely match those of the \(P_{\infty}\) and \(\tau\) curves (as indeed they should, since they are mathematically equivalent). However, in my opinion the \(P_{\infty}\) and \(\tau\) curves are more intuitive and easier to interpret in a qualitative manner.

Alpha and Beta During a Spike

The transition rate constants (\(\alpha\) and \(\beta\)) lie at the heart of the HH model, but it is impossible (at least for me) to go intuitively from a knowledge of the voltage-dependency of their values, to an understanding of how this leads to the generation of the spike. However, visualizing their values during a spike can be of at least some help in achieving a level of post hoc understanding.

The bottom axis in the Results display shows the values of the opening (\(\alpha\)) and closing (\(\beta\)) transition rate constants for the m gate (sodium activation) during a standard spike. The \(\alpha\) value actually looks a lot like the spike waveform itself, while the \(\beta\) value looks like an inverted version of the spike, but is much more distorted.

Note in the Kinetics dialog that alpha m (curving upwards to the right) shows a more-or-less linear voltage-dependency above about -45 mV (the cross-over point of the two curves). This explains why the waveform of (\(\alpha\)) in the Results view follows the shape of the spike quite closely. The voltage-dependency of beta m, on the other hand, is in the opposite direction and is much more non-linear. This is why the waveform of \(\beta\) is like an inverted spike, but more distorted. Since \(\alpha\) reflects the rate at which shut activation gates open and \(\beta\) reflects the rate at which they shut, it makes sense that the former follows the waveform of the spike, while the latter follows an inverted version of that waveform.

Neurosim does not directly display the time constants for the gate probability values during a simulation run, but given the values of \(\alpha\) and \(\beta\) these can easily be calculated from equation \eqref{eq:eqGateTimeConst}.

Optional task: Plot a graph of the spike waveform against time, and include the time constants for the gate activation/inactivation variables m, h and n on the same graph.

You should see that the sodium activation variable (m) has a relatively fast time constant at all times in the simulation, but that sodium inactivation and potassium activation have short(ish) time constants around the peak of the spike, but long time constants before and after the spike. This fits with the discussion of the time constants under voltage clamp conditions given earlier.

Channel Voltage-Dependency

Potassium Channels

So far we have looked at how the gates within a channel respond to voltage. How about the channel as a whole?

We can now see the potassium conductance (middle axis) as well as the probability of an n gate being open (bottom axis). In the HH model there are 4 n gates in each potassium channel, and they all have to be open for the channel to be open. Therefore, if the maximum potassium conductance that would be achieved if all the gates in all the channels were open is gK max, then the actual conductance gK is:

\begin{equation} g_{{\small K}} = n^{4} g_{{\small K \, max}} \label{eq:eqN4} \end{equation}

Let us check that this all works out.

When you look at the numbers you will see, perhaps surprisingly, the value of n is quite large (0.32) even at the resting potential of -70 mV. This means that about 1/3 of the n gates are open even at rest. However, the overall potassium conductance is still small, because 0.324 is only 0.01. The maximum potassium conductance (all gates and all channels open) is 36 mS.cm-2, as can be seen in the Membrane Properties group in the Setup view, so the actual conductance of the voltage-dependent potassium channels at rest is only about 0.3 mS.cm-2. However, this is about the same as the leakage conductance, which means that the 1% of voltage-dependent K channels that are open at rest contribute about as much to the resting potential as do the leakage channels themselves.

Take-home message: Even at the resting potential, open voltage-dependent potassium channels make a significant contribution to determining the membrane potential.

Sodium Channels

We are now looking at sodium conductance (middle axis) and the m and h probabilities (bottom axis). The sodium conductance is calculated from an equation similar to \eqref{eq:eqN4}, but adjusted for the different gate types:

\begin{equation} g_{{\small Na}} = m^{3}h g_{{\small Na\, max}} \label{eq:eqM3H} \end{equation}

Some key features are worth noting in the Results.

At the resting potential (-70 mV) on the left of the display, the value of h (green trace) is about 0.6. This means that even at rest, about 40% of sodium channels are inactivated. These channels are unavailable for recruitment by depolarization: even if all their m gates open, they will still not switch into the channel-open condition. However, these channels can be made available by initially hyperpolarizing the membrane before depolarizing it, so that the value of h increases above 0.6. The slow time constant of h in this voltage range means that its value will remain elevated for some time after release from hyperpolarization, and depolarization occurring in this time will be more effective in opening sodium channels. This is one of the mechanisms underlying rebound excitation demonstrated previously.

There is a narrow spike of increased sodium conductance at the start of the large depolarizing step. This reflects the “window of opportunity” mentioned earlier. The m gates open almost instantly because of their short time constants, and the h gates close rapidly, but not so rapidly as the m gates open. There is therefore a brief period of increased sodium conductance.

The peak conductance is only about 26 mS.cm-2, out of an available maximum of 120 mS.cm-2 (see the Membrane Properties in the Setup view). So this depolarization is only opening about 1/5 of the sodium channels present in the membrane.

Finally, note the slow recovery of the h value after the large depolarizing step. This means that it will be some time before the normal complement of sodium channels are available for spike generation, and this, along with the elevated value of n (potassium channel activation gate) observed earlier, underlies the refractory period.

 

Putting it all together: finding the membrane potential

Solving the HH equation (i.e. calculating the value of the membrane potential as time passes) requires solving 4 differential equations simultaneously. Three equations are like equations \eqref{eq:eqDerivAlphaBetaP} or \eqref{eq:eqDerivInfTauP} and use the instantaneously voltage-dependent factors to describe the rate of change in the gating variables:

\begin{align} \label{eq:eqGrandHH1} \frac{dm}{dt} &= \alpha_{m(V)}(1-m)-\beta_{m(V)} &= \frac{m_{\infty(V)} - m}{\tau_{m(V)}} \\[1.5ex] \frac{dh}{dt} &= \alpha_{h(V)}(1-h)-\beta_{h(V)} &= \frac{h_{\infty(V)} - h}{\tau_{h(V)}} \\[1.5ex] \frac{dn}{dt} &= \alpha_{n(V)}(1-n)-\beta_{n(V)} &= \frac{n_{\infty(V)} - n}{\tau_{n(V)}} \end{align}

Integrating these equations yields the moment-by-moment values for m, h and n in the HH model. Once these values have been calculated, they are put into equations \eqref{eq:eqN4} and \eqref{eq:eqM3H} to calculate the potassium and sodium channel conductances. Once the conductances are known, the driving force equation is used to calculate the current flowing through each ion channel type (sodium, potassium and the fixed-conductance leakage channels). So the end point of this is that we know the currents through the ion channels, and we know the stimulus current because we control that.

The fourth equation depends on Kirchoff’s rule, which states that the sum of currents entering and leaving any point in a circuit is zero. Since the only pathsAssuming the cell is isopotential and that there is no internal current flow to a different region such as an axon or dendrites. If there is such lateral flow, you need to build a compartmental model. into or out of the cell are through the membrane capacitance or the membrane conductance, their currents must be equal and opposite in order to achieve a total of zeroBut see compartmental models for a slightly more complex situation.. The current charging a capacitor is directly proportional to the rate of change of voltage over the capacitor, as given by equation (I = C dV/dt), so the rate of change of voltage multiplied by the constant of proportionality (the capacitance) must be equal to the total current through the ion channels, plus any stimulus current. This leads directly to the fourth equation:

\begin{equation} \begin{aligned} C\frac{dV}{dt} = I_{stim} &- g_{Na\,max}m^{3}h(V-E_{Na}) \\ &-g_{K\,max}n^{4} (V-E_{K}) \\[1.5ex] &-g_{leak}(V-E_{leak}) \end{aligned} \label{eq:eqGrandHH} \end{equation}

Integrating this equation is what actually calculates the membrane potential V at each moment in time.

Capacitive current (spike)

The bottom trace is the capacitive current that actually causes the changes in voltage. It is the inverted sum of the sodium, potassium and leakage currents minus the stimulus current (the minus is because the stimulus current flows in the opposite direction to the channel currents).

Note that there is a large positive (upward) spike in capacitive current during the rising phase of the action potential.

Question: What do you think might be special about the spike waveform at this moment in time?

Question: Why is the capacitive current zero at the peak of the action potential?

On the falling phase of the action potential there is a more sustained negative plateau of current. This is smaller than the positive spike because in an action potential, the rate of fall of voltage is slower than the rate of rise.

Note that there is a small negative trough in the capacitive current late in the falling phase of the action potential (at about 4.1 ms into the experiment). This may not have any biological significance, but it is an emergent feature of the HH model.

Question: What feature in the waveform corresponds to this peak?

Task: Check that the capacitive current really is proportional to the derivative of the membrane potential, as claimed in \eqref{eq:eqGrandHH}.

Use the Copy text drop-down option of the Copy button to copy the text values of the traces to the clipboard.  Start Excel (or a similar program) and paste in the text values. Construct a formula to calculate the negative time derivative of the membrane potential. A simple adjacent-difference method is fine (see below). Confirm that the derivative is equal to the capacitive current at each time point throughout the experiment by constructing an X-Y scatterplot, with time as the X value and the derivative and total current as two data series on the Y axis. By happy coincidence, with the default value for capacitance, the units should end up the same and the lines should superimpose. (There might be slight variation because the derivative calculation is rather crude, but the values should be very close.)

Excel derivative calculation
Excel Derivative. Calculate adjacent difference derivative in Excel.

Non-Ohmic Channels: The GHK Current Equation

A key simplifying assumption of the HH equations, and one that is followed in most published research models, is that the current flowing through ion channels can be described by a modified version of Ohms Law - the driving force equation (\(I=g(V_{m} - V_{eq})\)), and that the core channel conductance when all gates within it are open has a fixed value (gmax). This implies that any changes in channel conductance are due to the opening and closing of the voltage-dependent gates within it.

However, this is usually not exactly true. If there is an imbalance in the concentration of the carrier ion on either side of the membrane, then the core channel conductance will almost certainly rectify – it will have a higher conductance in the “downhill” direction of the concentration gradient than in the “uphill” direction. This follows from basic physical chemistry (there are more ions available to carry current from high to low concentration than vice versa), and Hodgkin and Huxley were well aware of it. They were able to ignore it because the effect is quite small in the voltage and concentration ranges with which they were dealing. However, if the concentration ratio is large, the effect can become significant, and it may be better to use a more sophisticated model of channel conductance based on permeability. This involves the Goldman-Hodgkin-Katz (GHK) currentThe GHK current equation underlies the GHK voltage equation, which we have come across before in the calculation of the resting potential. equation, which is an equation based on electro-diffusion theory that calculates the current flow across a membrane given the voltage, permeability and ion concentrations:

\begin{equation} i = PV \frac{z^2F^2}{RT} \left(\frac{[in]\mu - [out]}{\mu - 1} \right)\\[3ex] \textsf{where } \mu = e^{{\Large\frac{VzF}{RT}}} \label{eq:eqGHKCurrent} \end{equation}

In this equation i is the current density carried by the ion in question, [in] and [out] are the intracellular and extracellular concentrations of the ion respectively, P is the membrane permeabilityPermeability, unlike conductance, is a fixed value that is independent of ion concentration. of that ion, V is the membrane potential, z is the charge (valency) of the ion, R is the gas constant (8.314 J K-1 mol-1), F is Faraday’s number (96485 C mol-1), and T is the absolute temperature (294°K).

In Neurosim, the equation assumes that the user-specified value of P is relative to 1 cm2 of membrane, and during the calculation it is scaled by neuron surface area to produce an absolute current. Calculating the GHK equation is considerably more time consuming than calculating the standard driving force equation, which is why computational models avoid using it if possible.

Rectification

It is clear from the equation \eqref{eq:eqGHKCurrent} that for a constant voltage, the trans-membrane current is dependent on the absolute ion concentrations on either side of the membrane, and the difference between those concentrations. This means that the channel does not obey Ohm’s Law.

The initial result shows the membrane potential change induced by a negative current pulse, and at this stage it looks much like the response of a standard RC circuit.

First, note that despite the name of the dialog, this channel is not voltage dependent. It has a singleNeurosim requires that identified channels have at least one gate in them. This also means that if required the channel can have conventional gate-derived voltage dependence, as well as the voltage dependency conferred through using the GHK current equation. gate in it, but that gate is fixed in the open state (the time-infinity value of the activation variable is 1, and the time constant is 0). So we are observing the core properties of the channel, unaffected by any voltage-dependent gates.

Next, note that the Use GHK current option is selected (towards the bottom-left of the dialog). The carrier ion is unspecified, but it is positively charged and monovalent (z = +1), and the intracellular concentration is much higher than the extracellular concentration.

As you would expect, the smaller stimulus produces a smaller voltage deflection.

It should be clear that the largest positive current stimulus produces a smaller voltage deflection than the equivalent-sized negative current stimulus. Remember, for a given current, a smaller voltage response implies a larger conductance, so the result shows that it is “easier” for a positive current to push the positively-charged carrier ion out of the cell (down its concentration gradient), than for a same-sized negative current to pull it into the cell (up its concentration gradient). This is rectification.

Task: Construct a graph showing the relation between stimulus current and voltage response (a classic I/V curve).

If the response was Ohmic as in a classic RC circuit, the plot should be a straight line, but in fact it curves upwards on the right. A curved I/V plot is a diagnostic feature indicating rectification.

GHK and Leakage Currents

In this experiment, the current is passing through two channels – the GHK channel, and the standard leakage channels, the latter of which in Neurosim are Ohmic in their response. To separate these effects, we need to switch to voltage clamp mode.

The upper axis in the Results view shows the clamp potential, with a step change from a holding potential of -90 mV to a clamp potential of -140 mV. The lower axis shows two current tracesIn a real experiment, separating these two currents would probably require some drug that could specifically block the GHK channel, and then doing a subtraction experiment measuring the total current with and without the drug. But the beauty of simulation is that we can just ask the computer to generate separate traces for the two currents!. The upper trace (red) is the current through the GHK channel, the lower trace (green) shows the current through the leakage channel.

Note that there is a rather small change in the GHK current, but a bigger change in the leakage current.

As the voltage increases, the current step in the GHK channel (red trace) also increases. This again indicates rectification in the GHK channel - equal-sized voltage changes would produce equal-sized current changes in an Ohmic conductance. The latter can be seen in the leakage current (green trace), where each current step is of equal size.

Task: As before, construct an I/V plot of the voltage versus current changes, but this time switch the Y axis selection between the leakage current (I Leak) and the GHK channel current (I Chan 1), and compare the two plots.

The leakage I/V plot should be linear, crossing the 0 current intercept at -60 mV. This is consistent with the leakage equilibrium potential of -60 mV (you can open the Passive Properties dialog to check the settings if you wish). The GHK I/V plot is curved, indicating rectification. It does not cross the current intercept because the GHK equilibrium potential is more negative than any voltage in the experiment. (It is actually -192.6 mV, as can be seen in the Channel Properties dialog which uses the ion concentration ratio to calculate the Nernst potential.)

Take-home message: Real ion channels are often rectifying – they conduct current more easily in one direction than the other. This arises naturally if there is a difference in the concentration of the carrier ion on either side of the membrane. However, the effect is usually small, and can be ignored in most simulations.

The GHK current equation provides a more sophisticated description of channel behaviour than the standard Ohmic model. However, it too has limitations. An assumption underlying the equation is that ions move completely independently of each other, whereas in fact many channels act like narrow tunnels in which ions cannot move past each other. It also assumes that the membrane is homogeneous, with a uniform electric field and a uniform ion concentration gradient within the membrane, which is probably not true for proteinaceous channel pores. So although the GHK equation might be (somewhat) more accurate than the standard driving force equation, it is certainly not the last word on the matter!

Conductance

In a GHK channel the permeability is specified, but not the conductance. Neurosim calculates and displays an effective conductance by “reverse engineering” the driving force equation. This yields what is known as the chord conductance (see Thompson, 1986 for a detailed discussion of chord and slope conductances), which is analogous to the conductance used in the HH equations. However, the conductance value can be inaccurate if the membrane potential is close to the equilibrium potential because then the calculation involves dividing a very small current by a very small driving force. If the membrane potential is actually at the equilibrium potential, the conductance cannot be calculated due to a divide-by-zero singularity (avoided in the simulation by offsetting the membrane potential slightly).

 

Morris-Lecar: Reduced Kinetics

So far we have been looking at the Hodgkin-Huxley model and various complex elaborations developed from it. However, a very influential model has been one that reduced the complexity - the Morris-Lecar model (Morris and Lecar, 1981).

The simulation implements the model with the kinetic parameters specified in the original paper, which were based on muscle fibres in a barnacle. A stimulating current is applied to reproduce the behaviour shown in Figure 6 of that paper.

The spikes do not look very “convincing” in terms of what we expect nerve spikes to look like, but the important point about this model is that it produces spike-like behaviour from an extremely simple basis. There are only two channels, one with a depolarized equilibrium potential, nominally a calcium channel, and the other with a hyperpolarized equilibrium potential, nominally a potassium channel. Each channel has just a single gate, and neither channel inactivates. The spike-like behaviour arises purely from the activation kinetics of these two channels.

The calcium activation curve has a mid point at 0 mV (note the v – 0 component in the user defined equation – the 0 is obviously unnecessary here, but it indicates where the mid-point parameter would be placed if it were non-zero), which is very depolarized. But so is the resting potential of the cell (-50 mV), and it takes a depolarizing stimulus to activate spikes.

The activation curve moves slightly to the right (note v – 10 in the equation). This means that the potassium channels open at an even more depolarized membrane potential than the calcium channels – they will open after them as the membrane depolarizes in response to the stimulus.

Note from the graph Y-scale value that the time constant peaks at 10 ms. Also note the factor 0.1 which multiplies the cosh in the equation.

Now the maximum time constant is only 1 ms (the 0.1 cosh factor in the equation has disappeared). So not only does the calcium channel activate at a lower voltage than the potassium channel, but also it activates much more rapidly. This is an example of a so-called fast-slow model, where there are just 2 simple channel types, but with a significant difference in the time-scales at which they operate.

Equations

The user-defined equations involve terms which may be unfamiliar – the hyperbolic tangent and cosine functions (tanh, cosh). However, the equations could be re-written in exponential form to give similar results.

The time-infinity activation curve does not change at all with the parameters chosen for the pre-defined equation. You can switch back and forth between the Pre-defined and User-defined equations with no change.

Further reduction

Solving the Morris-Lecar model requires integrating 3 differential equations: one for m, one for n, and one for V. This is in contrast to the full HH model, which requires solving 4 differential equations (\eqref{eq:eqGrandHH1} - \eqref{eq:eqGrandHH}) (the additional one is for h). However, the calcium channel gate time constant is very short. What happens if it is removed altogether?

Note that the system still responds to the stimulus with a train of action potentials. The frequency and shape of the spikes is somewhat different, but essentially the system behaviour is unchanged when the short time constant is eliminated completely.

By removing the time-dependence of m we have made it an instantaneous function of voltage, and its value no longer has to be calculated by integrating a differential equation. It can be calculated explicitly just by solving the Infinity value equation at the required voltage. This means that the reduced Morris-Lecar model only requires solving two differential equations.

Significance

The significance of the reduced Morris-Lecar model lies mainly in its mathematical tractability as a dynamical system with just two differential equations. Such systems are amenable to a technique called phase-plane analysis, and there have been numerous papers published which use this to explore the properties of this system. This is discussed further in the next section.

Take-home message: The Morris-Lecar model uses the HH formalism but with fewer parameters. It produces spikes with just two single-gate non-inactivating voltage-dependent channel types, a fast depolarizing channel and a slow hyperpolarizing channel. The reduced form of the model involves only 2 differential equations, and is highly amenable to mathematical analysis.

 

Phase Plane Analysis

The original Hodgkin-Huxley model is a dynamical system described by a set of four coupled differential equations with state variables m, h, n and v \eqref{eq:eqGrandHH1} - \eqref{eq:eqGrandHH}. These variables change with time, and normal analysis and display involves plotting them (or dependent variables such as conductance or current) against time itself. That is what most of the displays in Neurosim show. However, if the variables are plotted against each other rather than against time, the result is a plot in phase space.

A plot with 4 dimensions is virtually impossible to display graphically, and equally difficult to analyse. However, if the system can be described by differential equations involving just two variables, then it is easy to make a phase plane plot on a 2D surface like a computer monitor (or, in the old days, a sheet of paper ...).

Why would one want to do that?

First, there is a very well developed set of mathematical procedures that can analyse the properties of just two coupled differential equations. This analysis can tell whether the system is stable, multi-stable, prone to regular oscillation, chaotic, explosively unstable etc. The other advantage of a phase plane plot is simply that it gives an alternative view of the shape of the system response, which can emphasise features that are not immediately obvious in a normal time-series plot.

The mathematical analysis of phase plots is beyond the remit of this tutorial (and the expertise of this author), but Neurosim can produce phase plane plots that give a qualitative flavour of the methodology.

The key point about this model is that it is driven by just two differential equations. There is an equation that specifies the voltage that operates on a fast time scale, and another that specifies the open probability of the potassium channel that operates on a slow time scale.

As the simulation progresses, you should see the plot trace a clockwise trajectory around the graph, with each spike generating an individual loop. (If you have selected too large a slow-down factor and you get bored, you can reduce it on-the-fly during the simulation.) While the simulation runs, the plot is black-and-white, but when it terminates it is colour-coded by time. This can be seen more clearly in a 3-D graph.

As the plot rotates, you should see that the plot trajectory forms a helix, with the blue end being the start of the simulation (time 0) and the yellow end being the end (time 800 ms).

So what does all this tell us? The first and probably most important point is that the repeating loops in the plot form what is called a limit cycle, and limit cycles are characteristic of oscillating dynamical systems. The limit cycle can be understood by following through the various stages as follows.

To further examine this, what happens if there are no spikes, or if spikes fail?

The stimulus is now below threshold so no spikes are generated. The membrane potential rises and then falls as the stimulus occurs, but the open probability of the n gates does not significantly change in this voltage range. Consequently, the phase plot is just a vertical line, with no sign of a limit cycle.

There can also be spike failure due to excessive stimulation. To see this most clearly we need to remove the noise.

The Results view shows spikes failing due to massive activation of potassium channels. The phase plot shows the limit cycle spiralling inwards to a single point, called a stable focus.

The 3-D plot trajectory has a conical shape that tapers to a straight line in the Z (time) axis, reflecting the focus in the 2-D plot

Take-home message: Dynamical systems with stable oscillations show a limit cycle (a circular trajectory) in the phase plane. If the oscillations are not stable but fade out, the limit cycle will spiral in towards a stable focus.

Phase Plots of Experimental Data

The main use of phase plots in computational neuroscience is for the mathematical analysis of model systems, where the analyst has access to the state variables of the underlying equations. However, in experimental neuroscience, phase plots are also quite commonly used to present an alternative visualization of a system response. In these cases, it is usual to plot the slope of the membrane potential against the potential itself (i.e. dV/dt against V).

We don't have real experimental data in Neurosim, so we will just use an HH spike instead.

A single spike shows in the Results view, and the Phase Plot shows a limit cycle similar to what we saw previously.

Task: Bearing in mind that the phase trajectory moves in a clockwise direction, work through how the phase plot arises from the voltage profile of the spike in the Results view. Here are some questions to help you focus:

Questions:

Now let us further increase the extracellular potassium concentration.

The limit cycle rapidly collapses into a stable focus. Which is pleasing visually, but would, of course, be disastrous in a biological system.

 

On to Synapses ...