Technical note: Markov channel simulation¶

This document lays out some of the basics of Markov models, and provides a formal mathematical notation.

Important note about terminology¶

The term "Markov model" is used in this document to describe what is known in physics and chemistry as a kinetic scheme or master equation. It is subtly different from the way "Markov model" is used in probabilistics, which is the more common form in single-channel analysis (following the work of e.q. Colquhoun et al.). In particular, we use the form

$$\frac{dp}{dt} = A p$$

as used in the master equation rather than the form

$$\frac{dp}{dt} = p Q$$

as used in Q-matrix notation. To convert, note that $A = Q^{T}$.

Single-channel models¶

The fundamental assumption is that a channel is always in one of $n$ states $S_i$. Writing $x(t)$ for the channel's state at time $t$, we can express this as $x(t) \in \{S_1, S_2, \ldots, S_n\}$. The channel can transition instantaneoulsy between states, with rates characterised by $r_{ij}$ with $i,j \in \{1,2,\ldots,n\}$. This is represented schematically using notation similar to that for chemical reactions, for example:

$$\begin{equation*} S_1 \underset{r_{21}}{\stackrel{r_{12}}{\rightleftharpoons}} S_2 \underset{r_{32}}{\stackrel{r_{23}}{\rightleftharpoons}} S_3 \end{equation*}$$

At any time $t$, the chance of finding the channel in state $S_i$ is $s_i(t)$, such that $s_i \in [0,1]$ and $\sum s_i = 1$.

$$\begin{equation} s_i(t) = P[x(t) = S_i] \end{equation}$$

