problems with cvode, at_time and something else

NMODL and the Channel Builder.
Post Reply
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

problems with cvode, at_time and something else

Post by patoorio »

Hi,

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
}
The 'if' block after at_time is in order to force the recalculation of alpha and beta (and therefore the next event time) if there is an abrupt change of voltage. I couldn't figure out a better way to do it, but it must be done in some way.
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:
Failed assertion
te > tstop_ || te <= t0_
at line 615 of file /home/Hines..../cvodeobj.cpp
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.
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
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: problems with cvode, at_time and something else

Post by hines »

The use of at_time(tev) requires that tev not be inside the current integration interval with the variable time step method. This is not guaranteed when the "if" statement is present.

Events are better handled by a NET_RECEIVE block (you can use the WATCH statement to watch for trigger conditions) but NET_RECEIVE is restricted to POINT_PROCESS or ARTIFICIAL_CELL and cannot be used with a SUFFIX (density) model.
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

Thanks a lot for the answer.
It's not a big deal to make it a point process, but it's been hard to find good documentation about NET_RECEIVE blocks and WATCH (by 'good' documentation I mean an explanation of the commands in the way it is found for hoc commands in NEURON's help files).
Where can I find a good example of usage that may be useful in this context?
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: problems with cvode, at_time and something else

Post by ted »

NET_RECEIVE is discussed in chapter 10 of The NEURON Book, and in
Hines, M.L. and Carnevale, N.T. Discrete event simulation in the NEURON environment. Neurocomputing 58-60:1117-1122, 2004
which is available as a preprint from the WWW site's "Papers about NEURON" page
http://www.neuron.yale.edu/neuron/bib/nrnpubs.html
For examples of working code, see the source code for ExpSyn, Exp2Syn, IntFire2, IntFire2, and IntFire4, which are in nrn-n.n/src/nrnoc of the expanded nrn-n.n.tar.gz, or c:\nrnnn\src\nrnoc if you're using MSWin. Also see these NEURON implementations of models of short-term use-dependent synapatic plasticity:
http://senselab.med.yale.edu/modeldb/Sh ... model=3264
http://senselab.med.yale.edu/modeldb/Sh ... model=3815

There is some discussion of WATCH, and a link to an example of STDP that uses it, in this "What's new" http://www.neuron.yale.edu/neuron/news/ ... 23_06.html
WATCH (var > thresh) flagvalue
is used in a NET_RECEIVE block to specify a condition in the postsynaptic cell that will generate a self-event with latency 0 and a specified flag value. Generally, WATCH is used to make NEURON monitor a variable for threshold crossing, and generates a self event with the specified flag value when the threshold is crossed. If the postsynaptic cell is a biophysical model cell, var is usually local membrane potential (or cai or some other concentration); if the postsynaptic cell is an artificial spiking cell, var is one of that cell's state variables. But WATCH could in principle be anything, such as the total number of spikes that a cell has fired, or perhaps even t (time) (although I haven't seen it used to monitor time).
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

Thanks!
I'll try to figure it out and if I succeed I'll post the code.
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

Hi,
This is my current iteration:

Code: Select all

TITLE Mock stochastic voltage-dependent sodium channel
: 
: Written by Patricio Orio, Ago 2009
:

NEURON {
	POINT_PROCESS stNav2
	USEION na READ ena WRITE ina
	RANGE N, No, Po, ina, next_ev, vlast
}

UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
}

PARAMETER {
	N  = 100
	gu = 20e-6	(micromho)
	V0 = -40	(mV)
	z = 0.09	(/mV)
	tau = 10	(ms)
}

ASSIGNED {
	area	(micron2)
	ina		(nanoamp)
	v		(mV)
	ena		(mV)
	No
	Po
	dt		(ms)
	alpha	(/ms)
	beta	(/ms)
	next_ev	(ms)
	vlast	(mV)
}

INITIAL {
	Po = 1/(1+exp(-z*(v - V0)))
	alpha = Po/tau
	beta = (1-Po)/tau
	No = floor(Po*N + 0.5)
	net_send(0,2)
}

BREAKPOINT {
	ina = gu * No * (v-ena)
}

NET_RECEIVE (w) {
	if (flag == 1 ) { 
		if (scop_random() >= No*beta/((N-No)*alpha+No*beta)) {No = No + 1} else {No = No -1}
		Po = 1/(1+exp(-z*(v - V0)))
		alpha = Po/tau
		beta = (1-Po)/tau
		next_ev = - log(1-scop_random())/(No*beta+(N-No)*alpha)
		net_send(next_ev,1)
	}
	if (flag == 2) {
		Po = 1/(1+exp(-z*(v - V0)))
		alpha = Po/tau
		beta = (1-Po)/tau
		next_ev = - log(1-scop_random())/(No*beta+(N-No)*alpha)
		net_send(next_ev,1)
		vlast = v
		WATCH (fabs(v-vlast) > 5) 2
	}
}
Basically it works, but there is a problem when using variable dt: the transition rate gets much higher than expected. It also happens with constant dt when selecting a dt 100 times smaller. None of that happens with the previous mechanism that uses at_time(). It also happens regardless of the value of tau.
The problem has to be with WATCH, because if I comment that line (therefore an event with flag 2 is generated only at the beginning) the mechanism becomes 'dt-independent' as it should be. I have tried to debug what may be going on, but I don't know how to monitor events during runtime. However I discovered something: the 'single step' button of the RunControl panel always advances 0.025 ms, regardless of the actual dt value.
If you want to try this mechanism, simply put it inside a compartment and add an SEClamp. Then give some depolarizing pulse and plot No; with the current parameters I suggest something like -20 mV for 100 ms.
Thanks for the help!

