In the late 1940’s and early 1950’s Alan Hodgkin and Andrew Huxley elucidated the biophysical underpinnings of nerve excitation. This tutorial walk the user through the steps necessary to reproduce and understand key aspects of their Nobel Prize-winning work using the MATLAB simulation environment. TheTutorial on CellML, OpenCOR & the Physiome Model Repository shows you how to implement the Hodgkin-Huxley model using the OpenCOR model authoring and simulation tool.

### Video Tutorial

### Step 1: Simulating the dynamics of the potassium conductance

Potassium conductivity $g_K(t)$ is simulated as proportional to the forth power of a gating variable, $n(t)$:

$$g_K=n^4 (\bar{g_K} ) $$

where $\bar{g_K}$ is the maximum conductance of potassium in the cell, occurring when $n = 1$. The dynamics of potassium conductivity are simulated to respond to changes in membrane potential via a first-order process

$$dn/dt=α_n (1-n)- β_n n$$

where $α_n$ and $β_n$ are the rates of opening and closing of the channel. The dependency of the opening and closing rates on voltage are modeled by functions that match the observed data:

$$α_n=0.01 \frac{(10-v)}{\exp\left(\frac{10-v}{10}\right)-1}$$

$$β_n=0.125 \exp \left(\frac{-v}{80} \right)$$

where $v$ is the membrane voltage (inside cell minus outside cell) and measured relative to the resting potential of around -60 mV. The units in the above expression are in mV. A comparison of the above equations for opening and closing rates to the original data derived from Hodgkin and Huxley’s experiments is illustrated below:

Estimated opening and closing rate constants for the potassium conductance as functions of membrane voltage are show on left. These data show that as the membrane depolarizes ($v$ becomes more positive) the opening rate increases and the closing rate decreases. At high values of $v$ the channels open. At low voltages the channels close. The model captures this important feature of the nerve cell conductivities through the above equations for $α_n$ and $β_n$ that effectively capture the data measured by Hodgkin and Huxley.

The time-dependent behavior of the potassium conductance can be simulated using a program that integrates the equation for $dn/dt$. Simulation in MATLAB requires a function that returns the computed time derivative $dn/dt$ as a function of $n(t)$ and $v$. The MATLAB function dXdT_n.m has the following syntax:

function [f] = dXdT_n(~,n,v)

% FUNCTION dXdT_N% Inputs: t - time (milliseconds)

% x - vector of state variables {n}

% V - applied voltage (mV)

% Outputs: f - vector of time derivatives

% {dn/dt}% alphas and betas:

a_n = 0.01*(10-v)./(exp((10-v)/10)-1);

b_n = 0.125*exp(-v/80);% Computing derivative dn/dt:f(1,:) = a_n*(1-n) - b_n*n;

This function can be integrated using MATLAB to simulate a voltage-clamp experiment. For example, we can simulate the response to a voltage-clamp experiment, with voltage set to $v = 100$ mV, and initial condition $n = 0$ with the following script:

v = 100; xo = 0; g_K = 35;[t,x] = ode23s(@dXdT_n,[0 12],xo,[],v);g = g_K*(x(:,1).^4);plot(t,g)xlabel('t (ms)');ylabel('g_{K}');set(gca,'fontsize',16)

which gives the output illustrated to the right above. (See Tutorial Enzyme Kinetics with MATLAB 2 for more on how to simulate differential equations using MATLAB.)

Hodgkin and Huxley’s full analysis of the potassium ion conductivity can be reproduced using the script HH_potassium_current.m, which includes data extracted from the original publication (A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol*., 117:500–544, 1952), and uses the function dXdT_n.m to simulate the conductivity response to a series of voltage-clamp experiments. This script generates the plot below, which may be compared to Figure 3 in Hodgkin and Huxley (*J. Physiol*., 117:500–544, 1952).

### Step 2: Simulating the dynamics of the sodium conductance

Sodium conductivity $g_{Na}(t)$ is simulated as governed by two gating variables, $m(t)$ and $h(t)$:

$$g_Na=m^3h(\bar{g_{Na}}) $$

where $\bar{g_{Na}}$ is the maximum conductance of sodium in the cell, occurring when $m = 1$ and $h = 1$. As for potassium conductivity, the dynamics of sodium conductivity are simulated to respond to changes in membrane potential via first-order processes

$$dm/dt=α_m (1-m)- β_m m$$

$$dh/dt=α_h (1-h)- β_hh$$

The difference is that since there are two gating variables there are two opening rates α_m and α_h and two closing rates $β_m$ and $β_h$. The equations that capture the voltage dependency of these variables are

$$α_m=0.1 \frac{25-v}{\exp\left( \frac{25-v}{10} \right) -1}$$

$$β_m=4 \exp\left( \frac{-v}{18} \right)$$

$$α_h=0.07 \exp \left( \frac{-v}{20} \right)$$

$$β_h=\frac{1}{\exp\left(\frac{30-v}{10}\right)+1}$$

Again, v is the membrane voltage (inside cell minus outside cell) and measured relative to the resting potential of around -60 mV. The units in the above expression are in mV. Comparison of the above equations for opening and closing rates to the original data derived from Hodgkin and Huxley’s experiments is illustrated below:

These data reveal that the m and h variables operate on substantially different time scales. Like the potassium conductance, the sodium conductance opening rate increases with increasing v. However, the rates are about 10 times higher than the rates for the potassium conductance. For this reason, the channels associated with the sodium conductance have been termed fast sodium channels.

The h gating variable shows qualitatively opposite behavior to that of the m gating variable. It varies more slowly and tends to close in response to an increase in v.

The time-dependent behavior of the sodium conductance can be simulated using a program that integrates the equations for dm⁄dt and dm⁄dt: dXdT_mh.m has the following syntax:

