WATCH and net_send

[sec6]

Is there any known incompatability between the WATCH statement and the net_send function?

I have a NET_RECEIVE block containing both a WATCH and several net_send calls. The WATCH fails to to put an event in the queue when the parameter it's WATCH-ing crosses threshold. This behavior ceases when I comment out the net_send calls. In fact, I can make the behavior come & go by changing the value of the "delay" argument to net_send.

However, I'm having the devil of a time producing a short piece of code that reliably reproduces the problem: there's something stochastic about it. For example, I reliably make the bug appear & dispear, repeatedly, by varying the delay parameter between two values, then, when I think I've gotten rid of it for good, I move on & start editing some other part of the code, and the problem reappears. (Obviously, I'm recompiling after each edit to the sourcecode.) I suspect I'm dealing with a race condition.

With a problem like this, it's not good enough to make the problem go away: I have to understand what's causing it. Otherwise, Murphy's law garantees it'll reappear mysteriously after the code's passed its unit tests and is in production.

I've written a pair of .mod and .hoc files which exercise the mechanism repeatedly, with a range of delay parameters, until the WATCH fails. This *seems* to replicate the problem reliably.
There's a range variable "watch_time" which takes value of 0 or nonzero. If it's nonzero, the mechanism looks for a threshold-crossing of the "t" variable (time). This does *not* bring out the bug. If "watch_time" is set to zero, the mechanism looks for a threshold crossing of the "v" variable (membrane potential). This *does* bring out the problem.

Here's the .mod file ("ex_watch.mod"):

Code: Select all

NEURON {
	POINT_PROCESS ex_watch
	RANGE out_var,watch_time,delay
}

PARAMETER {
	t_transition = 50 (ms)
	v_transition = 0 (mV)
}

ASSIGNED {
	watch_time
	out_var
	delay
	v
	:t :Do I really need to declare t as an ASSIGNED variable?  Apparently not. (Built-in NMODL global)
}

