How to use the ChannelKinetics
Mathematica   Packages

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

12.12.1994

Summary

The ChannelKinetics Mathematica Packages implement parts of "Relaxation and fluctuations of membrane currents that flow through drug-operated channels" by D. Colquhoun & A.G. Hawkes, Proc. R. Soc. Lond. B 199, 231-262 (1977).

DiagramToQMatrix builds Colquhoun's Q-matrix from a transition diagram, prompting for nonzero transition rate constants. It automatically checks cyclic diagrams for microscopical reversibility.


ColquhounHawkes contains various functions for

- computing the occupancy of the different states of a channel using the Runge-Kutta algorithm, which allows for arbitrary transmitter concentration functions of time,

- for plotting the mean current as a function of time,
- its variance, coefficient of variation and third central moment,
- and for computing the 20%-80% risetime of the current.


ColquhounHawkesMC is an analog of ColquhounHawkes based on Monte Carlo simulation of the channel state transitions. It is restricted to square pulse transmitter concentration functions.

Additionally, it contains a function that simulates experiments for determinimg the peak EPSC amplitude distribution for an ensemble of channels. This includes the effects of gaussian lowpass filtering and averaging over a certain interval of time during the peak.

The most recent version of ColquhounHawkesMC, version 1.2, installs and uses some C routines collected in the file "chmc". This file must reside in the same directory as ColquhounHawkesMC, as specified in Appendix A. If this is not possible on certain computer systems, you can install "chmc" manually before you load the package. Alternatively, you might want to change ColquhounHawkesMC's Mathematica code portion that is responsible for installing the C routines. It is located directly below the definition of the error messages for the exported functions.See Appendix B for information on how to build "chmc" on your computer system.

The Transition Diagram To Q-Matrix Conversion Package: DiagramToQMatrix

This package defines a single function with the same name, DiagramToQMatrix[].
We turn right away to a simple example, because that is the best way to show what DiagramToQMatrix[] does.

This loads the package. (If you type a line like this, do not change spelling, capitalization, nor the "`" character.)

<<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 you for the nonzero transition rates. By convention, their unit is ms^-1, and concentration-dependent rates contain a factor c[t].
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]

Now we can always refer to this Q-matrix as q.

q

For example, it might be convenient to display it with its entries nicely arranged in a rectangular array.

MatrixForm[q]

Often there will be several different calculations that use the same Q-matrix. If we store q in its input form, we can reuse it later simply by evaluating this input form, without rebuilding it from the transition diagram, reentering all the rates by hand. Also, the DiagramToQMatrix package is not needed in this case.
We accomplish this by invoking InputForm[] and changing the "cell style" of the result from output to input (in the Style/Cell Style menu). One might also rearrange the terms a bit and prepend a q = to the matrix.

InputForm[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 *)

For more information on the structure of the arguments that are passed to DiagramToQMatrix[], type

?DiagramToQMatrix

This is a precise but not very readable description, so it might be useful to look at the example diagram above in order to see what it means. The rectangularity of the diagram is essential, i.e., each line must have the same number of commas in it.

The Exact Solution Package: ColquhounHawkes

This loads the package.

<<ChannelKinetics`ColquhounHawkes`

You can get a list of the functions defined in this package by typing

Help[]

?SquarePulseConcentration

We construct a square pulse transmitter concentration function with a pulse length of 1 millisecond and a pulse height of 1 mM. We call it mySquare, for example.

mySquare = SquarePulseConcentration[t, 0.001, 1]

Copy q and stateConductivityList from the section above, and newly add in an initialOccupancy such that at t = 0, all channels are in state C1.

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 *)

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

initialOccupancy = {1, 0, 0, 0};

The order of the entries of the occupancy vector initialOccupancy corresponds to the order of the entries of stateConductivityList. For example, this setting of initialOccupancy means that the initial occupancy of state C1 is unity, and the occupancies of the states C2, C3 and O are all zero. The same ordering is used for the occupancy vector returned by the Occupancy[] function:

?Occupancy

The actual computation of the occupancies as a function of time is done by this function.

myOcc = Occupancy[
q, initialOccupancy, c[t] -> mySquare, {t, 10.1}];

The resulting occupancy vector is represented by InterpolatingFunction[] objects. They provide a spline interpolation of the occupancies of the states of the channel in the time range they have been computed in (here, from 0 to 10.1 milliseconds). 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. As a simple example, we check for the initial conditions to be satisfied:

initialOccupancy == Chop[myOcc /. t -> 0]

We now proceed to generating some plots from myOcc. First, the occupancies of the different states themselves.

?OccupancyPlot

This makes a color plot.

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

The color plot looks nice on the screen, but on a black-and-white printer, dashed lines might be easier to distinguish. The dashing definition is given as an option to OccupancyPlot[]:

