Modeling the kinetics of ligand-gated ion channels

Arnd Roth
Abteilung Zellphysiologie
Max-Planck-Institut fuer Medizinische Forschung
Postfach 10 38 20, D-69028 Heidelberg, Germany
roth@sunny.mpimf-heidelberg.mpg.de

1994 European Mathematica Conference, Oxford, 20.9.1994

Introduction

Picture

Implementation

subunit[t_, u_] :=
{0.2 u + Sin[u]^2 + 0.75 +
  (1 + 0.2 Sin[2 u]) Cos[t + Pi / 4] Sqrt[Cos[u]],
 0.2 u + Sin[u]^2 + 0.75 +
  (1 + 0.2 Sin[2 u]) Sin[t] Sqrt[Cos[u]],
 3 Sin[u]}

rotation[phi_] := {{Cos[phi],  Sin[phi], 0},
                   {-Sin[phi], Cos[phi], 0},
                   {0,         0,        1}}

channel =
ParametricPlot3D[Evaluate[
                 {subunit[t, u],
                  rotation[2 Pi / 5] . subunit[t, u],
                  rotation[4 Pi / 5] . subunit[t, u],
                  rotation[6 Pi / 5] . subunit[t, u],
                  rotation[8 Pi / 5] . subunit[t, u]}],
                 {t, 0, 2 Pi},
                 {u, -Pi / 2 + 10^-9, Pi / 2 - 10^-9},
                 DisplayFunction -> Identity];

membrane1 = Graphics3D[{{Hue[0.65], PointSize[0.03],
Table[ym = y + Mod[x, 1] / 2;
      {Point[{x, ym, 0}],
       Line[{{x - 0.05, ym, 0}, {x - 0.05, ym, -0.85}}],
       Line[{{x + 0.05, ym, 0}, {x + 0.05, ym, -0.85}}]},
      {x, -3, 3, 0.5}, {y, 1.5, 3, 0.5}]
                        },
                        {Hue[0.55], PointSize[0.03],
Table[ym = y + Mod[x, 1] / 2;
      {Point[{x, ym, -2}],
       Line[{{x - 0.05, ym, -2}, {x - 0.05, ym, -1.15}}],
       Line[{{x + 0.05, ym, -2}, {x + 0.05, ym, -1.15}}]},
      {x, -3, 3, 0.5}, {y, 1.5, 3, 0.5}]
                        }}];

membrane2 = Graphics3D[{{Hue[0.65], PointSize[0.03],
Table[ym = y + Mod[x, 1] / 2;
      {Point[{x, ym, 0}],
       Line[{{x - 0.05, ym, 0}, {x - 0.05, ym, -0.85}}],
       Line[{{x + 0.05, ym, 0}, {x + 0.05, ym, -0.85}}]},
      {x, -3, -1.5, 0.5}, {y, -3, 1, 0.5}]
                        },
                        {Hue[0.55], PointSize[0.03],
Table[ym = y + Mod[x, 1] / 2;
      {Point[{x, ym, -2}],
       Line[{{x - 0.05, ym, -2}, {x - 0.05, ym, -1.15}}],
       Line[{{x + 0.05, ym, -2}, {x + 0.05, ym, -1.15}}]},
      {x, -3, -1.5, 0.5}, {y, -3, 1, 0.5}]
                        }}];

Picture

Ligand-gated ion channels are made up of large proteins spanning cell membranes.1 There are many different types of channels; many of them are thought to be composed of five subunits. The following is a tentative picture of such a pentamer embedded in a lipid bilayer. Top: extracellular side, bottom: intracellular side.

Show[channel, membrane1, membrane2, Axes -> False,
     DisplayFunction -> $DisplayFunction]

Function

