This document lays out some of the basics of Markov models, and provides a formal mathematical notation.
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
dpdt=Apas used in the master equation rather than the form
dpdt=pQas used in Q-matrix notation. To convert, note that A=QT.
The fundamental assumption is that a channel is always in one of n states Si. Writing x(t) for the channel's state at time t, we can express this as x(t)∈{S1,S2,…,Sn}. The channel can transition instantaneoulsy between states, with rates characterised by rij with i,j∈{1,2,…,n}. This is represented schematically using notation similar to that for chemical reactions, for example:
S1r12⇌r21S2r23⇌r32S3At any time t, the chance of finding the channel in state Si is si(t), such that si∈[0,1] and ∑si=1.
si(t)=P[x(t)=Si]Transition rates rij are strictly non-negative, and depend on the membrane potential V and a set of model parameters 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 Si→Sj implies the existence of a transition Sj→Si. 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 fij(V(t),p_).
rij(t)={fij(V(t),p_)if i≠j0if i=jThe choice of functions fij 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 gmax, and that the channel only conducts a single species of ion. Using an Ohmic driving force (V−E) we find:
I(t)=gmax⋅g(x(t))⋅(V(t)−E)g(x(t))={1if x(t)∈Sconducting0otherwiseThe meaning of the transition rates is defined by postulating that, for a channel in state Si at time t, the probability of a transition Si→Sj occurring in the next infinitesimal interval dt is approximately equal to rij(t)dt.
More accurately:
P[x(t+dt)=Sj|x(t)=Si]=rij(t)dt+o(dt),i≠jWhere o(dt) is an error term that vanishes with decreasing dt (but is otherwise unbounded).
The chance of a change Si→Sj occurring during the interval (t,t+dt] is equal to the chance of the transition multiplied by the probability of being in state Si at time t:
P[x(t)=Si∩x(t+dt)=Sj]=si(t)rij(t)dt+o(dt),i≠jWe can now write an equation for the change in a probability si(t) during the interval dt, by summing the probability of changing to Si from any connected state Sj and subtracting the probability of changing from Si to any state Sj.
si(t+dt)−si(t)=n∑j=1sj(t)rji(t)dt−n∑j=1si(t)rij(t)dt+o(dt)From which we can derive the ordinary differential equation (ODE) form:
ddtsi(t)=n∑j=1(sjrji−sirij)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
S1r12⇌r21S2r23⇌r32S3we find
˙s1=(s2r21−s1r12)˙s2=−(s2r21−s1r12)+(s3r32−s2r23)˙s3=−(s3r32−s2r23)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 fij to obtain a linearised system with time-independent rates rij. We can then define a matrix
R=[rij]and write the system of ODEs as
˙s_=As_=(RT−diag(R1_))s_Here 1_ is a column-vector of n ones and diag is the operator described below.
For the three-state example:
ddt[s1s2s3]=[−r12r210r12−r21−r23r320r23−r32][s1s2s3]=As_and
A=[0r210r120r320r230]+[r12000r21+r23000r32]=RT−diag([r12r21+r23r32])=RT−diag(R1_)diag()
operator¶The diag() operator, when applied to an (n×1) vector u_ creates an (n×n) diagonal matrix D such that Dii=ui.
diag[3514]=(3000050000100004)When applied to an (n×n) matrix A, the diag() operator returns a vector containing only the elements on the diagonal.
diag(7002221863191765)=[7215]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 ∑si=1, we can remove an arbitrary state from the system, solve it, and then calculate the removed state as
si=1−∑j≠isj.
Given a matrix A such that ˙s_=As_, the rate constant matrix R can be obtained using
RT=A−diag(diag(A))Biophysical constraint: If loops then sum both ways is equal, otherwise energy consumption / generation.
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''.
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 Si to equal si exactly.
The conductance can now be defined as the fraction of channels in a conducting state:
G(V(t),t)=∑conductingsiFor the combined current, we find
I=Gmax⋅G(V(t),t)⋅(V(t)−E)Assuming there are m identical (but not synchronised!) channels, the maximum conductance can, in principle, be broken down as
Gmax=m⋅gmaxwhere gmax is the single channel conductance.
However, since both these quantities are typically unknown when a model is created, Gmax 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:
Gmax=αm⋅gmaxAs 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.
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 z_(t) such that zi(t) is the number of channels in state Si at time t.
The initial state of the system is written as z_(t=t0)=z_0
We number every transition Si→Sj so that we can refer to a transition Si→Sj as Rk using the appropriate index k.
Rates are defined as: λk(t,z)=rijzi(t)
The sum of all λk is written as λ(t,z)=∑λij(t,z)
The algorithm proceeds in the following steps:
TODO
In some situations, the membrane potential V(t) is fixed to some (piecewise) constant value and the transition rates rij 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.
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.