printerOption =
PlotStyle -> {
               {},
               {Dashing[{0.05, 0.05}]},
               {Dashing[{0.01, 0.01}]},
               {Dashing[{0.01, 0.05}]},
               {Dashing[{0.01, 0.02, 0.05, 0.02}]}
             };

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

Next, we plot the current. It is negative by convention.

?CurrentPlot

n = 10;    (* number of available channels *)
u = 0.07; (* Volts, voltage across the membrane *)

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

This gives a current plot covering a smaller time interval.

CurrentPlot[myOcc, stateConductivityList,
            n, u, {t, 1}]

?CurrentRisetime

Any implementation of CurrentRisetime[] must rely on numerical root finding. This requires the routine to guess what starting values for the secant method will lead it to converge in the right place. Such guesswork is fallible if the current function is a sufficiently pathological one. Our example, though, is well-behaved so there should be no problems. (You can always verify the result of CurrentRisetime[] by looking at the current plot.)

CurrentRisetime[myOcc, stateConductivityList, t] *
milliseconds

?CurrentVariancePlot

This plots the variance of the current, in pA^2. Note that the argument sequence to most of the plotting functions is the same. We take the n and u from the current plot above.

CurrentVariancePlot[myOcc, stateConductivityList,
                    n, u, {t, 10}]

?CoefficientOfVariationPlot

CoefficientOfVariationPlot[myOcc, stateConductivityList,
                           n, u, {t, 10}]

?CurrentThirdMomentPlot

CurrentThirdMomentPlot[myOcc, stateConductivityList,
                       n, u, {t, 10}]

First Latencies

?FirstLatencies

myFL = FirstLatencies[q, stateConductivityList,
                      initialOccupancy,
                      c[t] -> 0.001, {t, 1000}];

Plot[myFL, {t, 0, 1}, AxesLabel -> {"t/ms",
     "PDF"}, PlotLabel -> "First Latencies"]

Check the normalization:

NIntegrate[myFL, {t, 0, 1000}, MaxRecursion -> 15]

Double Pulses

All this was for a single square pulse transmitter concentration function, 1ms long, 1 mM high. But the Occupancy[] function allows for arbitrary transmitter concentration functions - exponentials, gaussians, multiple pulses, and others. So let us try a double pulse calculation.

Construct a double square pulse transmitter concentration function with a pulse length of 1 millisecond, a pulse height of 1 mM, and a pulse distance of 8 milliseconds. We call it myDoublePulse, for example.

myDoublePulse =
SquarePulseConcentration[t, 0.001, 1] +
SquarePulseConcentration[t - 9, 0.001, 1]

See what the concentration function looks like. Plot[] takes many options whose description in Wolfram is worth looking at (p. 137-143).

Plot[myDoublePulse, {t, -1, 11},
     PlotRange -> All,
     AxesOrigin -> {-1, 0},
     AxesLabel -> {"t/ms", "c[t]/M"},
     PlotLabel -> "Transmitter Concentration"]

Compute the occupancy for this double pulse concentration function in the time interval from 0 to 20.1 milliseconds.

myOcc = Occupancy[
q, initialOccupancy, c[t] -> myDoublePulse, {t, 20.1}];

But there is a possible pitfall. The stepsize adaption facility of the Runge-Kutta algorithm employed by Occupancy[] effectively assumes the transmitter concentration function to vary smoothly with time. Although the square pulse functions are essentially non-smooth, usually there is no problem at all. But in the case of two pulses, if the pulse distance is too great, the algorithm may miss the second pulse. If this happens, you must increase the accuracy demands on the algorithm using the options PrecisionGoal, AccuracyGoal, and MaxSteps. (PrecisionGoal -> n is an option to some numerical functions that states that the result of a computation should have n significant digits. AccuracyGoal -> n states that the result of a computation should have n significant digits at the right of the decimal point. MaxSteps -> limit causes Occupancy[] to terminate if it takes more than limit steps in t.) Specifying Occupancy[..., {t, 20.1}, PrecisionGoal -> 8, AccuracyGoal -> 8, MaxSteps -> 2000] or similar values should force the algorithm to use a stepsize small enough not to miss the second pulse. In Mathematica version 2.2, there is a new option MaxStepSize that can also be used to enforce small enough steps.

The new occupancy plot looks like this:

OccupancyPlot[myOcc, stateConductivityList,
              {t, 20}, printerOption]

We see that, even with no increased setting of AccuracyGoal and PrecisionGoal, the second pulse has not been missed in this case.

Open channel (and other) blockers can be modeled easily:

Needs["ChannelKinetics`DiagramToQMatrix`"]

Add a blocked state B to the transition diagram:

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

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

