Because the timing of individual openings and the validity for small copy numbers mattered in my application, I ended up implementing a Gillespie algorithm in a NMODL mechanism. It reproduces the kinetics of the mean field variant in the average of many reproductions and also for many channels. Surprisingly, it does not limit the simulation time. Probably, because even with 100 channels, no transition occurs in most of the time steps. Therefore, after the first random number is drawn, the update is over and in which case the mechanism is faster than any implementation of a differential equation solver or any mechanism that always updates all channels separately.

Because it runs orders of magnitude faster than some of the brute-force implementations around and because I am not aware of a Gillespie algorithm implemented in NMODL, I thought it might could serve as a useful reference, even though the code is specific for two-state models and needs some extension for a general case.

Code: Select all

```
: an ion channel with gating that depends on an external variable other than voltage
: this could be membrane stretch or light or something else. Here it is called "displacement"
: the external variable has to be played in via a vector.play in the hoc
: stochastic gating is implemented via Gillespie algorithm,
: performance should be checked when Ntotal - number of channels - increases
: for this two state implementation, 10 to 1000 channels ran as fast as the mean field implementation
: 05/2019 andreas neef
NEURON {
SUFFIX stomet
GLOBAL displacement
NONSPECIFIC_CURRENT icat
RANGE gbar, Ntotal, N_open, updateCount
}
UNITS {
(mV) = (millivolt)
(cm2)= (centimeter2)
(mA) = (milliamp)
(S) = (siemens)
}
PARAMETER {
dt (ms)
v (mV)
rate (1/ms)
time (ms)
updateCount
Ntotal
N_open
closingRate (1/ms)
openingRate (1/ms)
rateCO (1/ms)
rateTOTAL (1/ms)
}
CONSTANT {
: parameters estimated from figures in Beurg et al. 2014 for outer hair cell
: mechanotransduction channels
maxrate = 1.69 (1/ms)
s_half_o= 204 (nm)
s_half_c= 15 (nm)
k_o = 50.6 (1/nm)
k_c = 51.4 (1/nm)
}
ASSIGNED {
displacement (nm)
icat (mA/cm2)
gbar (S/cm2)
}
STATE {C O}
BREAKPOINT {
closingRate = maxrate/(exp( (displacement - s_half_c)/k_c) + 1)
openingRate = maxrate/(exp(-(displacement - s_half_o)/k_o) + 1)
rateCO = C*openingRate
rateTOTAL = rateCO + O*closingRate
SOLVE kinetics METHOD sparse
icat = v * O * gbar / Ntotal
}
INITIAL {
C = Ntotal - N_open
O = N_open
time = 0
closingRate = maxrate/(exp( (displacement - s_half_c)/k_c) + 1)
openingRate = maxrate/(exp(-(displacement - s_half_o)/k_o) + 1)
rateCO= C*openingRate
: the rate with which a transition from C to O occurs:
: occupancy of C times opening rate
rateTOTAL = rateCO + O*closingRate
: the rate at which ANY transition occurs
}
KINETIC kinetics {
~ C <-> O ( 0 , 0) : just here in case the call "SOLVE"
: requires a recipe like this
time = updateTime()
: presumably, most updates end here, because the update time is beyond dt/2
: only for Markov models this is not a problem, because the process has no memory
: one dt/2 later the transition probabilities are unchanged
WHILE ( time < dt/2) : using dt/2 instead of dt, because the block is executed twice per dt update
{
updateStates()
updateRates()
time = time + updateTime()
}
updateCount = updateCount + 1
: checks how often the entire update
: has been run
: just there to debug
: it turns out, the kinetics is call TWICE per dt (as described in the NMODL instructions)
: therefore, the limit in the WHILE clause is now set to dt/2
: in principle the updateCount could now be removed, but it is reassuring to be able to check whether
: a simulation of T/dt time steps produced indeed 2*T/dt update steps
N_open = O
}
PROCEDURE updateStates(){
: this determines (stochastically)
: which transition (out of all possible ones) actually occured
: and updates the states
LOCAL ran
ran=scop_random()
ran=ran*rateTOTAL
IF (ran < rateCO) : an opening event has occured
: in the general case, here not only an upper bound has to be tested, but multiple comparisons with several brakets
: have to be made. Ideally, sort the brakets first, so that the most likely transition is tested first
: the sorting overhead might also kill performance, so test whether its worth it
{
C = C - 1
O = O + 1
}
: this is the place to add more states
ELSE : something has occured, so realize the last option
{
C = C + 1
O = O - 1
}
}
PROCEDURE updateRates() {
rateCO = C*openingRate
rateTOTAL = rateCO + O*closingRate
}
FUNCTION updateTime(){
: this draws the time intervall until the next event
: and returns it
LOCAL ran
ran=scop_random()
updateTime = -log(ran)/rateTOTAL
}
: it should be possible to use the random number from exponential distribution instead
: but I have zero documentation
: best guess is this line:
: updateTime = exprand()/rateTOTAL
```