h.run vs h.continuerun(h.tstop)

The basics of how to develop, test, and use models.
Post Reply
jseo977
Posts: 3
Joined: Wed Sep 25, 2024 1:47 am

h.run vs h.continuerun(h.tstop)

Post by jseo977 »

Hello all,

Apologies if this is a very novice question. I am just getting started but struggling to understand the difference between h.run and h.continuerun(h.tstop), after explicitly defining h.tstop in both situations.

I found that h.continuerun(h.tstop) applies my h.finitialize(-60 * mV), but h.run() seems to ignore this and sets it to the default -65 mV.
However, h.run() applies the extracellular potential field (that I precomputed in COMSOL) to the axon, but h.continuerun(h.tstop) doesn't do this.
This makes it very annoying as I cannot get both constraints in the same simulation.

Further, when I use h.run(), I seem to get a remnant of the previous simulation at the beginning of my code. What could be the reason for this?
jseo977
Posts: 3
Joined: Wed Sep 25, 2024 1:47 am

Re: h.run vs h.continuerun(h.tstop)

Post by jseo977 »

For reference, my code looks like this:

Cell 1:

Code: Select all

# Create a simple nerve (axon) model
axon = h.Section(name='axon')
axon.L = 2000
axon.diam = 10 
axon.nseg = 101

h.load_file('stdrun.hoc')

axon.insert('hh')
axon.insert('extracellular')

#Set up recording variables
v_rec_0_33 = h.Vector().record(axon(0.33)._ref_v)
v_rec_0_5 = h.Vector().record(axon(0.5)._ref_v)
v_rec_0_66 = h.Vector().record(axon(0.66)._ref_v)
t_rec = h.Vector().record(h._ref_t)
Cell 2:

Code: Select all

# Initialize simulation
h.finitialize(-60 * mV)

# Simulation parameters (ms)
stim_start = 20
stim_duration = 1
h.tstop = 50

time_points = np.arange(stim_start, stim_start + stim_duration, h.dt)
time_points = np.insert(time_points, 0, stim_start - h.dt)
time_points = np.append(time_points, stim_start + stim_duration)

spatio_temporal_wave_values = np.tile(axons_interpolated_extracellular_voltages[3], (len(time_points), 1)).T
spatio_temporal_wave_values = np.insert(spatio_temporal_wave_values, 0, 0)
spatio_temporal_wave_values = np.append(spatio_temporal_wave_values, 0)
stim_t_vec = h.Vector(time_points)
for i, seg in enumerate(axon):
    gvec = h.Vector(spatio_temporal_wave_values[i])
    gvec.play(seg._ref_e_extracellular, stim_t_vec)

#h.continuerun(h.tstop)
h.run()

# Plot results
plt.figure(figsize=(13, 5))
plt.plot(list(t_rec), list(v_rec_0_33))
plt.plot(list(t_rec), list(v_rec_0_5))
plt.plot(list(t_rec), list(v_rec_0_66))

plt.legend(['0.33', '0.5', '0.66'])
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')

plt.axvline(x=stim_start, color='r', linestyle='--', label=f'Stimulus Start ({stim_start} ms)', alpha=0.2)
plt.axvline(x=stim_start + stim_duration, color='g', linestyle='--', label=f'Stimulus End ({stim_start + stim_duration} ms)', alpha=0.2)

plt.title('Membrane potential in response to extracellular field stimulation')
plt.show()
I have been running cell 1 just once, and then changing parameters in the beginning of cell 2. I found that when I run both cell 1 and cell 2 more than once, I get an error where t_rec is only of length 1.
ted
Site Admin
Posts: 6358
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: h.run vs h.continuerun(h.tstop)

Post by ted »

Decide whether you want to use h.run() or the sequence
h.finiitalize
h.continuerun

If you prefer the latter, it's probably best to avoid inserting any other statements between them.

If you prefer to call h.run(), there's no point in calling h.finiitalize, because h.run() executes hoc statements that include
finitialize(v_init)
and v_init has its own default value (-65) that may not be what you want. So if you prefer to call h.run(), first assign the desired value to h.v_init.