initialOccupancy = {1, 0, 0, 0, 0};

qB = DiagramToQMatrix[diagram, stateConductivityList]

Now we can always refer to this Q-matrix as qB.

qB

MatrixForm[qB]

qB = {{-100000.*c[t], 100000.*c[t], 0, 0, 0},
      {50., -64., 4., 10., 0},
      {0, 0.02, -0.02, 0, 0},
      {0, 0.5, 0, -0.5 - 50000.*b[t], 50000.*b[t]},
      {0, 0, 0, 0.01, -0.01}};

We construct two square pulse transmitter concentration functions: one with a pulse length of 1 millisecond and a pulse height of 1 mM (glutamate), and another with a pulse length of 10.1 ms and a height of 0.01 mM (blocker).

glutamate = SquarePulseConcentration[t, 0.001, 1]

blocker = SquarePulseConcentration[t, 0.00001, 10.1]

The actual computation of the occupancies as a function of time is done by this function.

myOcc = Occupancy[qB, initialOccupancy,
                  {c[t] -> glutamate, b[t] -> blocker},
                  {t, 10.1}];

This makes a color plot.

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

Next, we plot the current. It is negative by convention.

n = 10;    (* number of available channels *)
u = 0.07; (* Volts, voltage across the membrane *)

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

CurrentRisetime[myOcc, stateConductivityList, t] *
milliseconds

The Monte Carlo Simulation Package: ColquhounHawkesMC

Load the package. You should start a fresh Mathematica session before loading this package.

<<ChannelKinetics`ColquhounHawkesMC`

You can get a list of the functions defined in this package by typing

Help[]

Copy q and stateConductivityList from the section above, and newly add in the variable initialState such that at t = 0, all channels are in state C1. Note that initialState is not the same type of expression as initialOccupancy.

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 *)

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

initialState = C1;

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 you want to make sure that you always get the same sequence of pseudorandom numbers, you can explicitly give an integer seed for the random number generator, using CSeedRandom[].

?CSeedRandom

An example call.

CSeedRandom[6328]

CSeedRandom[] gives a seed to the random number generator called by the C routines. It is different from Mathematica's built-in SeedRandom[] which resets the seed for Mathematica's own random number generator.

Now we turn to doing some simulations. First, we try simulation of the current flowing through a single channel.

?SingleTrace

Prepare some variables for the simulation (alternatively, you might wish to give these values directly as arguments to SingleTrace[]):

height = 0.001; (* M *)
duration = 1;   (* ms *)
u = 0.07;       (* V *)
tmax = 10;      (* ms *)
tstep = 0.2;    (* ms *)

Invoke SingleTrace[].

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

This is a list of numbers that represents the single channel current at sample points with distance tstep milliseconds. We plot it using the built-in Mathematica function ListPlot[].

ListPlot[myTrace]

We can make this look better by giving some options to ListPlot[]:

ListPlot[myTrace, PlotJoined -> True, PlotRange -> All,
         AspectRatio -> 1 / 5,
         AxesLabel -> {"t/samples", "I/pA"},
         PlotLabel -> "Single Channel Current"]

A smaller setting for tstep will make the transitions steeper.

tstep = 0.01; (* ms *)
myTrace =
SingleTrace[q, stateConductivityList, initialState,
            height, duration, u, {t, tmax, tstep}];
ListPlot[myTrace, PlotJoined -> True, PlotRange -> All,
         AspectRatio -> 1 / 5,
         AxesLabel -> {"t/samples", "I/pA"},
         PlotLabel -> "Single Channel Current"]

We can easily change the unit of the x-axis to be milliseconds.

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

Using TraceTable[], we generate a list of traces.

?TraceTable

n = 3; (* 3 traces *)

myTraces =
TraceTable[q, stateConductivityList, initialState,
           height, duration, n, u, {t, tmax, tstep}];

This error message is not always useful.

We map ListPlot[] on every element of the trace list.

