A model of ion channel gating and current: Introducing CellML units¶
A good example of a model based on a first order equation is the one used by Hodgkin and Huxley [AAF52] to describe the gating behaviour of an ion channel (see also next three sections). Before we describe the gating behaviour of an ion channel, however, we need to explain the concepts of the ‘Nernst potential’ and channel conductance.
An ion channel is a protein or protein complex embedded in the bilipid membrane surrounding a cell and containing a pore through which an ion \(Y^{+}\) (or \(Y^{}\)) can pass when the channel is open. If the concentration of this ion is \(\left\lbrack Y^{+} \right\rbrack_{o}\) outside the cell and \(\left\lbrack Y^{+} \right\rbrack_{i}\) inside the cell, the force driving an ion through the pore is calculated from the change in entropy.
Entropy \(S\) (\(J.K^{1}\)) is a measure of the number of microstates available to a system, as defined by Boltzmann’s equation \(S = k_{B}\text{lnW}\), where \(W\) is the number of ways of arranging a given distribution of microstates of a system and \(k_{B}\) is Boltzmann’s constant1. The driving force for ion movement is the dispersal of energy into a more probable distribution (see Fig. 12; cf. the second law of thermodynamics2).
The energy change \(\Delta q\) associated with this change of entropy \(\Delta S\) at temperature \(T\) is \(\Delta q = T\Delta S\) (J).
For a given volume of fluid the number of microstates \(W\) available to a solute (and hence the entropy of the solute) at a high concentration is less than that for a low concentration3. The energy difference driving ion movement from a high ion concentration \(\left\lbrack Y^{+} \right\rbrack_{i}\) (lower entropy) to a lower ion concentration \(\left\lbrack Y^{+} \right\rbrack_{o}\) (higher entropy) is therefore
\(\Delta q = T\Delta S = k_{B}T\left( \ln{\left\lbrack Y^{+} \right\rbrack_{o}  \ln\left\lbrack Y^{+} \right\rbrack_{i}} \right) = k_{B}T\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\) (\(J.ion^{1}\))
or
\(\Delta Q = RT\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\) (\(J.mol^{1}\)).
\(R = k_{B}N_{A} \approx 1.34x10^{23} (J.K^{1}) \text{x} 6.02x10^{23} (mol^{1}) \approx 8.4 (J.mol^{1}K^{1})\) is the ‘universal gas constant’4. At 25°C (298K), \(\text{RT} \approx 2.5 kJ.mol^{1}\).
Every positively charged ion that crosses the membrane raises the potential difference and produces an electrostatic driving force that opposes the entropic force (see Fig. 13). To move an electron of charge e (\(\approx 1.6x10^{19}\) C) through a voltage change of \(\Delta\phi\) (V) requires energy \(e\Delta\phi\) (J) and therefore the energy needed to move an ion \(Y^{+}\) of valence z=1 (the number of charges per ion) through a voltage change of \(\Delta\phi\) is \(\text{ze}\Delta\phi\) (\(J.ion^{1}\)) or \(\text{ze}N_{A}\Delta\phi\) (\(J.mol^{1}\)). Using Faraday’s constant \(F = eN_{A}\), where \(F \approx 0.96x10^{5} C.mol^{1}\), the change in energy density at the macroscopic scale is \(\text{zF}\Delta\phi\) (\(J.mol^{1}\)).
No further movement of ions takes place when the force for entropy driven ion movement exactly equals the opposing electrostatic driving force associated with charge movement:
\(\text{zF}\Delta\phi = \text{RT}\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\) (\(J.mol^{1}\)) or \(\Delta\phi = E_{Y} = \frac{\text{RT}}{\text{zF}}\ln\frac{\left\lbrack Y^{+} \right\rbrack_{o}}{\left\lbrack Y^{+} \right\rbrack_{i}}\) (\(J.C^{1}\) or V)
where \(E_{Y}\) is the ‘equilibrium’ or ‘Nernst’ potential for \(Y^{+}\). At 25°C (298K), \(\frac{\text{RT}}{F} = \frac{2.5x10^{3}\ }{0.96x10^{5}} (J.C^{1}) \approx 25mV\).
For an open channel the electrochemical current flow is driven by the open channel conductance \({\overset{\overline{}}{g}}_{Y}\) times the difference between the transmembrane voltage \(V\) and the Nernst potential for that ion:
\({\overset{\overline{}}{i}}_{Y}\mathbf{=}{\overset{\overline{}}{g}}_{Y}\left( V  E_{Y} \right)\).
This defines a linear currentvoltage relation (‘Ohms law’) as shown in Fig. 14. The gates to be discussed below modify this open channel conductance.
To describe the time dependent transition between the closed and open states of the channel, Hodgkin and Huxley introduced the idea of channel gates that control the passage of ions through a membrane ion channel. If the fraction of gates that are open is y, the fraction of gates that are closed is \(1y\), and a first order ODE can be used to describe the transition between the two states (see Fig. 15):
\(\frac{\text{dy}}{\text{dt}} = \alpha_{y}\left( 1  y \right)  \beta_{y}\text{.y}\)
where \(\alpha_{y}\)is the opening rate and \(\beta_{y}\) is the closing rate.
The solution to this ODE is
\(y = \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}} + Ae^{ \left( \alpha_{y} + \beta_{y} \right)t}\)
The constant \(A\) can be interpreted as \(A = y\left( 0 \right)  \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}}\) as in the previous example and, with \(y\left( 0 \right) = 0\) (i.e. all gates initially shut), the solution looks like Fig. 16(a).
The experimental data obtained by Hodgkin and Huxley for the squid axon, however, indicated that the initial current flow began more slowly (Fig. 16(b)) and they modelled this by assuming that the ion channel had \(\gamma\) gates in series so that conduction would only occur when all gates were at least partially open. Since \(y\) is the probability of a gate being open, \(y^{\gamma}\) is the probability of all \(\gamma\) gates being open (since they are assumed to be independent) and the current through the channel is
\(i_{Y} = {\overset{\overline{}}{i}}_{Y}y^{\gamma} = y^{\gamma}{\overset{\overline{}}{g}}_{Y}\left( V  E_{Y} \right)\)
where \({\overset{\overline{}}{i}}_{Y}{= \overset{\overline{}}{g}}_{Y}\left( V  E_{Y} \right)\) is the steady state current through the open gate.
We can represent this in OpenCOR with a simple extension of the first order ODE model, but in developing this model we will also demonstrate the way in which CellML deals with units.
Note that the decision to deal with units in CellML, rather than just ignoring them or insisting that all equations are represented in dimensionless form, was made in order to be able to check the physical consistency of all terms in each equation5.
There are seven base physical quantities defined by the International d’Unités (SI)6. These are (with their SI units):
length (meter or m)
time (second or s)
amount of substance (mole)
temperature (K)
mass (kilogram or kg)
current (amp or A)
luminous intensity (candela).
All other units are derived from these seven. Additional derived units that CellML defines intrinsically (with their dependence on previously defined units) are: Hz (\(s^{1}\)); Newton, N (\(kg.m.s^{1}\)); Joule, J (\(N.m\)); Pascal, Pa (\(N.m^{2}\)); Watt, W (\(J.s^{1}\)); Volt, V (\(W.A^{1}\)); Siemen, S (\(A.V^{1}\)); Ohm, \(\Omega\) (\(V.A^{1}\)); Coulomb, C (\(s.A\)); Farad, F (\(C.V^{1}\)); Weber, Wb (\(V.s\)); and Henry, H (\(Wb.A^{1}\)). Multiples and fractions of these are defined as follows:
Prefix 
deca 
hecto 
kilo 
mega 
giga 
tera 
peta 
exa 
zetta 
yotta 