When they bind neurotransmitter molecules such as glutamate, they may undergo a conformation change and reversibly form a pore several Angstroms in diameter that is permeable for ions like Na+, K+ or Ca2+. The currents flowing through a single channel are in the order of a few picoamperes, but even these small currents can be measured using the patch clamp technique, at a time resolution of a fraction of a millisecond.2 The plot below may illustrate what is typically seen in such an experiment (it shows simulated data, though).

Models

Briefly, the purpose of kinetic models of ion channels is to describe the distribution of the open and closed times (for the case of a single channel under study) or the time evolution of the total current flowing through a group of independent channels caused by a transient in the transmitter concentration.
The experimental findings for the open time distribution of single channels (exponential) and for the spectral density of the fluctuations of current flowing through an ensemble of channels (Lorentzian) are consistent with the assumption that the kinetics of these channels can be modeled as a Markov process. At any given time, the molecule is supposed to be in one of a small number of states such as "singly liganded", "doubly liganded", "open" or "desensitized", the transition rate constants between these states being constant or linear in the transmitter concentration (for those transitions involving the binding of a transmitter molecule to the receptor).3
We turn to a simple example model assuming four states and one transmitter binding site. The states are as follows:

C1 or 1: unliganded or "resting" state
C2 or 2: singly liganded state
C3 or 3: desensitized state (singly liganded)
O   or 4: open, conducting state (singly liganded)
kC2C3 = q23, for example, is the transition rate constant (time-1) for transitions from state C2 to C3.

The time evolution of the occupancies p1, p2, p3, p4 of the states 1, 2, 3, 4 is described by the following first order scheme:

This system of equations can be put in matrix form by defining the transition rate matrix Q according to

We then have

where S is a matrix whose columns are the eigenvectors of QT, L is a diagonal matrix of the eigenvalues of QT and p0 is some initial occupancy of the states of the channel. This matrix exponential solution4 exists only if Q is constant. The mean current flowing through N channels and the variance of the current fluctuations are given by

where V is the potential difference across the membrane, Vrev is the reversal potential (at which the current-voltage relation crosses zero) and g is a vector whose entries are the single channel conductivities of the respective states. That is, we assume Ohmic behavior of the channels.
The nice thing about this formalism is that it can be immediately translated to Mathematica code. However, setting up the Q-matrix for a given model is a bit tedious, so we use a little Mathematica package (which is available on MathSource)5 to facilitate this process.

<<ChannelKinetics`DiagramToQMatrix`

Prepare the following two expressions that will be used as arguments to DiagramToQMatrix[].

diagram = { (* must be rectangular *)
            {C1, connected, C2       , connected, O},
            {  ,          , connected,          ,  },
            {  ,          , C3       ,          ,  }
          };

stateConductivityList = { (* conductivity in pS *)
                          {C1, 0},
                          {C2, 0},
                          {C3, 0},
                          {O,  8}
                        };

Now we call DiagramToQMatrix[] on them and assign the resulting Q-matrix to be the value of q. DiagramToQMatrix[] prompts for the nonzero transition rates. By convention, their unit is ms-1, and concentration-dependent rates contain a factor c[t] (unit: M, mol/liter). The following prompts and user input lines occurred during execution of the example below:

          k[C1, C2] = ? 100000 c[t]
          k[C2, C1] = ? 50
          k[C2, O] = ?  10
          k[O, C2] = ?  0.5
          k[C2, C3] = ? 4
          k[C3, C2] = ? 0.02

q = DiagramToQMatrix[diagram, stateConductivityList]

This example Q-matrix will be used several times in the later sections, so let us display it again in a more readable manner:

MatrixForm[q]

q = {
      {-100000.*c[t], 100000.*c[t], 0, 0}, 
      {50., -64., 4., 10.}, 
      {0, 0.02, -0.02, 0}, 
      {0, 0.5, 0, -0.5}
    }; (* c[t] in M *)

Define the other necessary constants.

n = 25; (* number of available channels *)
v = 0.07; (* membrane potential - reversal pot., Volts *)
g = Last[Transpose[stateConductivityList]]; (* pS *)
p0 = {1, 0, 0, 0}; (* initial occupancy *)

