Markov models of ionic channels

Currently, Brian is designed mainly for Hodgkin-Huxley type models. It is possible to write a Markov model of ionic channel, but it is neither convenient nor very accurate. Here I am only considering deterministic (not stochastic) models.

Syntax

A Markov model consists of states and transition rates between them, which are typically voltage-dependent, but could also depend on concentrations (ions or ligands). Here I only consider models that are expressed as a graph between states. For example, one might consider the following two-state model:

C <-> O

with two transition rates (in Hz). To each state corresponds a variable between 0 and 1 representing the probability that the channel is in that state. Therefore I propose the following syntax:

C->O, exp(a*v)

O->C, exp(b*v)

That is, each line defines a transition between two states, with its transition rate in Hz. The full graph together with the list of states is obtained when all lines are gathered. We may replace the comma with some other delimiter (for example space). Then the corresponding variables C and O are created. For example, there could be a membrane equation as follows:

dv/dt = g*O*(E-v) : volt

Gathering

First, all transitions must be gathered. We obtain a graph of transitions. There could be several channels. These could simply be separated by separating the graph into connected components. For each channel, there are equations equivalent to the list of transitions. For example, in the case of the two-channel model:

dO/dt = exp(a*v)*C - exp(b*v)*O : 1

dC/dt = exp(b*v)*O - exp(a*v)*C : 1

These are actually a bit redundant since O+C=1. This is a general property that the sum of all states is 1. This will be important for integration.

Integration

A simple option would be to generate the equations and let Brian deal with them. However, it is quite likely that the sum of all state variables will drift from 1. One simple way to solve this problem is to generate equations for n-1 variables, and add a static equation for the last one, e.g.:

C = 1-O

Then ideally, we would like to have an appropriate integration method. We notice that the equations are linear with respect to the state variables:

dS/dt = M(v).S

Therefore, a good approach would be to use exact linear integration with fixed v. It is rather unlikely that this can be done symbolically at initialization, and therefore it seems more appropriate to use the exponential matrix to solve it:

S(t+dt) = exp(M(v).dt).S(t)

It might be possible that the structure of the matrix helps. It has the following property: sum(M(i,j),i=1..,) = 0 for every j.

Perspective

A natural generalization is with binding models. For example, a calcium-inactivated channel would look like this:

C <-> O

O + [Ca] <-> I

Then we would also need to introduce a syntax for concentrations. It needs a little bit of thought.