Event delivery (CVode.event) not impacting simulation time

Anything that doesn't fit elsewhere.
Post Reply
ngreiner
Posts: 11
Joined: Tue Feb 28, 2017 8:12 am

Event delivery (CVode.event) not impacting simulation time

Post by ngreiner »

Dear all,

I am encountering an issue with the use of custom event deliveries using the CVode.event method.
My problem, reproducible with the MWE below, is that the delivery of events does not impact the progress of the simulation time.
Specifically, when I use the variable time step method with maxstep=0.4, and try to deliver events at t=5.0 and t=5.2, these events seem to effectively be delivered (because they are printing stuff in the console), but when I inspect the recorded time vector of the simulation, I see that time went from t=4.84113817206 directly to t=5.24113817206.

The code below

Code: Select all

from neuron import h
h.load_file('stdrun.hoc')

# Create section
s = h.Section()
s.insert('hh')
# s.insert('extracellular')

# Record time
t = h.Vector()
t.record(h._ref_t)

# Simulation parameters
h.tstop = 10.0
h.cvode_active(1)
h.cvode.maxstep(0.4)
h.finitialize()

# Prepare and push events in the queue
def switch_on_eext():
	print 'switching on'
	print h.t
	print ''
	# s.e_extracellular = 1.0

def switch_off_eext():
	print 'switching off'
	print h.t
	print ''
	# s.e_extracellular = 0.0

h.cvode.event(5.0, switch_on_eext)
h.cvode.event(5.2, switch_off_eext)

# Run simulation
h.cvode.solve(h.tstop)

# Print time steps
for ti in t.to_python():
	print ti
produces the following output on my machine:
NEURON -- VERSION 7.5 master (c693a84) 2017-12-07
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2016
See http://neuron.yale.edu/neuron/credits

switching on
5.0

switching off
5.2

0.0
0.220569086028
0.441138172057
0.841138172057
1.24113817206
1.64113817206
2.04113817206
2.44113817206
2.84113817206
3.24113817206
3.64113817206
4.04113817206
4.44113817206
4.84113817206
5.24113817206
5.64113817206
6.04113817206
6.44113817206
6.84113817206
7.24113817206
7.64113817206
8.04113817206
8.44113817206
8.84113817206
9.24113817206
9.64113817206
10.0
It is problematic because I would like to use these events to turn on/off the extracellular batteries of my model, which constitute abrupt changes requiring the time step to get adapted (reduced) locally around these events.

This problem can be circumvented by the use of FInitializeHandler objects, but this solution does not allow to dynamically create new events in the course of a simulation.

I suppose there is something wrong in my code which leads to the observed results, but I don’t know what it is.

If you have a suggestion, it is more than welcome.

With regards,

Nathan
ngreiner
Posts: 11
Joined: Tue Feb 28, 2017 8:12 am

Re: Event delivery (CVode.event) not impacting simulation time

Post by ngreiner »

Hi again,

After reading this post: viewtopic.php?t=3034, I found the means to force the integrator to stop at the times of the specified events by inserting a cvode.re_init() call in the py_callable called by my events.

Code: Select all

def switch_on_eext():
	print 'switching on'
	print h.t
	print ''
	s.e_extracellular = 1.0
	h.cvode.re_init()
However, I now have an issue, which is also reported in the mentioned post (viewtopic.php?t=3034), whereby the integrator prints warnings to the console as:
IDA initialization failure, weighted norm of residual=1.31588
I am not sure, but I believe this is because the changes (in my case in e_extraceullar) occurring during the delivered events are too large and abrupt.

If this is the case, do you know a way to avoid this issue by more appropriately implementing the switching on and off of an external stimulation pulse ?

I think that vectors specifying ramp-up and ramp-down of e_extracellular could be used, but:
(i) the ramping-up or down would be arbitrary, so how to choose it ?
(ii) it seems a more cumbersome solution than the delivery of events at specified times
(iii) would it allow to dynamically configure new extracellular stimulation pulses during simulation ?

Thank you in advance for your help,

Nathan
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Event delivery (CVode.event) not impacting simulation time

Post by ted »