Write the "program" [it might be even simpler... this one is optimized for speed. Once it computes expQt for a particular drug concentration, it remembers the exponential for the next time it is needed.]:

ls[cc_] := ls[cc] = Eigensystem[Transpose[q] /. c[t] -> cc]
l[c_] := l[c] = First[ls[c]] (* eigenvalues *)
s[c_] := s[c] = Last[ls[c]]  (* rows of s are eigenvectors *)

expQt[c_] :=
expQt[c]   = Chop[Transpose[
                            Inverse[s[c]] .
                            DiagonalMatrix[E^(l[c] t)] .
                            s[c]
                  ]
             ]

p   := expQt[0.001] . p0 (* occupancy for c = 0.001 M *)
i   := n v g . p
var := n v^2 ((g^2) . p - (g . p)^2)

That's all!  Now we can plot the occupancy, the mean current, the variance of the current fluctuations etc.

Plot[Evaluate[p], {t, 0, 2}, PlotStyle -> {{}, {Hue[0]},
     {Hue[0.4]}, {Hue[0.8]}}, AxesLabel -> {"t/ms", "p"}]

The current is negative by convention.

Plot[-i, {t, 0, 10}, AxesLabel -> {"t/ms", "I/pA"}]

Plot[var, {t, 0, 10}, AxesLabel -> {"t/ms",
     "var[I]/pA^2"}, PlotRange -> All]

The Exact Solution Package: ColquhounHawkes

Also available on MathSource,5 this package contains a set of functions for numerically integrating the occupancy equation (using NDSolve[], which is very effective on this kind of differential equation) and displaying the results. Since the matrix exponential solution exists only for constant or at least piecewise constant transmitter concentration functions of time, and since this condition is violated in many experiments, it is useful to choose numerical integration as the default method.

<<ChannelKinetics`ColquhounHawkes`

This gives the occupancy predicted for an exponentially decaying transmitter concentration.

myOcc = Occupancy[q, p0, c[t] -> 0.001 E^-t, {t, 10.1}];

The resulting occupancy vector is represented by InterpolatingFunction[] objects. The advantage of this approach is that for a given Q-matrix and transmitter concentration time course, the occupancy equation has to be solved only once, while several (occupancy, current, variance, c.v. etc.) plots can be generated from it. This can save a lot of time. Furthermore, since the solution is stored in a standard Mathematica format, it can be conveniently accessed for further analysis beyond the plots.

OccupancyPlot[myOcc, stateConductivityList, {t, 10}]

This is the mean current. It is negative by convention. The interval of +-1 standard deviation is shown in gray.

CurrentPlot[myOcc, stateConductivityList, n, v, {t, 10}]

One of the observables typically measured in experiments is the time the current takes to rise from 20% to 80% of its peak value. Here, we compute it using FindRoot[].

?CurrentRisetime

CurrentRisetime[myOcc, stateConductivityList, t] *
milliseconds

The Monte Carlo Simulation Package: ColquhounHawkesMC

Monte Carlo simulation is an alternative way of "solving" the occupancy equation. The history of sojourns in the various states is probabilistically simulated for a single channel or a group of individual channels.
The following command loads the MC package. It uses some C routines via MathLink in order to be fast.5

<<ChannelKinetics`ColquhounHawkesMC`