function [f] = dXdT_mh(~,x,v)

% FUNCTION dXdT_MH

% Inputs: t - time (milliseconds)

% x - vector of state variables {m,h}

% V - applied voltage (mV)

%

% Outputs: f - vector of time derivatives% {dm/dt,dh/dt}

% state variables

m = x(1);

h = x(2);% alphas and betas:

a_m = 0.1*(25-v)/(exp((25-v)/10)-1);

b_m = 4*exp(-v/18);

a_h = 0.07*exp(-v/20);

b_h = 1 ./ (exp((30-v)/10) + 1);% Computing derivatives:

f(1,:) = a_m*(1-m) - b_m*m;

f(2,:) = a_h*(1-h) - b_h*h;

This function can be integrated using MATLAB to simulate a voltage-clamp experiment. For example, we can simulate the response to a voltage-clamp experiment, with voltage set to $v = 100$ mV, and initial condition $n = 0$ with the following script:

v = 100; xo = [0 1]; g_Na = 35;

[t,x] = ode23s(@dXdT_mh,[0 12],xo,[],v);

g = g_Na*(x(:,1).^3).*x(:,2);

plot(t,g)

xlabel('t (ms)');

ylabel('g_{Na}');

set(gca,'fontsize',16)

which gives the output illustrated to the right. Here we use the initial condition $m(0) = 0$, $h(0) = 1$, which means that the conductivity is initially zero. Because voltage is clamped at $v = 100$ mV, the ‘fast’ conductivity associated with the m gate rapidly increases, and the slower h gate eventually closes, resulting in the biphasic transient captured by the model.

Hodgkin and Huxley’s analysis of the sodium ion conductivity can be reproduced using the script HH_sodium_current.m, which includes data extracted from the original publication (A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 117:500–544, 1952), and uses the function dXdT_mh.m to simulate the conductivity response to a series of voltage-clamp experiments. This script generates the plot below, which may be compared to Figure 6 in Hodgkin and Huxley ( J. Physiol., 117:500–544, 1952).

Note that there is a slight mis-match between data and simulation. This is because these simulations use the functions for $α_m$, $α_h$, $β_m$, and $β_h$ that are associated with the best fit to values extracted from the ensemble of axons characterized by Hodgkin and Huxley, while the data on the right represent individual recordings.

### Step 3: Putting it all together and simulating the action potential

The overall model is represented by the circuit model on the right, where $V_{Na}$ and $V_K$ are the Nernst potentials for sodium and potassium. Since voltage is measured as inside potential minus outside potential, for a typical cell with a value of $V_K$ of approximately −70 mV, the membrane potential and the Nernst potential work in opposition when the $V \approx −70$ mV. The sodium concentration, however, is typically higher on the outside of the cell, and $V_Na$ may be in the range of +50 mV. Thus when $V = −70$ mV, the thermodynamic driving force for Na$^+$ influx is $−(V − V_{Na}) \approx +120$ mV.

In the model for the action potential, membrane voltages are expressed relative to resting potential: $v = V – V_o$, where $V$ is the absolute membrane potential (inside minus outside potential) and $V_o = -56$ mV is the resting potential. The governing equation for the circuit model is

$$C_m dv/dt=-g_{Na} (v-v_{Na} )-g_K (v-v_K )-g_L (v-v_L )+I_{app}$$

where the membrane capacitance has the value $C_m = 1\times 10^{-6}$ μF cm$^{-2}$. The term $g_L (v-v_L )$ represents a leak current and $I_{app}$ is the current externally injected into the cell. The Nernst potentials for the squid giant axon prep are $v_{Na} = V_{Na} - V_o = 115$ mV, $v_K = V_K - V_o = -12$ mV, and $v_L = 10.6$ mV.

Combining the equations for the gating variables with the equation for membrane potential, we can write a function for the full combined model:

function [f,varargout] = dXdT_HH(~,x,I_app)

% FUNCTION dXdT_HH

% Inputs: t - time (milliseconds)

% x - vector of state variables {v,m,n,h}

% I_app - applied current (microA cm^{-2})

%

% Outputs: f - vector of time derivatives

% {dv/dt,dm/dt,dn/dt,dh/dt}% Resting potentials, conductivities, and capacitance:

V_Na = 115;

V_K = -12;

V_L = 10.6;

g_Na = 120;

g_K = 36;

g_L = 0.3;

C_m = 1e-6;

% State Variables:

v = x(1);

m = x(2);

n = x(3);

h = x(4);

% alphas and betas:

a_m = 0.1*(25-v)/(exp((25-v)/10)-1);

b_m = 4*exp(-v/18);

a_h = 0.07*exp(-v/20);

b_h = 1 ./ (exp((30-v)/10) + 1);

a_n = 0.01*(10-v)./(exp((10-v)/10)-1);

b_n = 0.125*exp(-v/80);

% Computing currents:

I_Na = (m^3)*h*g_Na*(v-V_Na);

I_K = (n^4)*g_K*(v-V_K);

I_L = g_L*(v-V_L);

% Computing derivatives:

f(1) = (-I_Na - I_K - I_L + I_app)/C_m;

f(2,:) = a_m*(1-m) - b_m*m;

f(3) = a_n*(1-n) - b_n*n;

f(4) = a_h*(1-h) - b_h*h;

% Outputting the conductivities

varargout{1} = [(m^3)*h*g_Na (n^4)*g_K g_L];

The model may be simulated with the script HodHux.m to generate the following output.

Simulated action potential from the Hodgkin-Huxley model.

### References

Hodgkin AL and AF Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 117:500–544, 1952

Beard DA. Biosimulation. Chapter 8: Cellular Electrophysiology. Cambridge University Press, 2012

www.cambridge.org/biosim