PS: I've tested it with both NEURON 6.2 and 7.0
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: problems with cvode, at_time and something else

Post by hines »

I believe the problem is that every time the voltage changes from vlast by more than 5mV you start a whole
new series of flag=1 streams that never ends. This is occurring with the fixed step method as well but
you might not notice except for going slower and slower as time increases.
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

hines wrote:I believe the problem is that every time the voltage changes from vlast by more than 5mV you start a whole
new series of flag=1 streams that never ends.
And why is that? I can't see it in the code. Did I miss something in the way WATCH works?
hines wrote:This is occurring with the fixed step method as well but
you might not notice except for going slower and slower as time increases.
Yes indeed!, however it only happens when the WATCH statement is there
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: problems with cvode, at_time and something else

Post by hines »

A flag=2 event indicates (from the watch statement) that the voltage has changed
by 5mV from the previous flag=2 event.
A flag=2 event sends a flag=1 event. And a flag=1 event sends a flag=1 event.

The number of flag=1 events on the queue is therefore equal to the total number of
flag=2 events you have received.
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

hines wrote:A flag=2 event sends a flag=1 event. And a flag=1 event sends a flag=1 event.
Right, but all events are sent with a proper next_ev time delay, isn't it? This is what seems not to be working when WATCH is present.
Is there a way to monitor event generation, delays, and arrivals while the simulation is running? (the way I can do with assigned variables, for instance)
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: problems with cvode, at_time and something else

Post by hines »

You can add printf statements to the mod file like
printf("%d %g v=%g next=%g\n", flag, t, v, next_ev)
or at the hoc level
cvode.debug_event(1)

However, to a first approximation, I don't see why you are sending a flag=1 event from a flag=2 event
(modulo initialization which perhaps should be a flag=3 event). Shouldn't a flag=2 event only change
parameters and not calcuate a next_ev time. And if you are worried about not getting the effect of
changed parameters til after the next event (not a large numerical problem if many events per 5mv
interval) you can experiment with net_move(newtime) which is limited to one selfevent being on the queue
at a time.
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

Great!
I don't know why I assumed that the previous flag=1 event was deleted or something by the new one. So now I understand that I was creating and creating events at each voltage change.
The fixed (and probably final) code looks like this:

Code: Select all

TITLE Mock stochastic voltage-dependent sodium channel
: 
: Written by Patricio Orio, Ago 2009
:

NEURON {
	POINT_PROCESS stNav2
	USEION na READ ena WRITE ina
	RANGE N, No, Po, ina, next_ev, vlast
}

UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
}

PARAMETER {
	N  = 100
	gu = 20e-6	(micromho)
	V0 = -40	(mV)
	z = 0.09	(/mV)
	tau = 10	(ms)
}

ASSIGNED {
	area	(micron2)
	ina		(nanoamp)
	v		(mV)
	ena		(mV)
	No
	Po
	dt		(ms)
	alpha	(/ms)
	beta	(/ms)
	next_ev	(ms)
	vlast	(mV)
}

INITIAL {
	Po = 1/(1+exp(-z*(v - V0)))
	alpha = Po/tau
	beta = (1-Po)/tau
	No = floor(Po*N + 0.5)
	net_send(0,0)
}

BREAKPOINT {
	ina = gu * No * (v-ena)
}

NET_RECEIVE (w) {
	Po = 1/(1+exp(-z*(v - V0)))
	alpha = Po/tau
	beta = (1-Po)/tau
	next_ev = - log(1-scop_random())/(No*beta+(N-No)*alpha)
	if (flag == 0 ) {
		net_send(next_ev,1)
		vlast = v
		WATCH (fabs(v-vlast) > 5) 2
	} 
	if (flag == 1 ) {
		if (scop_random() >= No*beta/((N-No)*alpha+No*beta)) {No=No+1} else {No=No-1}
		net_send(next_ev,1)
	}
	if (flag == 2) {
		net_move(t + next_ev)
		vlast = v
		WATCH (fabs(v-vlast) > 5) 2
	}
}
Apparently works perfectly, with both variable and fixed dt. The value for me is that now I can modify alpha and beta (or Po and tau) to make them dependent on any variable I want, or having any function I can imagine.
I'll appreciate any suggestions about testing the mechanism.

Thanks a lot Mike for your patience.
Regards.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: problems with cvode, at_time and something else

Post by hines »

Yes. I believe that is conceptually correct. Though it does not matter here, I generally avoid the use of flag=0 for initialization
since external events have a flag=0. Just don't connect a NetCon to this. Or else if a flag=0 comes from the outside ignore
it with something like
if (flag == 0) {
:ignore
}else{
:self event handlers
}
}
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: problems with cvode, at_time and something else

Post by patoorio »

Good advice. Thanks again.
Post Reply