Regardless of your choice, make sure that all of these things have already been done:
execution of all model specification code
complete specification of instrumentation e.g. vector record, vector play, attachment of signal sources like IClamp, SEClamp, or synaptic mechanisms

Finally, keep your code properly organized. If you're going load stdrun.hoc, do that before doing any specification of biological properties e.g. creation of sections, parameter assignments.
jseo977
Posts: 3
Joined: Wed Sep 25, 2024 1:47 am

Re: h.run vs h.continuerun(h.tstop)

Post by jseo977 »

Thank you so much for your prompt response.

I have now updated my code to below (in a Jupyter notebook), and notice that when I run the cell more than once in the same kernel, I get an error that says

Code: Select all

plt.plot(list(t_rec), list(v_rec_0_33))
^
ValueError: x and y must have same first dimension, but have shapes (0,) and (2001,)
It seems like with my new code, the t_rec doesn't record or update the times after the first run. This must be something done under the hood that I am not fully understanding. Could you provide some guidance? Apologies for the inconvenience.

Another issue I find with this implementation is that the extracellular potentials don't get applied to the segment/axon, so I cannot see any perturbation to the membrane voltage. I see that due to the inclusion of the hh mechanism that the neuron's membrane voltage does still return to ~ -65 mV after setting h.v_init to -60 mV.

Code: Select all

h.load_file('stdrun.hoc')

# Create a simple nerve (axon) model
axon = h.Section(name='axon')
axon.L = 6000 * mm
axon.diam = 10 * um
axon.nseg = nseg

axon.insert('hh')
axon.insert('extracellular')

#Set up recording variables
v_rec_0_33 = h.Vector().record(axon(0.33)._ref_v)
v_rec_0_5 = h.Vector().record(axon(0.5)._ref_v)
v_rec_0_66 = h.Vector().record(axon(0.66)._ref_v)
t_rec = h.Vector().record(h._ref_t)

# Simulation parameters (ms)
stim_start = 20
stim_duration = 10
h.tstop = 50 * ms
h.v_init = -60 * mV

axon_type = 4

spatio_temporal_wave_values = time_varying_v_interpolation[:, axon_type, :]*1000
stim_t_vec = h.Vector(np.array(v_times) + stim_start)
for i, seg in enumerate(axon):
    gvec = h.Vector(spatio_temporal_wave_values.T[i])
    gvec.play(seg._ref_e_extracellular, stim_t_vec)

h.run()

# Plot results
plt.figure(figsize=(13, 5))
plt.plot(list(t_rec), list(v_rec_0_33))
plt.plot(list(t_rec), list(v_rec_0_5))
plt.plot(list(t_rec), list(v_rec_0_66))

plt.legend(['0.33', '0.5', '0.66'])
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential (mV)')

plt.axvline(x=stim_start, color='r', linestyle='--', label=f'Stimulus Start ({stim_start} ms)', alpha=0.2)
plt.axvline(x=stim_start + stim_duration, color='g', linestyle='--', label=f'Stimulus End ({stim_start + stim_duration} ms)', alpha=0.2)

plt.title('Membrane potential in response to extracellular field stimulation')
plt.show()
ted
Site Admin
Posts: 6358
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: h.run vs h.continuerun(h.tstop)

Post by ted »

when I run the cell more than once in the same kernel, I get an error
That's strange. What do you mean "run more than once in the same kernel"? just call h.run() again, and then execute the "plot results" code? That should work if you did it in plain python, without using Jupyter.
the extracellular potentials don't get applied to the segment/axon, so I cannot see any perturbation to the membrane voltage
This is a big problem, so it's important to be very sure of what's happening. What do the plots of each segment's e_extracellular vs. time look like? And do they look like the plot of gvec vs. stim_t_vec?
the neuron's membrane voltage does still return to ~ -65 mV after setting h.v_init to -60 mV.
Yes, the HH mechanism's resting potential is -65 mV, so no matter what you do to membrane potential at t==0, it will always eventually go back to -65 mV. There are ways to force a model with the HH mechanism to have a different resting potential, but that would affect the excitability of the model.
Post Reply