h.continuerun() behavior different when cvode is active

[Robert Claar]

Hello,

I have spent the last hour or so trying to find an old forum post relevant to my question, but can't seem to locate it.

I need to run my model for a certain amount of time for it to reach stable outputs. I also am generating a lot of data from my simulations so I am "chunking" my simulation, meaning, for example, if I run a simulation for an hour I will break it up into 5 minute chunks and at the end of each chunk will write the data and then resize my recording vectors using h.frecord_init(). Since I have no interest in any of the data before the model behavior stabilizes, I have set up my code in the following way:

Code: Select all

h.finitialize()
...
h.t = -stabilize_time

num_time_chunks = (stabilize_time + simulation_time) / chunk_time

for i in range(num_time_chunks):
	h.continuerun(h.t + chunk_time)
	
	if h.t >= 0:
		#write data
	# resize recording vectors
	h.frecord_init()

This obviously isn't my actual code, I just wanted it give you an idea of how it's structured. When I use a fixed time step this works perfectly. However, when I use cvode it breaks. For example I was running a simple test where

Code: Select all

stabilize_time = 1000
simulation_time = 1000
chunk_time = 250

and therefore the first chunk should run from -1000 to -750, but after the first call to continuerun() h.t becomes just slightly bigger than 0 and the data starts to be written to disk. I have fixed this by not using negative values for time so it's not causing me any more trouble per se, but I am just curious why this would not work for negative h.t values.

[ramcdougal]

Anytime you change something in a CVode simulation after initialization, you must run: h.CVode().re_init() to let CVode know. This is true for both time and state variables.

Code: Select all

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

soma = h.Section("soma")

t = h.Vector().record(h._ref_t)
h.CVode().active(True)

h.finitialize(-65)
h.t = -10
h.CVode().re_init()
h.continuerun(5)

print(list(t))

That said, you can avoid needing to do this by using an h.FInitializeHandler and initializing to the right time from the beginning.

Code: Select all

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

soma = h.Section("soma")

t = h.Vector().record(h._ref_t)
h.CVode().active(True)

def my_init():
    h.t = -10

fih = h.FInitializeHandler(my_init)

h.finitialize(-65)
h.continuerun(5)

print(list(t))

Since our cell is entirely unchanging, CVode will only take two steps in the interval (these are what will appear in the vector), and h.t will be advanced way beyond the interval. If you want variable step integration to end at a specific time, you need to declare an event there.

Code: Select all

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

soma = h.Section("soma")

t = h.Vector().record(h._ref_t)
h.CVode().active(True)

def my_init():
    h.t = -10

fih = h.FInitializeHandler(my_init)

TSTOP = 5

h.finitialize(-65)
h.CVode().event(TSTOP)
h.continuerun(TSTOP)
assert(h.t == TSTOP)

[Robert Claar]

That makes sense, thank you!