Transition rates $r_{ij}$ are strictly non-negative, and depend on the membrane potential $V$ and a set of model parameters $\underline{p}$. While $V$ may change over time, the rates themselves are not explicitly time-dependent. We assume that all transitions are reversible, so that the existence of any transition $S_i \rightarrow S_j$ implies the existence of a transition $S_j \rightarrow S_i$. The transition rate from a state to itself is defined as zero. All other rates are either zero (if the states aren't connected) or given by some model-specific function $f_{ij}\big(V(t), \underline{p}\big)$.

$$\begin{equation} r_{ij}(t) = \begin{cases} f_{ij}\left(V(t), \underline{p}\right) & \text{if $i \neq $j} \\ 0 & \text{if $i = j$} \end{cases} \end{equation}$$

The choice of functions $f_{ij}$ changes per model and is sometimes \emph{ad-hoc} and sometimes based on physical assumptions. Similarly, the number of states and connections is chosen freely by the modeler.

One or multiple of the model's states allow it to conduct ionic currents. Here, we assume that all conducting states have the same fixed conductance $g_{max}$, and that the channel only conducts a single species of ion. Using an Ohmic driving force $(V - E)$ we find:

$$\begin{eqnarray} I(t) & = & g_\text{max} \cdot g(x(t)) \cdot (V(t) - E) \\ g(x(t)) & = & \begin{cases} 1 & \text{if } x(t) \in S_\text{conducting} \\ 0 & \text{otherwise} \end{cases} \end{eqnarray}$$

State transitions¶

The meaning of the transition rates is defined by postulating that, for a channel in state $S_i$ at time $t$, the probability of a transition $S_i \rightarrow S_j$ occurring in the next infinitesimal interval $dt$ is approximately equal to $r_{ij}(t) dt$.

More accurately:

$$\begin{equation} P\Big[x(t + dt) = S_j \,|\, x(t) = S_i\Big] = r_{ij}(t) dt + o(dt) ,\quad\quad i \neq j \end{equation}$$

Where $o(dt)$ is an error term that vanishes with decreasing $dt$ (but is otherwise unbounded).

The chance of a change $S_i \rightarrow S_j$ occurring during the interval $(t, t + dt]$ is equal to the chance of the transition multiplied by the probability of being in state $S_i$ at time t:

$$\begin{equation} P\Big[x(t) = S_i \cap x(t + dt) = S_j \Big] = s_i(t)r_{ij}(t) dt + o(dt) ,\quad\quad i \neq j \end{equation}$$

We can now write an equation for the change in a probability $s_i(t)$ during the interval $dt$, by summing the probability of changing to $S_i$ from any connected state $S_j$ and subtracting the probability of changing from $S_i$ to any state $S_j$.

$$\begin{equation} s_i(t + dt) - s_i(t) = \sum_{j=1}^n{s_j(t)r_{ji}(t)dt} - \sum_{j=1}^n{s_i(t)r_{ij}(t)dt} + o(dt) \end{equation}$$

From which we can derive the ordinary differential equation (ODE) form:

$$\begin{equation} \label{eq:derived-ode-form} \frac{d}{dt}s_i(t) = \sum_{j=1}^n({s_jr_{ji} - s_ir_{ij}}) \end{equation}$$

Formulation as an ODE¶

Using the equation above we can write a system of ordinary differential equations (ODEs) that describe the probability of the system being in each state over time.

For the three-state example

$$\begin{equation*} S_1 \underset{r_{21}}{\stackrel{r_{12}}{\rightleftharpoons}} S_2 \underset{r_{32}}{\stackrel{r_{23}}{\rightleftharpoons}} S_3 \end{equation*}$$

we find

$$\begin{equation} \begin{array}{ccccc} \dot{s}_{1} & = & \,\,\,\,\,\left(s_{2}r_{21}-s_{1}r_{12}\right)\\ \dot{s}_{2} & = & -\left(s_{2}r_{21}-s_{1}r_{12}\right) & + & \,\,\,\,\,\left(s_{3}r_{32}-s_{2}r_{23}\right)\\ \dot{s}_{3} & = & & & -\left(s_{3}r_{32}-s_{2}r_{23}\right) \end{array} \end{equation}$$

The transition rates are typically calculated using non-linear functions of $V$.

However, for situations where $V$ is (piecewise) constant, we can evaluate the functions $f_{ij}$ to obtain a linearised system with time-independent rates $r_{ij}$. We can then define a matrix

$$\begin{equation} R = [r_{ij}] \end{equation}$$

and write the system of ODEs as

$$\begin{equation} \dot{\underline{s}} = A \underline{s} = \left( R^{T} - \text{diag}\left(R\underline{1}\right) \right) \underline{s} \end{equation}$$

Here $\underline{1}$ is a column-vector of $n$ ones and $\text{diag}$ is the operator described below.

For the three-state example:

$$\begin{equation} \frac{d}{dt} \left[\begin{array}{c} s_1 \\ s_2 \\ s_3 \\ \end{array}\right] = \left[\begin{array}{ccc} -r_{12} & r_{21} & 0 \\ r_{12} & -r_{21} - r_{23} & r_{32} \\ 0 & r_{23} & -r_{32} \\ \end{array} \right] \left[ \begin{array}{c} s_1 \\ s_2 \\ s_3 \\ \end{array} \right] = A\underline{s} \end{equation}$$

and

$$\begin{align*} A &= \left[\begin{array}{ccc} 0 & r_{21} & 0 \\ r_{12} & 0 & r_{32} \\ 0 & r_{23} & 0 \\ \end{array}\right] + \left[\begin{array}{ccc} r_{12} & 0 & 0 \\ 0 & r_{21} + r_{23} & 0 \\ 0 & 0 & r_{32} \\ \end{array}\right] \\ &= R^{T} - \text{diag}\left(\left[\begin{array}{ccc} r_{12} \\ r_{21} + r_{23} \\ r_{32} \\ \end{array}\right]\right) \\ &= R^{T} - \text{diag}\left(R\underline{1}\right) \\ \end{align*}$$

The diag() operator¶

The $\text{diag}()$ operator, when applied to an $(n \times 1)$ vector $\underline{u}$ creates an $(n \times n)$ diagonal matrix $D$ such that $D_{ii} = u_i$.

$$\begin{equation} \text{diag}\left[ \begin{array}{c} 3\\5\\1\\4 \end{array} \right] = \left( \begin{array}{cccc} 3 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 4 \\ \end{array} \right) \end{equation}$$

When applied to an $(n \times n)$ matrix $A$, the $\text{diag}()$ operator returns a vector containing only the elements on the diagonal.

$$\begin{equation} \text{diag}\left( \begin{array}{cccc} 7 & 0 & 0 & 2 \\ 2 & 2 & 1 & 8 \\ 6 & 3 & 1 & 9 \\ 1 & 7 & 6 & 5 \\ \end{array} \right) = \left[ \begin{array}{c} 7\\2\\1\\5 \end{array} \right] \end{equation}$$

$A$ is an indefinite matrix¶

Note that the sum of the values in each column of $A$ is always zero. This means that the system is, by definition, indefinite. However, since $\sum{s_i}=1$, we can remove an arbitrary state from the system, solve it, and then calculate the removed state as

$$\begin{equation} s_i = 1 - \sum_{j \neq i}{s_j} \end{equation}$$

.

Obtaining $R$ from $A$¶

Given a matrix $A$ such that $\dot{\underline{s}} = A \underline{s}$, the rate constant matrix $R$ can be obtained using

$$\begin{equation} R^{T} = A - \text{diag}\left(\text{diag}\left(A\right)\right) \end{equation}$$

Loops: Microscopic reversibility¶

Biophysical constraint: If loops then sum both ways is equal, otherwise energy consumption / generation.

Solving the ODE¶

For proof, see "Kattrin Arning (2009) Mathematical Modelling and Simulation of Ion Channels". For pairwise different eigenvalues, the solution ``converges to the equilibrium solution as a sum of n-1 exponentials''.

Modeling whole-cell current¶

The same system of ODEs used to describe single ion channels can be used to describe the aggregated "whole-cell" current created when multiple identical ion channels are measured at the same time. To do this, we assume that the number of channels in the cell is large enough for the fraction of channels in each state $S_i$ to equal $s_i$ exactly.

The conductance can now be defined as the fraction of channels in a conducting state:

$$\begin{equation} G\big(V(t), t\big) = \sum_\text{conducting} s_i \end{equation}$$

For the combined current, we find

$$\begin{equation} I = G_\text{max} \cdot G\big(V(t), t\big) \cdot \big(V(t) - E\big) \end{equation}$$

Assuming there are $m$ identical (but not synchronised!) channels, the maximum conductance can, in principle, be broken down as

$$\begin{equation*} G_\text{max} = m \cdot g_\text{max} \end{equation*}$$

where $g_\text{max}$ is the single channel conductance.

However, since both these quantities are typically unknown when a model is created, $G_\text{max}$ is commonly set to scale the current to an appropriate size. This may lead to the inclusion of an unknown scaling constant in its definition, that should really have been applied to the transition rates:

$$\begin{equation*} G_\text{max} = \alpha m \cdot g_\text{max} \end{equation*}$$

As a result, it is very dangerous to draw conclusions about the number of channels or single channel conductance from a current model's "maximum conductance" constant.

Discrete simulation¶

Gillespie's algorithm can be used to simulate a group of $m$ channels with stochastic behavior, under the assumption of a constant membrane potential $V$.

First, we define a state vector $\underline{z}(t)$ such that $z_i(t)$ is the number of channels in state $S_i$ at time $t$.

The initial state of the system is written as

$$\underline{z}(t=t_0) = \underline{z}_0$$

We number every transition $S_i \rightarrow S_j$ so that we can refer to a transition $S_i \rightarrow S_j$ as $R_k$ using the appropriate index $k$.

Rates are defined as:

$$\lambda_{k}(t, z) = r_{ij} z_i(t)$$

The sum of all $\lambda_{k}$ is written as

$$\lambda(t, z) = \sum{\lambda_{ij}(t, z)}$$

The algorithm proceeds in the following steps:

1. Initialise the system to $t = t_0$ and $\underline{z} = \underline{z}_0$
2. Calculate the rates $\lambda_{k}$ and the sum $\lambda$.
3. Determine the time $\tau$ until the next transition (of any kind) by sampling from an exponential distribution with mean $1/\lambda$.
4. Determine the type of transition by drawing a random number $r_2$ from a uniform distribution on the interval $[0, \lambda)$. Choose transition $T_1$ if $0 \leq r_2 < \lambda_1$, choose $T_2$ if $\lambda_1 \leq r_2 < \lambda_2$ and so on.
5. Update the state by moving one channel through the selected transition. Increase the time variable by $\tau$.
6. Repeat steps 2 to 5 until $t > t_\text{end}$.

Isolating Markov models from cell ephys models¶

TODO

Analytical solution¶

In some situations, the membrane potential $V(t)$ is fixed to some (piecewise) constant value and the transition rates $r_{ij}$ become constant. In such cases, for example when simulating a voltage-clamp experiment, the whole-cell currents can be calculated by solving the system of ODEs using eigenvalue decomposition.

Whole-cell simulation¶

If $V(t)$ cannot be fixed, for example when simulating an action potential or an imperfect (i.e. realistic) voltage-clamp experiment, the system of ODEs to solve is non-linear. In these cases a numerical quadrature method is used, for example the forward Euler method, an adaptive Runge-Kutta method or advanced methods such as CVODE.