WATCH accuracy

[stil]

Hello,

I am trying to build a threshold-detection mechanism, with a state reinitialization, embedded in a NMODL.
after trying several things using at_time() and state_discontinuity(,) without success, I realized that the event-delivery system could help me.

For test purposes, i've build a minimal NMODL that should reproduce what i expect (Suppose that XX results from a process, complicated enough so that i cannot determine analytically the events dates):

Code: Select all

NEURON {POINT_PROCESS test}
STATE {XX}
INITIAL {
	XX = 1.0
	net_send(0,1)
}
BREAKPOINT {SOLVE states METHOD cnexp}
DERIVATIVE states{XX' = -XX/10.0}
NET_RECEIVE(w){
	if (flag == 1){
		WATCH(XX<0.1) 2
	}
	if(flag == 2){ : internal event : reinit XX occur now
		XX = 1.0
	}
}

with the following hoc to drive it :

Code: Select all

// PRE PART
create CELL

objref truc
truc = new test(.5)

load_file("nrngui.hoc")
objref cvode
cvode = new CVode(1)
cvode_active(1)

I see that the threshold detection works as it it supposed to, but i was hoping that a mechanism would have determined the exact date of the threshold crossing so that XX would never be < 0.1.
Is there a way to ensure that? Is there something more i should know (dealing with cvode_event ?)

Thanks for your help !

[ted]

As a rule, fixed time step solutions must overshoot threshold crossings except in the rare instance in which the variable becomes equal to the threshold precisely at a solution time. If you insist that some variable never exceed a particular limit, then in the BREAKPOINT block you'll have to test for overshoot and prevent it there, e.g.
if (foo>thresh) {
foo = thresh
}

You might prefer to create a NetCon that monitors the variable, and use the NetCon's record method to force execution of some hoc statement that does what you want.

[stil]

Thanks Ted, i tried your workaround, and it works fine with constant step-size.
But, i would like to use variable stepsize, and it seems that the statements in the BREAKPOINT block were not evaluated. For example, this modification in the code has no effect on the results (i.e i have an overshoot):

Code: Select all

NEURON {POINT_PROCESS test}
STATE {XX}
INITIAL {
	XX = 1.0
	net_send(0,1)
}

BREAKPOINT {
SOLVE states METHOD derivimplicit
if(XX<0.1){
	XX = 0.1
}
}

DERIVATIVE states{
XX' = -XX/10.0
}

NET_RECEIVE(w){
	if (flag == 1){
		WATCH(XX<0.10001) 2
	}
	if(flag == 2){ : internal event : a release must occur now
		net_event(t)
		XX = 1.0
	}
}

I found that it "worked" when i moved the if(XX<0.1){XX = 0.1} in the DERIVATIVE block :

Code: Select all

DERIVATIVE states{
XX' = -XX/10.0
if(XX<0.1){
	XX = 0.1
}
}

In that case, i have a WARNING at compilation-time :
WARNING: XX (a STATE) is assigned a value in a DERIVATIVE block.
Multistep integrators (such as Runge) may not work correctly.

Is it normal that if statements within the BREAKPOINT block are not evaluated ? Did i miss something ?

Thanks again !

[ted]

i would like to use variable stepsize, and it seems that the statements in the BREAKPOINT block were not evaluated. For example, this modification in the code has no effect on the results (i.e i have an overshoot)

Zip up just enough code (hoc, ses, and mod files) for me to reproduce the phenomenon, and I will see what can be done to fix it.

i have a WARNING at compilation-time :
WARNING: XX (a STATE) is assigned a value in a DERIVATIVE block.
Multistep integrators (such as Runge) may not work correctly.

Probably not a good idea to do that, then.

[stil]

Thanks Ted, i sent it to you.

Isn't it related to the fact that the BREAKPOINT block is evaluated twice per timestep with cvode ?
It's just an intuition, i still cannot realize what happens actually ...

[ted]

Isn't it related to the fact that the BREAKPOINT block is evaluated twice per timestep with cvode ?

No (BREAKPOINT is evaluated more than twice per advance with cvode) and no (it is not related to the number of times BREAKPOINT code is evaluated). I'll look at your code and get back to you. Traveling today so there will be some delay.

[ted]

So there are three issues to be discussed.

First, WATCH is perfectly accurate. It is the fact of doing numerical integration with a digital computer that forces the discretization of time that allows the watched variable to overshoot a threshold.

Second, regarding assignment to state variables in a BREAKPOINT block--the warning message didn't say it couldn't be done. It said that the result could be failure of a multistep integrator to work properly. Multistep integrators use information about the temporal evolution of state variables in the past in order to make a better prediction of their future values. This requires the variables and their first few derivatives to be continuous (the derivative order at which discontinuity is allowable depends on the details of the integration method). It was naive of me to suggest preventing overshoot by forcing XX to a particular value. That will work with NEURON's default integrator (implicit Euler) but it may not work with adaptive integration.

Third, the final issue is what to do about the threshold overshoot. If you have a mechanism with dynamics so simple that there is an analytical solution--as is the case for the example you sent--then given the current state of the mechanism at any time, its future behavior is completely predictable, and you can discover the time at which it will intersect threshold without having to do numerical integration.* All computation can then be done inside the NET_RECEIVE block and there is no need for a BREAKPOINT block at all. And if you use the ARTIFICIAL_CELL keword in the NEURON block, e.g.

Code: Select all

NEURON {
  ARTIFICIAL_CELL Test
  RANGE tau
}

and you can treat it the same way you would an IntFire or a NetStim. I didn't see anything in the example you sent that indicated that it has to be a point process.

For that matter, if all you want is a sequence of events at predetermined intervals, why not use a NetStim? The NetCon class's record() method allows events to trigger execution of arbitrary hoc statements. Or you could just use a combination of an FInitializeHandler and cvode.event to force execution of hoc statements at regular intervals.

*--depending on the details of the analytical solution, you might have to resort to Newton's method to discover when the solution crosses threshold, but that can be done to any reasonable degree of precision inside a NET_RECEIVE block much more quickly than the time required to integrate over a series of time steps. Granted, in this case your state variable will still over- or undershoot the threshold by a small amount, but you can decide how small the error is.