This is a topic I take to the surface from time to time. I want to simulate stochastic single channels from nmodl code instead of using the channel builder and the associated hoc code (I'm not going to discuss the reasons of my choice here).
I have successfuly written the mod file for this simulation using the 'classic' algorithm: at each time step a uniform random number is drawn for each of the N channels in the section, to decide wether the channel will change its state. This is very time consuming and with several hunderds of channels things can get very slow.
Now I'm trying to do it with the optimized algorithm that the Channel Builder uses: first, a random number will give the time of the next transition (and until then nothing is calculated!) and then a second random number will decide which kind of transition will occur.
This is the best that I have got. Notice that I'm using a distributed mechanism instead of a point process. At the beggining, the number of channels (N) is calculated based on the maximal conductance, the unit conductance and the area.
Code: Select all
TITLE Mock stochastic voltage-dependent sodium channel
:
: Written by Patricio Orio, Ago 2009
:
NEURON {
SUFFIX stNav
USEION na READ ena WRITE ina
RANGE gmax, N, No, Po, vlast, prev_ev, next_ev
}
UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
}
PARAMETER {
gmax = 0.00025 (mho/cm2)
V0 = -40 (mV)
z = 0.09 (/mV)
tau = 10 (ms)
gu = 20e-12 (mho)
}
ASSIGNED {
area (micron2)
ina (mA/cm2)
v (mV)
ena (mV)
i (mA/cm2)
N
No
Po
dt (ms)
alpha (/ms)
beta (/ms)
prev_ev (ms)
next_ev (ms)
vlast (mV)
}
INITIAL {
Po = 1/(1+exp(-z*(v - V0)))
alpha = Po/tau
beta = (1-Po)/tau
N = floor((1e-8) * gmax * area / gu + 0.5)
No = floor(Po*N + 0.5)
prev_ev = 0
next_ev= -log(1-scop_random())/(No*beta + (N-No)*alpha)
}
BREAKPOINT {
at_time(next_ev)
if ( fabs(v-vlast) > 5) { : if the voltage changes
prev_ev = t : more than 5 mV
recalc_ev() : alpha and beta have to be
} : recalculated
if (t >= next_ev) {
transitions()
}
ina = gmax * (v-ena) * No/N
}
PROCEDURE transitions() {
LOCAL delta
delta = -1
if (scop_random() >= No*beta/((N-No)*alpha+No*beta)) {delta = 1}
No = No + delta
prev_ev = next_ev
recalc_ev()
}
PROCEDURE recalc_ev() {
Po = 1/(1+exp(-z*(v - V0)))
alpha = Po/tau
beta = (1-Po)/tau
next_ev = prev_ev - log(1-scop_random())/(No*beta+(N-No)*alpha)
vlast = v
}
This code works perfectly (or at least it seems) with constant time step. However, when variable time step is enabled, it crashes with the message:
If either the at_time() statement OR the if block just after it is removed, neuron doesn't crash. But as you can imagine both things are needed for some reason.Failed assertion
te > tstop_ || te <= t0_
at line 615 of file /home/Hines..../cvodeobj.cpp
Why is this crash occurring? and why is that it only occurs when both the at_time and the if block are present?
I have read that at_time is deprecated, that it should be replaced with event handling instructions; however all I have found about events is related to network simulations and this is not the case.
Any help will be appreciated.
Regards