STATE {
	dummy_state :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

INITIAL {
	out_var = 0
	dummy_state = 0
	net_send(0,99)
}

BREAKPOINT {
	SOLVE states METHOD cnexp :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

DERIVATIVE states {
	dummy_state' = 0 :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

NET_RECEIVE (dummy_weight) {
	if (watch_time == 1) {
		WATCH (t > t_transition) 1
	} ELSE {
		WATCH (v > v_transition) 1
	}
	IF (flag == 1) {
		out_var = out_var + 1
	} ELSE {
		IF (flag == 0) {
			net_send(delay,-1)
		} ELSE {
			IF (flag == -1) {
				:do nothing
			} ELSE {
				IF (flag == 99) {
					:do nothing
				} ELSE {
					:This should never happen.
				}
			}
		}
	}
}

Here's the .hoc file to exercise it (name is "ex_ex_watch.hoc"):

Code: Select all

//BEGIN Define simpleCell class

begintemplate simpleCell
public celldef, position
public soma

proc celldef() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
}

create soma

proc topol() { local i
  basic_shape()
}
proc basic_shape() {
  soma {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(15, 0, 0, 1)}
}

proc position () {
	soma {
		netX = $1
		netY = $2
		netZ = $3
	}
}

objref all
proc subsets() { local i
  objref all
  all = new SectionList()
    soma all.append()

}
proc geom() {
  soma {  L = 30  diam = 30  }
}
proc geom_nseg() {
  soma area(.5) // make sure diam reflects 3d points
}
proc biophys() {
  soma {
    Ra = 80
    cm = 1
    insert hh
      gnabar_hh = 0.12
      gkbar_hh = 0.036
      gl_hh = 0.0003
      el_hh = -54.3
  }
}
access soma
endtemplate simpleCell

//END Define simpleCell class

//Create first cell
objref cell1
cell1 = new simpleCell()
cell1.celldef()

//Create second cell
objref cell2
cell2 = new simpleCell()
cell2.celldef()

//Insert ex_watch mechanism into cell1
objref mech1
cell1.soma mech1 = new ex_watch(0.5)

//Instrument cell2 with a current-injecting electrode
objref trode
cell2.soma trode = new Ipulse2(0.5)
trode.del = 50
trode.amp = 0.5
trode.dur = 1.0
trode.num = 1
trode.per = 50

//Connect the mechanism in cell1 to cell2 via NetCon object
objref nc
cell2.soma nc = new NetCon( &v(0.5), mech1, 20, 0, 1)  //threshold 20, delay 0, weight 1

//Load the standard run system
load_file("stdrun.hoc")

tstop = trode.del * 2

//mech1.watch_time = 1
mech1.watch_time = 0

proc testit() {
incr = 0.1
final = 25
	for i = 1, (final/incr) {
		mech1.delay = i*incr
		run()
		if (mech1.out_var == 0) {
			print "WATCH failed at the following delay:"
			print mech1.delay
			break
		}
	}
}

////Create graphs to monitor range variables 
objref g1
g1 = new Graph()
addplot(g1,0)
g1.size(0,tstop,0,5)
g1.addvar("mech1.out_var")

//Create graphs to monitor range variables 
objref g2
g2 = new Graph()
addplot(g2,0)
g2.size(0,tstop,-70,50)
g2.addvar("cell2.soma.v(0.5)")

Here's a transcript of an interactive hoc session demonstrating this:

Code: Select all

	
[58] ./i686/special
NEURON -- Release 6.0.4 (1819) 2007-07-20
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2007
See http://www.neuron.yale.edu/credits.html

loading membrane mechanisms from /Users/sec/neuron/current/i686/.libs/libnrnmech.so
Additional mechanisms from files
 ipulse2.mod ex_watch.mod
oc>load_file("ex_ex_watch.hoc")
	0 
	0 
	1 
	1 
	1 
	1 
	1 
	1 
oc>mech1.watch_time=0
oc>testit()
WATCH failed at the following delay:
0.1 
oc>mech1.watch_time=1
oc>testit()
oc>

The .hoc file produces two plots. It'll run quicker without them, but they help reassure me that I haven't introduced a new bug in my testing code.

When all is said & done, despite their frustrations, NEURON & NMODL are the greatest thing since sliced bread.

[sec6]

I just discovered that

A WATCH statement only needs to be executed once, during simulation initialization.
http://www.neuron.yale.edu/phpBB2/viewt ... 16cf9a54da

However, this is not the source of my problem, because the following, revised version of ex_watch.mod exhibits exactly the same behavior.

Code: Select all

NEURON {
	POINT_PROCESS ex_watch
	RANGE out_var,watch_time,delay
}

PARAMETER {
	t_transition = 50 (ms)
	v_transition = 0 (mV)
}

ASSIGNED {
	watch_time
	out_var
	delay
	v
	:t :Do I really need to declare t as an ASSIGNED variable?  Apparently not. (Built-in NMODL global)
}

STATE {
	dummy_state :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

INITIAL {
	out_var = 0
	dummy_state = 0
	net_send(0,99)
}

BREAKPOINT {
	SOLVE states METHOD cnexp :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

DERIVATIVE states {
	dummy_state' = 0 :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?
}

NET_RECEIVE (dummy_weight) {
	IF (flag == 1) {
		out_var = out_var + 1
	} ELSE {
		IF (flag == 0) {
			net_send(delay,-1)
		} ELSE {
			IF (flag == -1) {
				:do nothing
			} ELSE {
				IF (flag == 99) {
					IF (watch_time == 1) {
						WATCH (t > t_transition) 1
					} ELSE {
						WATCH (v > v_transition) 1
					}
				} ELSE {
					:This should never happen.
				}
			}
		}
	}
}

[ted]

Quite apart from whether there is a bug, or what its cause and fix may be, what are you
trying to accomplish? i.e. what is the underlying purpose?

[sec6]

Shouval et al '02 PNAS 99(16):10831

Synaptic plasticity model based on the notion that small rise in intraspinous Ca depresses, and large rise potentiates (e.g. D'Alcantara et al '03 Eur J Neurosci 17(12):2521, though the model's at a higher level of abstraction than that).

It relies, for detection of pre-post spike coincidence on temporal overlap between glutamate binding to NMDA receptors and a backpropagating dendritic afterpotential.

The WATCH statement, plus one of the net_send calls initiates the afterpotential at a specified delay following a somatic spike in the neuron in which the mechanism is inserted.

The number of NMDA receptors in a single spine is actually rather small (as you probably know, 'though I didn't until recently). Shouval & Kalantzis '05 J Neurophys 93(2):1069 showed that the resulting stochastic variation in NMDA-R-mediated Ca current results in *qualitative* changes in the behavior of their model (dissapearance of long-latency pre-post LTD).

For the deterministic version of the model, I've managed to get agreement between output of Dr. Shouval's C & Matlab code (a copy of which he generously provided me) and my own implementation in NEURON & NMODL. I'm now introducing random variation in the Ca current, and that's where I'm encountering difficulty.

Shouval & Kalantzis implement that variation in a way that won't serve my purposes. I'm trying to implement it by modelling NMDA receptor closing* as a Markov process, using net_send, to put channel closing events in the event queue at random intervals drawn from an appropriate distribution. This seems to work OK (that is, if I set the number of channels to be very large, I get a calcium-current vs. time curve identical to the detrministic model (which, in turn, is identical to my gold standard, the original Matlab/C implementation by Shouval.)) However, those net_send calls seem to be messing up the WATCH (or else I'm totally misunderstanding something -- certaily a very real possibility.)

[*Just closing, not opening, and only a three-state model, because that's what corresponds to Shouval's model. Later on, I might try to implement something more realistic.]

Sorry for omitting all this from the original post. I actually thought I'd make the question easier to understand by paring away extraneous detail.

[ted]

:t :Do I really need to declare t as an ASSIGNED variable?

t does not have to be declared.

dummy_state :Can I perhaps omit the STATE, BREAKPOINT, and DERIVATIVE blocks altogether?

No need for them since the STATE variable plays no real role in the mechanism.

[ted]

Shouval & Kalantzis implement that variation in a way that won't serve my purposes. I'm trying to implement it by modelling NMDA receptor closing* as a Markov process, using net_send, to put channel closing events in the event queue at random intervals drawn from an appropriate distribution. This seems to work OK (that is, if I set the number of channels to be very large, I get a calcium-current vs. time curve identical to the detrministic model (which, in turn, is identical to my gold standard, the original Matlab/C implementation by Shouval.)) However, those net_send calls seem to be messing up the WATCH (or else I'm totally misunderstanding something -- certaily a very real possibility.)

[*Just closing, not opening, and only a three-state model, because that's what corresponds to Shouval's model. Later on, I might try to implement something more realistic.]

Ligand-gated point processes with stochastic kinetics described by kinetic schemes--
sounds like something the Channel Builder could do easily. The only aspect that might
be difficult is emulating voltage dependent Mg block.

[hines]

I'll look into this. But the simulation needs another mod file that implements Ipulse2.
Can you please send me a zip file with
everything needed to build the model to
michael dot hines at yale dot edu

[sec6]

Ligand-gated point processes with stochastic kinetics described by kinetic schemes--
sounds like something the Channel Builder could do easily. The only aspect that might
be difficult is emulating voltage dependent Mg block.[/quote]

Hmm, looks promising.
I'll start reading documentation/tutorials on the channel builder.
Meanwhile:

The channel builder has a GUI, right? So far, I've managed to avoid dealing with the NEURON GUI (I really prefer editing hoc & NMODL code directly). Am I obliged in this case to use the GUI? If I understand correctly, the channel builder produces NMODL sourcecode. Can I build my channel with the channel builder, then open up the output in a text editor? Will it be safe to alter the resulting NMODL sourcecode, or must I leave that untouched?


Thanks.

[sec6]

I'll look into this. But the simulation needs another mod file that implements Ipulse2.
Can you please send me a zip file with
everything needed to build the model to
michael dot hines at yale dot edu

Done.

Thanks.

[ted]

The channel builder has a GUI, right? So far, I've managed to avoid dealing with the NEURON GUI (I really prefer editing hoc & NMODL code directly).

"The proper tool for the proper use." Some things are much more easily done with the GUI,
some with hoc, and some are best done by combining the two. Ignore one or the other, and
you're working with a self-imposed handicap.

If I understand correctly, the channel builder produces NMODL sourcecode.

Not so.
If you read that somewhere on the Forum or NEURON's other WWW pages, please let me
know so I can correct whatever text is conveying that misunderstanding.

My speculation about whether your goal could be achieved with the Channel Builder remains
for the moment speculation. The implementation would have the following components:
--a point process with a NET_RECEIVE block would receive "synaptic activation events," and
respond by producing a transient change of a variable that represents concentration of a
"transmitter substance."
--a Channel Builder that uses a kinetic scheme to specify the open/closed state of a
ligand-gated channel, operating in "single channel mode" so that it emulates stochastic
gating. The value of the open state variable is visible to hoc.
--a point process that uses a POINTER to get the value of the Channel Builder's open
state variable, and multiplies that by a function that expresses the voltage dependence
of NMDA conductance. This doesn't have to be a second point process--the necessary
code could easily be part of the point process mentioned above

The only limitation I see is that a single point process - Channel Builder pair could not
serve as the recipient of multiple afferent streams, each with its own weight and activation
history. This is a problem only if you need to attach more than one of these "synapses"
to a single compartment (segment).

[sec6]

The business about the channel builder producing NMODL sourcecode was the result of a quick Google. If I find it again I'll make a note of the url. It's possible, though, that I just mis-read.

The implementation would have the following components:
--a point process with a NET_RECEIVE block would receive "synaptic activation events," and respond by producing a transient change of a variable that represents concentration of a "transmitter substance."
--a Channel Builder that uses a kinetic scheme to specify the open/closed state of a ligand-gated channel, operating in "single channel mode" so that it emulates stochastic gating. The value of the open state variable is visible to hoc.
--a point process that uses a POINTER to get the value of the Channel Builder's open state variable, and multiplies that by a function that expresses the voltage dependence of NMDA conductance. This doesn't have to be a second point process--the necessary code could easily be part of the point process mentioned above.

I might have to do it that way. However, that approach requires hoc code to "glue" together the channel builder and the point process. I was hoping to keep this self-contained so that I could simply substitute it into an existing network model in place of the current, non-plastic glutamate synapse mechanism.

The only limitation I see is that a single point process - Channel Builder pair could not serve as the recipient of multiple afferent streams, each with its own weight and activation history. This is a problem only if you need to attach more than one of these "synapses" to a single compartment (segment).

I do need to do that.
However, could I not have multiple point-process / Channel Builder pairs inserted into the same segment?

[ted]

I might have to do it that way. However, that approach requires hoc code to "glue" together the channel builder and the point process. I was hoping to keep this self-contained so that I could simply substitute it into an existing network model in place of the current, non-plastic glutamate synapse mechanism.

It might be possible to implement modularity that allows you to use the same basic code
framework, with "pluggable" substitution of such details, by a minor refactoring of your
program. Specifically, the bit of code that defines and places synaptic mechanisms
could be wrapped in procs that are factored out into a separate hoc file. You could
have several different versions of that hoc file, each of which has different internal
details, but a similar interface (identical proc name, similar or identical args). It might
also be helpful to define a new object class that pulls all these things together. Then
just by changing the arguments to one or more load_file() statements, you can bring
up models with the same general architecture, but different synaptic details. For an
example of something like this, see the NEURON implementation of test network models
used for
Brette R, Rudolph M, Carnevale T, Hines M, Beeman D, Bower JM, Diesmann M, Morrison A, et al. (2007)
Simulation of networks of spiking neurons:
A review of tools and strategies. J Comp Neurosci 23:349-98
http://senselab.med.yale.edu/modeldb/Sh ... rks\NEURON\

The only limitation I see is that a single point process - Channel Builder pair could not serve as the recipient of multiple afferent streams, each with its own weight and activation history. This is a problem only if you need to attach more than one of these "synapses" to a single compartment (segment).

I do need to do that.
However, could I not have multiple point-process / Channel Builder pairs inserted into the same segment?

Not with the strategy I had in mind. The point process that receives synaptic events
has to compute transmitter concentration, and the Channel Builder has to be able to
access this value. My idea was to use the extracellular concentration of a bogus "ionic
species" to represent transmitter concentration, because the Channel Builder can "read"
ion concentrations. A trivial workaround is to use an electrically insignificant passive
compartment as the "host" for each synapse that must be attached to the same location.
That adds an extra ODE per synapse. But maybe there's some other way to accomplish
this communication between mechanisms.

[sec6]

Thanks. Ultimately, I was able to "solve" the problem at the level of the NMODL code. (I missed the email notifying me you'd replied to my last post, and noticed it when I came here to post another question.)
I'll summarize for the benefit of anybody who encounters a similar problem and finds his/her way to this thread.

Michael Hines kindly agreed to look at my code, if I could write a minimal program that exhibited the problem. I wrote a "toy" NMODL program to exercise WATCH and net_send -- but it didn't exhibit the bug. I started adding in, bit by bit, the functionality of my original problem program, looking for the bit that would elicit the bug. Eventually, I'd rewritten the entire original program (but with "cleaner" syntax) and it worked buglessly. I'd be happier if I actually knew what I'd fixed, still: you can't always get what you want, but if you try sometimes, you get what you need.

[ted]

Thanks for describing the path that culminated in working code. That's the very approach
that I try to encourage people to use, because it so often leads to success.