Map[
    ListPlot[Transpose[{Range[0, tmax, tstep], #}],
             PlotJoined -> True, PlotRange -> {-0.6, 0},
             AspectRatio -> 1 / 5,
             AxesLabel -> {"t/ms", "I/pA"},
             PlotLabel -> "Single Channel Current"] &,
    myTraces];

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

?MeanCurrent

MeanCurrent[], as well as some other functions, take the same sequence of arguments as TraceTable[].
The following computes a mean of 200 single channel traces.

n = 200;
myMean = 
MeanCurrent[q, stateConductivityList, initialState,
            height, duration, n, u, {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"]

This gives the sum of 200 single channel traces.

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

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

Write the simulated EPSC to a file, with the numbers separated by the End-Of-Line character.

ColumnForm[myEPSC] >> file

The following three functions have close analogs in the ColquhounHawkes exact solution package. We won't discuss them here, but give their help texts.

?CurrentVariance

?CoefficientOfVariation

?CurrentThirdMoment

FilterTest[] is useful for determinimg how the chosen filter acts on a simulated EPSC. You can visually check the effects of different cutoff and sampling frequencies.

?FilterTest

A good choice for samplingfreq and for the cutoff frequency fc is 40 kHz and 3 kHz, respectively.

height = 0.001;    (* M *)
duration = 1;      (* ms *)
u = 0.07;          (* V *)
tmax = 10;         (* ms *)
samplingfreq = 40; (* kHz *)
fc = 3;            (* kHz *)

FilterTest[EnsembleSize -> 20,
           q, stateConductivityList, initialState,
           height, duration, u,
           {t, tmax, samplingfreq, fc}]

The construction of a histogram for the distribution of the peak EPSC amplitudes is done by

?PeakEPSCAmplitudeDistribution

EnsembleSize -> m defines the number of available channels, while n is the number of trials. In the following two simulations we use an integration interval of length 0.5 ms around the peak, that is, from (time_of_the_peak) - 0.25 ms to (time_of_the_peak) + 0.25 ms.

tint = 0.5;  (* ms *)

m100n200 = (* 100 available channels, 200 trials *)
PeakEPSCAmplitudeDistribution[EnsembleSize -> 100,
q, stateConductivityList, initialState,
height, duration, 200, u,
{t, tmax, samplingfreq, fc, tint}]

m200n121 = (* 200 available channels, 121 trials *)
PeakEPSCAmplitudeDistribution[EnsembleSize -> 200,
q, stateConductivityList, initialState,
height, duration, 121, u,
{t, tmax, samplingfreq, fc, tint}]

Once the distributions have been computed, you can plot them in various ways.

?DistributionPlot

DistributionPlot[m100n200]

DistributionPlot[m200n121]

DistributionPlot[m100n200, m200n121]

Conclusion

It is clear that the examples shown in this manual 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.

Appendix A: File Structure

ChannelKinetics consists of the following files:

Unix (Sun)

/usr/local/math/Packages/ChannelKinetics/ColquhounHawkes.m
/usr/local/math/Packages/ChannelKinetics/ColquhounHawkesMC.m
/usr/local/math/Packages/ChannelKinetics/DiagramToQMatrix.m
/usr/local/math/Packages/ChannelKinetics/NeuronInterface.m
/usr/local/math/Packages/ChannelKinetics/README.txt
/usr/local/math/Packages/ChannelKinetics/chmc
/usr/local/math/Packages/ChannelKinetics/chmc.tm

MacOS (Macintosh)

:Mathematica 2.2:Packages:ChannelKinetics:
   ColquhounHawkes.m
   ColquhounHawkesMC.m
   DiagramToQMatrix.m
   NeuronInterface.m
   README.txt
   chmc
   chmc Folder:
      chmc.\pi
      chmc.tm
      chmc.tm.c
      chmc.\pi.rsrc
      chmc

Microsoft Windows (PC)

C:\WNMATH22\PACKAGES\CHANNELK\COLHAW.M
C:\WNMATH22\PACKAGES\CHANNELK\COLHAWMC.M
C:\WNMATH22\PACKAGES\CHANNELK\DIAGRAMT.M
C:\WNMATH22\PACKAGES\CHANNELK\NEURONIN.M
C:\WNMATH22\PACKAGES\CHANNELK\README.TXT
C:\WNMATH22\PACKAGES\CHANNELK\CHMC.EXE
C:\WNMATH22\PACKAGES\CHANNELK\CHMC.TM
C:\WNMATH22\PACKAGES\CHANNELK\CHMCTM.C

Appendix B: Building chmc

Unix (Sun)

Do a cd to the Packages directory. Then simply type mcc chmc.tm -o chmc -O. See the man page on mcc for more information. MathLink on Unix workstations is very easy to use. (It is much harder to compile chmc.tm on Macintosh or even Windows systems, but this is not my fault.)

MacOS (Macintosh)

See the ReadMe Notebook in your MathLink directory. It describes the process of building MathLink programs for the Macintosh. Use 32-bit integers if possible. If you are using Think C, this also requires recompiling the ANSI library using the 32-bit integers option.

Microsoft Windows (PC)

See the file RELNOTES.WRI that comes with your MathLink for Windows Developer's Kit. Use the "large memory model". Be sure to proceed step by step as described in RELNOTES.WRI. Unfortunately, you can't use the MathLink Developer's Kit that came with Mathematica 2.2.3 for Windows since it has a bug - you need to use a post-2.2.3 MathLink for Windows Developer's Kit.