Multiples 
Symbol 
da 
h 
k 
M 
G 
T 
P 
E 
Z 
Y 

Factor 
\(10^0\) 
\(10^{1}\) 
\(10^{2}\) 
\(10^{3}\) 
\(10^{6}\) 
\(10^{9}\) 
\(10^{12}\) 
\(10^{15}\) 
\(10^{18}\) 
\(10^{21}\) 
\(10^{24}\) 

Prefix 
deci 
centi 
milli 
micro 
nano 
pico 
femto 
atto 
zepto 
yocto 

Fractions 
Symbol 
d 
c 
m 
μ 
n 
p 
f 
a 
z 
y 

Factor 
\(10^{0}\) 
\(10^{1}\) 
\(10^{2}\) 
\(10^{3}\) 
\(10^{6}\) 
\(10^{9}\) 
\(10^{12}\) 
\(10^{15}\) 
\(10^{18}\) 
\(10^{21}\) 
\(10^{24}\) 
Units for this model, with multiples and fractions, are illustrated in the following CellML Text code:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32  def model first_order_model as
def unit millisec as
unit second {pref: milli};
enddef;
def unit per_millisec as
unit second {pref: milli, expo: 1};
enddef;
def unit millivolt as
unit volt {pref: milli};
enddef;
def unit microA_per_cm2 as
unit ampere {pref: micro};
unit metre {pref: centi, expo: 2};
enddef;
def unit milliS_per_cm2 as
unit siemens {pref: milli};
unit metre {pref: centi, expo: 2};
enddef;
def comp ion_channel as
var V: millivolt {init: 0};
var t: millisec {init: 0};
var y: dimensionless {init: 0};
var E_y: millivolt {init: 85};
var i_y: microA_per_cm2;
var g_y: milliS_per_cm2 {init: 36};
var gamma: dimensionless {init: 4};
var alpha_y: per_millisec {init: 1};
var beta_y: per_millisec {init: 2};
ode(y, t) = alpha_y*(1{dimensionless}y)beta_y*y;
i_y = g_y*pow(y, gamma)*(VE_y);
enddef;
enddef;

