Saving fixed timestep data using cvode variable timestep integration

[Robert Claar]

Hello,

I recently have had to do some filtering and additional processing of my model output that I didn't need to do before and now need to save my recorded data using a fixed timestep. However, I have a dramatic improvement in execution time using cvode compared to using a fixed timestep and I would like to preserve this efficiency if possible.

Reading the NEURON Vector().record() documentation it seems like my two options are providing a Dt value to determine the regular interval at which our data will be recorded

Code: Select all

vdest = vdest.record(var_reference, Dt)

or provide a time Vector, tvec, that will determine the sampling rate.

Code: Select all

vdest = vdest.record(var_reference, tvec)

Providing Dt seemed more desirable since we can just create the time vector using our desired Dt and run time.

In a Jupyter Notebook I used the magic command %%timeit to check the execution time of a simple simulation using a fixed timestep and using cvode along with providing the Dt argument.

Code: Select all

%%timeit
sec1 = h.Section(name="sec1")
sec1.insert(h.hh)
stim = h.IClamp(sec1(0))
stim.delay = 5
stim.dur = 1
stim.amp = 0.5
v_fixed = h.Vector().record(sec1(0.5)._ref_v, dt)
t_fixed = t = h.Vector().record(h._ref_t, dt)
h.cvode.active(True)
h.finitialize()
h.continuerun(run_time)

8.74 ms ± 18.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Code: Select all

%%timeit
sec2 = h.Section(name="sec2")
sec2.insert(h.hh)
stim = h.IClamp(sec2(0))
stim.delay = 5
stim.dur = 1
stim.amp = 0.5
h.cvode.active(False)
h.dt = 0.1
v_fixed = h.Vector().record(sec2(0.5)._ref_v)
t_fixed = t = h.Vector().record(h._ref_t)
h.finitialize()
h.continuerun(run_time)

11 ms ± 24.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

As a sanity check I also checked to see if providing a time Vector provided a benefit versus providing a Dt and it did not which was expected. So it seems my performance boost is slightly lost which, unless I am ignorant of something, makes sense because I'm assuming that my system of equations need to be solved at the provided time values based on Dt which essentially means I am using a fixed time step.

I also looked into h.cvode.record versus Vector().record() in case that would provide a performance boost but with the following code I got the following error and couldn't discern what to fix from the error message (likely syntax I'm guessing).

Code: Select all

%%timeit
sec4 = h.Section(name="sec4")
sec4.insert(h.hh)
stim = h.IClamp(sec4(0))
stim.delay = 5
stim.dur = 1
stim.amp = 0.5
h.cvode.active(True)
dt = 0.1
fixed_time_vec = np.arange(0,500 + dt, dt)
fixed_time_vec = h.Vector(fixed_time_vec)
v_fixed = h.Vector()
h.cvode.record(sec4(0.5)._ref_v, v_fixed, fixed_time_vec)
h.finitialize()
h.continuerun(run_time)

NEURON: Cvode.record  pointer not associated with currently accessed section
Use section ... (&var(x)...) intead of ...(&section.var(x)...)

 near line 0
 ^
        CVode[0].record(..., ..., ...)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[112], line 11
      9 fixed_time_vec = h.Vector(fixed_time_vec)
     10 v_fixed = h.Vector()
---> 11 h.cvode.record(sec4(0.5)._ref_v, v_fixed, fixed_time_vec)
     12 h.finitialize()
     13 h.continuerun(run_time)

RuntimeError: hocobj_call error

Another Idea I had is to keep using cvode and allow my data to not be recorded at regular intervals and then after the simulation resample my data using something like linear interpolation to get the data to have a regular sampling frequency.

If anyone could let me know if I am going about this in the best way that would be greatly appreciated.

[ted]

Thank you for an interesting post. This is what I'd do:

keep using cvode and allow my data to not be recorded at regular intervals and then after the simulation resample my data

This is best because it works, produces "correct" results, and saves your time.

NEURON's Vector class has an interpolate() method that will resample using linear interpolation at compiled code speed and can be called from hoc or Python. I am sure that Python has its own native method(s) for doing this.

For example, suppose a simulation has been run with adaptive integration, in which Vector.record was used to capture the time values and corresponding values of membrane potentiat at some point of interest to a pair of vectors called tvec and vvec, respectively. Also assume that the simulation ran from t==0 to some value tstop in ms that is evenly divisible by 0.1. Then

1. create a new pair of vectors called trs and vrs
2. initialize trs to hold the values of time at 0.1 ms intervals from 0 to tstop (that is, if tstop is 10 ms, trs will have 101 elements, and their values will be 0, 0.1, . . . 10)
3. finally, execute

vrs.interpolate(trs, tvec, vvec)

and vrs will contain the membrane potential values sampled at 0.1 ms, and the elements of trs will be the corresponding times. The hoc statements to do this are

objref trs, vrs
trs = new Vector()
trs.indgen(0, tstop, 0.1) // fill trs with values 0, 0.1, 0.2 . . . tstop
vrs.interpolate(trs, tvec, vvec)

Pythonistas are invited to add their own solutions. Extra points for code that is compact and easily understood.

[Robert Claar]

Thank you Ted for your quick reply!

In case anyone wants to just grab some code to do this:

NEURON Vectors

Code: Select all

t_resamp = h.Vector().indgen(start_time, end_time, Dt)
v_resamp = h.Vector().interpolate(t_resamp, t_orig, v_orig)

Python NumPy Arrays

Code: Select all

# You have to add Dt to end_time to include ent_time in your array
t_resamp = np.arange(start_time, end_time + Dt, Dt)
v_resamp = np.interp(t_resamp, t_orig, v_orig)

I did a quick search to try and find a native python method for linear interpolation to try and add it here but just found approaches for NumPy and SciPy (method scipy.interpolate.interp1d)