Event delivery (CVode.event) not impacting simulation time

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

Event delivery (CVode.event) not impacting simulation time

Post by ngreiner » Thu Oct 11, 2018 12:01 pm

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: 9
Joined: Tue Feb 28, 2017 8:12 am

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

Post by ngreiner » Fri Oct 12, 2018 9:28 am

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: 5372
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 » Fri Oct 12, 2018 12:00 pm

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.

Post Reply