The solution of these equations for the parameters indicated above is illustrated in Fig. 17.
The model of a gated ion channel presented here is used in the next two sections for the neural potassium and sodium channels and then in Section 11 for cardiac ion channels. The gates make the channel conductance time dependent and, as we will see in the next section, the experimentally observed voltage dependence of the gating rate constants \(\alpha_{y}\) and \(\beta_{y}\) means that the channel conductance (including the open channel conductance) is voltage dependent. For a partially open channel (\(y < 1\)), the steady state conductance is \(\left( y_{\infty} \right)^{\gamma}{.\overset{\overline{}}{g}}_{Y}\), where \(y_{\infty} = \frac{\alpha_{y}}{\alpha_{y} + \beta_{y}}\). Moreover the gating time constants \(\tau = \frac{1}{\alpha_{y} + \beta_{y}}\) are therefore also voltage dependent. Both of these voltage dependent factors of ion channel gating are important in explaining channel properties, as we show now for the neural potassium and sodium ion channels.
Footnotes
 1
The Brownian motion of individual molecules has energy \(k_{B}T\) (J), where the Boltzmann constant \(k_{B}\) is approximately \(1.34x10^{23}\) (\(J.K^{1}\)). At 25°C, or 298K, \(k_{B}T\) = \(4.10^{21}\) (J) is the minimum amount of energy to contain a ‘bit’ of information at that temperature.
 2
The first law of thermodynamics states that energy is conserved, and the second law (that natural processes are accompanied by an increase in entropy of the universe) deals with the distribution of energy in space.
 3
At infinitely high concentration the specified volume is jammed packed with solute and the entropy is zero.
 4
\(N_{A}\) is Avogadro’s number (\(6.023x10^{23}\)) and is the scaling factor between molecular and macroscopic processes. Boltzmann’s constant \(k_{B}\) and electron charge e operate at the atomic/molecular scale. Their effect at the physiological scale is via the universal gas constant \(R = k_{B}N_{A}\) and Faraday’s constant \(F = eN_{A}\).
 5
It is well accepted in engineering analysis that thinking about and dealing with units is a key aspect of modelling. Taking the ratio of dimensionally consistent terms provides nondimensional numbers which can be used to decide when a term in an equation can be omitted in the interests of modelling simplicity. We investigate this idea further in a later section.
 6