Your intuition that you probably need to implement finite rise and fall times is probably correct. Don't be too concerned about the exact details of the transients, in particular their duration--anything that is "very short" compared to the smaller of [duration of an action potential | duration of your stimulus pulse] will do. What's "very short"? Two orders of magnitude, more or less. Vector play with interpolation is the easiest way to do this. Read about it in the Programmer's Reference documentation of NEURON's Vector class. You only need one t, v data pair per break point in your driving function.

Examples:

This sequence of tuples
0, 0
1, 0
1, 1
2, 1
2, 0
3, 0
specifies a square pulse with amplitude 1 that starts at 1 ms and ends at 2 ms.

This sequence
0, 0
0.99, 0
1.01, 1
1.99, 1
2.01, 0
3, 0
specifies a trapezoidal pulse (rise and fall times are both 0.01 ms) with amplitude 1 that starts at 0.99 ms and ends at 2.01 ms.

Note that each adjacent pair of data points specifies a line segment that defines the slope of the function between those two points. Also note that the slope of the final straight line segment becomes the slope of the function once t passes the final data point. That's why both of these example sequences contains that final data point 3, 0. Omit that, and Vector play with interpolation will happily extrapolate the falling phase of your current pulse at a slope of -1 per 0.02 ms forever so that it soon becomes a ridiculously large negative number.

"I need to drive e_extracellular in multiple segments, each with a different pulse amplitude. Do I need a separate pair of Vectors for each segment?"

Only if the normalized e_extracellular waveforms are different, as would be the case if you cannot assume that the extracellular medium is purely resistive (i.e. can't assume that it is non-dispersive). If you don't have to take dispersion into account, you can use the xtra mechanism to couple a single waveform generating function to all of the e_extracellular instances in your model--see http://www.neuron.yale.edu/ftp/ted/neur ... nd_rec.zip.
ngreiner
Posts: 11
Joined: Tue Feb 28, 2017 8:12 am

Re: Event delivery (CVode.event) not impacting simulation time

Post by ngreiner »

Dear Ted,
Thank you for your reply and sorry for the delay of mine.
I ran a test using played vector with interpolation to drive my extracellular batteries instead of abrupt changes delivered using the event delivery system. It successfully got rid of the IDA initialization failures.
However, I naively created as many vectors as there were segments in my model, which I assume may become problematic for models with a lot of segments.
Following your suggestion, I looked into the xtra mechanism to see how I could use one single vector and use transfer resistances for each segment, and I have one last question about it: to pass the computed extracellular potential from the xtra mechanism to the extracellular mechanism, your scripts used the idiom 'setpointer ex_xtra(x), e_extracellular(x)' which is a HOC statement. Is there an equivalent statement in Python ?
(I have tried 'h.setpointer(a(x)._ref_ex_xtra, a(x)._ref_e_extracellular)' and 'a(x)._ref_ex_xtra.setpointer(a(x)._ref_e_extracelluler)', without success.)
Nathan
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Event delivery (CVode.event) not impacting simulation time

Post by ted »

I naively created as many vectors as there were segments in my model, which I assume may become problematic for models with a lot of segments.
Can eat a lot of storage and slow things down.
to pass the computed extracellular potential from the xtra mechanism to the extracellular mechanism, your scripts used the idiom 'setpointer ex_xtra(x), e_extracellular(x)' which is a HOC statement. Is there an equivalent statement in Python ?
One might expect this to be documented in the Python version of the Programmers' Reference. It is, but simply going to the Index and clicking on setpointer brings up the page "HOC Keywords" which describes itself as "a page about HOC syntax; it is not directly applicable to Python-based simulations." Sure enough, everything on that page is pure hoc, including the documentation of setpointer.

Instead I had to resort to the "Quick search" widget in the left column of the Programmers' Reference. Searching for instances of setpointer turned up four hits. The first hit links to the "HOC Keywords" page, but the other three link to information about using setpointer via python. I found each of them to contain potentially useful information, but for your purpose the essential one is to material on the "Python Language" page which includes an example of using setpointer to couple a density mechanism's range variable pointer to a driving variable.
Post Reply