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

[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]

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]

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]

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]

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.

[jseo977]

Hello Ted,

Apologies for the delayed response!
I plotted gvec for each segment over time and got the expected response (for my case, the extracellular potential at each point is scaled sinusoidally for one period).

However, I am not sure how to plot the segment's e_extracellular over time.
When I run `seg.e_extracellular` (after assigning `gvec.play(seg.e_extracellular, stim_t_vec)`, also tried `seg._ref_e_extracellular`), I get a single value. How can I get the time-varying e_extracellular for a segment?

Further, do you have any recommended resources to learn how to use xtra and extracellular mechanisms for Python? Maybe I am missing some very fundamental concepts about how Neuron('s extracellular mechanism) works.

Thanks!

[ted]

Assuming that the time course of seg's extracellular potential is specified by the vector pair gvec and stim_t_vec, then the correct syntax for Vector.play is

gvec.play(seg._ref_e_extracellular, stim_t_vec)

This should force seg.e_extracellular to follow the time course specified by the vector pair gvec and stim_t_vec.

To see what happens to the potential at the outer surface of seg, you need to see the time course of seg.vext[0]. You could use Vector.record to capture that to a vector, and then after a simulation plot the recorded values vs. the values in stim_t_vec.

You will find documentation of extracellular at nrn.readthedocs.io. The entry for python is at https://nrn.readthedocs.io/en/latest/ho ... ml#index-5

https://www.neuron.yale.edu/ftp/ted/neu ... sphere.zip presents examples of extracellular stimulation in hoc and Python. These examples exploit the assumptions that

1. the extracellular medium is linear and purely resistive
2. the normalized time courses of the current delivered by the extracellular electrodes are identical

These assumptions imply that
1. the normalized time course of extracellular potential is identical at all locations
and
2. extracellular potential at any point xyz is proportional to the stimulus current

Consequently, for monopolar stimulation, the extracellular potential at (x,y,z) at time t equals
Rx(x,y,z) i(t)
where i(t) is the time course of the stimulus current delivered by the stimulating electrode
and Rx(x,y,z) is the transfer resistance between that electrode and location (x,y,z).

The same formula holds for multipolar stimulation (2 or more stimulating electrodes), where i(t) is the time course of the stimulus current delivered by one electrode, and the value of Rx(x,y,z) depends on (x,y,z), the locations of all stimulus electrodes, the apportionment of currents over those electrodes (which remains constant), and the geometry and conductive properties of the medium.

How to discover the Rx values? For monopolar stimulation, set i to 1 mA and see what that does to extracellular potential (in mV) at all locations of interest. From those values, calculate the corresponding rx values in megohms as 1e-6 * (extracellular potential in mV)/(1 mA)

For multipolar stimulation, pick one of the electrodes to be the "special electrode". Then set the current it delivers to 1 mA, and set the currents delivered by the other electrodes to be whatever they should be according to how you designed your stimulus. Example: for bipolar stimulation, the "other" electrode should deliver a current of -1 mA; for tripolar stimulation you might want both "other" electrodes to deliver -0.5 mA (remember to preserve current balance). Finally, calculate the extracellular potentials in mV at all locations of interest in your model cell, and finally calculate the corresponding rx values (in megohms) for those locations as 1e-6 * (extracellular potential in mV)/(1 mA)