Simulations in ColquhounHawkesMC are done using a random number generator. (The routine used is due to L'Ecuyer, "Efficient and Portable Combined Random Number Generators", Commun. ACM 31, 742-774 (1988), implemented in C by Eberhard v. Kitzing.) When the package is loaded, it takes the time of day as the seed for the random number generator. Two different Mathematica sessions will therefore almost always give different sequences of random numbers and, therefore, different results of the simulations. If this is not desired, the user can explicitly give an integer seed for the random number generator, using CSeedRandom[].

CSeedRandom[6328]

The method of simulation is very simple: the channel starts in some initial state, and the durations of its sojourns are exponentially distributed according to

tau[k_?Positive] := -Log[Random[]] / k
tau[k_ /; k == 0] := Infinity

where k is a transition rate constant. If a channel state allows transitions to several target states, the respective lifetimes are computed and the transition associated with the shortest random lifetime is taken.
The following simulates the response of a single channel to a 1 ms, 1 mM pulse of transmitter concentration.

initialState = C1;
height = 0.001;   (* M *)
duration = 1;     (* ms *)
v = 0.07;         (* V *)
tmax = 10;        (* ms *)
tstep = 0.01;     (* ms *)

myTrace =
SingleTrace[q, stateConductivityList, initialState,
            height, duration, v, {t, tmax, tstep}];

ListPlot[Transpose[{Range[0, tmax, tstep], myTrace}],
         PlotJoined -> True, PlotRange -> All,
         AspectRatio -> 1 / 5,
         AxesLabel -> {"t/ms", "I/pA"},
         PlotLabel -> "Single Channel Current"]

Simulations of the current flowing through an ensemble of channels are done by

?MeanCurrent

The following computes a mean of 200 single channel traces.

n = 200;
myMean = 
MeanCurrent[q, stateConductivityList, initialState,
            height, duration, n, v, {t, tmax, tstep}];

It looks like this.

ListPlot[Transpose[{Range[0, tmax, tstep], myMean}],
         PlotJoined -> True, PlotRange -> All,
         AxesLabel -> {"t/ms", "I/pA"},
         PlotLabel -> "Average Current, 200 Channels"]

It is clear that the examples shown are not exhaustive. There are many possibilities for further analysis of the simulation results. Mathematica allows you to very quickly write functions for such purposes. If some of them turn out to execute too slowly, it is always possible to put their computation-intensive parts in a C routine. That way, you avoid both unnecessary long program coding and execution times, both of which must be considered when assessing the effectiveness of a program.

Semi-automatic model optimization

Procedure

We try to establish a set of states and estimate the transition rate constants for a given type of channel (here, the AMPA/kainate type glutamate receptor channels in CA3 cells of rat hippocampus) from the experimental measurement of the current response of a group of a few hundred channels to various types of steps and pulses in the transmitter concentration as a function of time. The procedure is as follows: we start with an initial guess of a model structure and its rate constants. For any input waveform of the transmitter concentration (pulses, steps and double pulses, the height, duration and distance of which is varied) used in experiment, we compute the response of this model and extract observables such as time constants and peak amplitudes that can be compared with experimental observables. For every observable, the difference between the value predicted by the model and the value measured by experiment is squared and divided by the variance of the experimental value. These terms are summed over all observables to yield the standard chi-square merit function. Next, we search for a set of rate constants that (at least locally) minimizes the chi-square function. This process is repeated for several different sets of initial values of the rate constants in order to explore the location of the possibly numerous local minima. Now, depending on the goodness of fit we either stop optimization, extend the structure of the model or try a completely new structure. In theory, one could automatically generate the set of all model structures containing up to n states and do a parameter optimization on each of them. However, this is far too costly. The only practical method is to look carefully at the model responses generated for the different experiments and guess what feature of the present model constrains further optimization most. Such considerations typically lead to the addition of states to the model. Since the number of free parameters (and therefore the number of states) should be kept at a minimum, we also try to remove states that have become obsolete during optimization. With the new model structure being established, we can start the chi-square optimization over again etc.

Mathematica methods

The computation of the model predictions needed for chi-square optimization and the optimization steps themselves are done in Mathematica, using both the matrix exponential and NDSolve[] solution of the occupancy equation, whichever is faster for the simulation of a given experiment. Time constants of the responses are typically obtained via NonlinearFit[], and peak amplitudes are found using FindMinimum[]. The optimization of the final chi-square function is also done by FindMinimum[], which appears to use a very clever conjugate direction set method for this multidimensional problem. Clearly, the computation of chi-square is very expensive (it takes about 5 to 10 seconds on a Sun Sparcstation 10, Mathematica 2.2), but FindMinimum[] does a good job in minimizing the number of evaluations of chi-square. Since chi-square is being computed using numerical procedures subject to numerical error, and since the outermost FindMinimum[] sometimes tries to compute numerical approximations to the local gradient, it is necessary that the PrecisionGoal and AccuracyGoal settings for the inner numerical procedures be higher than those of the outermost FindMinimum[]. It is definitely worth playing around with different ways of coding the model predictions in Mathematica because a speedup factor of about five is easily possible. For this problem, it means reducing the run time for the full optimization of a structure containing seven states from a few weeks to a few days. Further reduction may be expected from the availability of the new Compile[] facilities that will also accept vectors and matrices.

Results

The resulting model is6

The performance of this model is shown in the following table (format: observable, experimentally measured value, value predicted by the model):

Deactivation time constant, 2.5 ms, 2.5 ms.
Desensitization time constant, 11.3 ms, 11.1 ms.
Maximal open probability (3 mM glutamate), 0.71, 0.72.
Ratio of steady-state current to peak current, 0-5%, 2.9%.
Fraction of channels desensitized directly after a 1 ms pulse, 60 %, 61 %.
Time constant of recovery from desensitization, 48 ms, 48 ms.
Concentration-response EC50, 0.342 mM, 0.355 mM.
Concentration-response Hill coefficient, 1.3, 1.4.
Equilibrium desensitization EC50, 0.0096 mM, 0.0095 mM.
Equilibrium desensitization Hill coefficient, 1.2, 1.2.

For details see the reference.

Discussion

Clearly, we have a complicated protein that has a vast number of states, most of which are short-lived and cannot be seen by this type of experiment. Any model we construct this way must be assumed to be incomplete. It does, however, quite accurately represent the behavior of the channel type investigated for the range of conditions covered by the experiments. It can therefore serve as a building block in detailed simulations of the electrical activity of whole neurons, which are presently being conducted.

Identifiability testing

Associated with the problem of incompleteness, but in some cases amenable to solution using a computer algebra system such as Mathematica is the identifiability and distinguishability testing of such models. For example, if the model has few enough states, we can determine whether there are two (or more) different sets of rate constants k that lead to identical behavior of the system, or if any two different sets of parameters cause a different response.

"When two models M(k) and M(khat) with the same structure M are considered, we shall use the notation M(k) = M(khat) to express that they are output-indistinguishable, i.e. that they have exactly the same external behavior for all admissible inputs. M is structurally locally identifiable if, for almost any admissible value of khat, M(k) = M(khat) implies that k belongs to a countable set (most often finite). When this set reduces to a singleton k = khat, M is structurally globally identifiable."7 8

We apply the method of Raksanyi et al. to our example model. For this purpose we need a symbolic form of the Q-matrix.

q = DiagramToQMatrix[diagram, stateConductivityList]

Define the operators f0 and f1 used for generating the terms of the series expansion of the observables (in our case, p4 = popen) with respect to both time and input around the initial state:

f0 = kC2C1 p2 D[#, p1] +
     (-1 (kC2C1 + kC2C3 + kC2O) p2 +
      kC3C2 p3 + kOC2 p4) D[#, p2] +
     (kC2C3 p2 - kC3C2 p3) D[#, p3] +
     (kC2O p2 - kOC2 p4) D[#, p4] &

f1 = -kC1C2 p1 D[#, p1] + kC1C2 p1 D[#, p2] &

g = p4 (* observable: open probability *)

(* Test: p' = *)
f0[{p1, p2, p3, p4}] + c[t] f1[{p1, p2, p3, p4}]

Now do the same for the model M(khat):

f0hat = kC2C1hat p2hat D[#, p1hat] +
        (-1 (kC2C1hat + kC2C3hat + kC2Ohat) p2hat +
         kC3C2hat p3hat + kOC2hat p4hat) D[#, p2hat] +
        (kC2C3hat p2hat - kC3C2hat p3hat) D[#, p3hat] +
        (kC2Ohat p2hat - kOC2hat p4hat) D[#, p4hat] &

f1hat = -kC1C2hat p1hat D[#, p1hat] +
        kC1C2hat p1hat D[#, p2hat] &

ghat = p4hat

(* Test: phat' = *)
f0hat[{p1hat, p2hat, p3hat, p4hat}] +
c[t] f1hat[{p1hat, p2hat, p3hat, p4hat}]

At time t = 0, all channels are in state C1 whose occupancy is p1.

initialCondition = {p1 -> 1, p2 -> 0, p3 -> 0, p4 -> 0,
                    p1hat -> 1, p2hat -> 0, p3hat -> 0,
                    p4hat -> 0}

We now want to see what the conditions for M(k) = M(khat) are, i.e. given khat, what values of k satisfy this equation? If the only solution is k = khat, then M is identifiable. We will prove now that this is indeed the case.

equation1 =
f1 @ f0 @ g == f1hat @ f0hat @ ghat /. initialCondition

equation2 =
f1 @ f1 @ f0 @ g == f1hat @ f1hat @ f0hat @ ghat /.
initialCondition

equation3 =
f1 @ f0 @ f0 @ g == f1hat @ f0hat @ f0hat @ ghat /.
initialCondition

equation4 =
f1 @ f0 @ f0 @ f0 @ g ==
f1hat @ f0hat @ f0hat @ f0hat @ ghat /. initialCondition

equation11 =
f1 @ f0 @ f0 @ f0 @ f1 @ f0 @ f0 @ g ==
f1hat @ f0hat @ f0hat @ f0hat @ f1hat @
f0hat @ f0hat @ ghat /. initialCondition

equation12 =
f1 @ f1 @ f1 @ f0 @ f0 @ f1 @ f0 @ f0 @ f0 @ f1 @
f0 @ f0 @ g == f1hat @ f1hat @ f1hat @ f0hat @
f0hat @ f1hat @ f0hat @ f0hat @ f0hat @ f1hat @
f0hat @ f0hat @ ghat /. initialCondition

Solve[{equation1, equation2, equation3, equation4,
       equation11, equation12},
      {kC1C2, kC2C1, kC2O, kOC2, kC2C3, kC3C2}]

The only solution is k = khat, and therefore M is identifiable.

Acknowledgements

I thank Professors B. Sakmann and P. Jonas for teaching me ion channels, and I also thank S. Born for many useful discussions on chi-square optimization.

References

1) B. Hille, Ionic Channels of Excitable Membranes, 2nd ed., Sinauer, Sunderland, Mass., 1992.
2) B. Sakmann and E. Neher, Eds. Single-Channel Recording, Plenum, New York, 1983.
3) D. Colquhoun and A.G. Hawkes, Relaxation and fluctuations of membrane currents that flow through drug-operated channels, Proc. Roy. Soc. Lond. B 199, 231 (1977).
4) G. Strang, Linear Algebra and its Applications, 3rd ed., Harcourt Brace Jovanovich, San Diego, 1988.
5) A. Roth, ChannelKinetics: Packages for Modeling Ion Channel Kinetics, MathSource item No. 0205-782 (1993).
6) P. Jonas, G. Major and B. Sakmann, Quantal components of unitary EPSCs at the mossy fibre synapse on CA3 pyramidal cells of rat hippocampus, J. Physiol. 472, 615 (1993).
7) A. Raksanyi, Y. Lecourtier, E. Walter and A. Venot, Identifiability and Distinguishability Testing Via Computer Algebra, Math. Biosci. 77, 245 (1985).
8) E. Walter, Identifiability of State Space Models, Springer, Berlin, 1982.