## Stimulating a Cell

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Stimulating a Cell

Can anyone help me understand why this code elicits no response in the cell? I've confirmed it's excitable with h.IClamp(), and this code is an adaptation of Ted's from another post. I don't understand the problem.

Appreciate it!

Code: Select all

``````# plays a waveform into xtra's is
from neuron import h, gui
from matplotlib import pyplot as plt
import numpy as np

dend = h.Section(name='dend') # name='dend' for nice section names in NEURON's graphs
dend.insert('pas')
dend.insert('hh')
dend.insert('extracellular') # this mechanism's defaults are OK
dend.insert('xtra')
dend.nseg = 5 # so we can see variation over distance
for seg in dend:
seg.xraxial[0] = 1e9
seg.xraxial[1] = 1e9
seg.xg[0] = 1e9
seg.xg[1] = 1e9
seg.xc[0] = 0
seg.xc[1] = 0

seg.xtra._ref_im = seg._ref_i_membrane
seg.xtra._ref_ex = seg._ref_e_extracellular
seg.xtra.x = 0
seg.xtra.y = 0
seg.xtra.z = 0

# Vector.play with interpolation needs 6 points to define a rectangular pulse
# successive elements of tsvec and isvec specify the x,y coordinates
# of the pulse waveform's "breakpoints"
time = np.arange(6)
t = np.insert(time, 3,2)
t = np.insert(t,5,4)
mag = [0,0,0,1e4,1e4,1e4,0,0]
tsvec = h.Vector().from_python(t)
isvec = h.Vector().from_python(mag)

# now drive h.is_xtra
isvec.play(h._ref_is_xtra, tsvec, 1)

ti = h.Vector().record(h._ref_t)
v = h.Vector().record(dend(0.5)._ref_v)

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

plt.plot(tsvec,isvec)
plt.show()

plt.plot(ti, v)
plt.show()``````
ted
Posts: 6127
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Stimulating a Cell

How to figure out what's happening:

Reduce amplitude of stimulus current to 1. If that turns out to be too small, you can always increase it later.
Walk the model through a simulation, one step at a time, and examine plots of e_extracellular and v (membrane potential) along the length of the stimulated section to see how these variables evolve over time and space.

Unfortunately matplotlib which blocks simulation execution when its plots are displayed, so you'll have to use NEURON's range variable plots. A range variable plot shows the values of a range variable (like e_extracellular and v) along a user-specified path, and it can be linked to NEURON's run time system so that it updates automatically with each fadvance()). See

Code: Select all

``https://nrn.readthedocs.io/en/8.2.2/python/visualization/rvarplt.html#rangevarplot``
.

Also follow the hints at the end of this post.

To manually walk through a simulation, bring up a RunControl panel, click on Init, then click on
"Single Step" to advance by dt
"Continue til" to advance to a specified stop time
"Continue for" to advance by a specified number of ms

What you observe will give you something new and substantive (pertaining to extracellular stimulation) to think and ask questions about.

Hints:

1. Get rid of dend.insert('pas'). hh includes its own passive (leak) channel. Adding pas gives you something that deviates from the original HH model and has a resting potential that is not the same as the original HH model's resting potential.

2. Comment out all statements from h.finitialize to plt.show

3. Append your code for creating the range variable plots, then save your file.

3. At the system command prompt, execute

python3 test.py # assuming python3 is your python executable, and test.py is the name of the file that contains your python code

Then use the NEURON Main Menu toolbar to bring up a RunControl panel. Use its Init and other buttons as described above. If necessary, rescale the y axes of the range variable plots (easily done via the GUI or by using the Graph class's size() method as shown in the Programmer's Reference examples.
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

Thanks for the advice, Ted. Full disclosure: my familiarity with the GUI and the native plotting functions is limited. I made this effort, but the output just flat-lines. So, I'm not sure - is the problem in how I set up the plot, ran the simulation, or implemented the stimulus?

Code: Select all

``````# plays a waveform into xtra's is
from neuron import h, gui
from matplotlib import pyplot as plt
import numpy as np

dend = h.Section(name='dend') # name='dend' for nice section names in NEURON's graphs
dend.insert('hh')
dend.insert('extracellular') # this mechanism's defaults are OK
dend.insert('xtra')
dend.nseg = 5 # so we can see variation over distance
for seg in dend:
seg.xraxial[0] = 1e9
seg.xraxial[1] = 1e9
seg.xg[0] = 1e9
seg.xg[1] = 1e9
seg.xc[0] = 0
seg.xc[1] = 0
seg.e_extracellular = -65

#  h.setpointer(secA(0.5)._ref_v, 'V_pre', synAB) is 100% equivalent to: synAB._ref_V_pre = secA(0.5)._ref_v
# h.setpointer(dend(seg.x)._ref_e_extracellular, 'ex', dend(seg.x).xtra)
seg.xtra._ref_im = seg._ref_i_membrane
seg.xtra._ref_ex = seg._ref_e_extracellular
# seg.xtra.rx = 1
seg.xtra.x = 0
seg.xtra.y = 0
seg.xtra.z = 0

# Vector.play with interpolation needs 6 points to define a rectangular pulse
# successive elements of tsvec and isvec specify the x,y coordinates
# of the pulse waveform's "breakpoints"
time = np.arange(6)
t = np.insert(time, 3,2)
t = np.insert(t,5,4)
mag = [0,0,0,1e0,1e0,1e0,0,0]
tsvec = h.Vector().from_python(t)
isvec = h.Vector().from_python(mag)

# now drive h.is_xtra
isvec.play(h._ref_is_xtra, tsvec, 1)

ti = h.Vector().record(h._ref_t)
v = h.Vector().record(dend(0.5)._ref_v)

rvp = h.RangeVarPlot('v', dend(0), dend(1))
g = h.Graph()
g.size(0, 1002, 0, 2)
g.flush()

rvp.plot(g, 2, 3)

v_init = -65

h.finitialize(v_init)
h.fcurrent()

xvec = h.Vector()
yvec = h.Vector()
rvp.to_vector(yvec,xvec)

xvec.printf()
yvec.printf()``````
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

I think I'm not implementing the RangeVarPlot correctly; it also shows no activity using the same h.IClamp() instance which successfully excites the cell in a matplotlib graph. Is there additional documentation on RangeVarPlot I can study, if indeed this is the issue?
ted
Posts: 6127
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Stimulating a Cell

A lot of your code is good, but you're having to wrestle with many different, unfamiliar pieces. Also, not all documentation is as clear or complete as it should be. And, unfortunately debugging is complicated by the fact that, while matplotlib + NEURON is a powerful combination, it doesn't enable interactive simulations. So here's a working example with some discussion of what I left out and what I included.

First the code that defines the model cell.

Code: Select all

``````from neuron import h, gui
# not using numpy or matplotlib

##### model cell
dend = h.Section(name='dend')
# not specifying dend.L or diam; they will have their default values 100 and 500 um
dend.insert('hh')
``````
Set the stage for extracellular stimulation

Code: Select all

``````##### instrumentation
## stimulus first
dend.insert('extracellular') # this mechanism's default params are OK
dend.insert('xtra')
dend.nseg = 5 # so we can see variation over distance
for seg in dend:
seg.xtra._ref_im = seg._ref_i_membrane
seg.xtra._ref_ex = seg._ref_e_extracellular
## Vector play with interpolation
t =   [0, 1, 1, 2, 2, 3]
mag = [0, 0, 1, 1, 0, 0]
tsvec = h.Vector().from_python(t)
isvec = h.Vector().from_python(mag)
isvec.play(h._ref_is_xtra, tsvec, 1)
``````
Next the graphs, starting with membrane potential. Note how updating a graph at each advance is enabled by appending the graph to graphList[0].

Code: Select all

``````## dend(0.5).v vs. t
gv = h.Graph()
gv.addvar('dend(0.5).v', dend(0.5)._ref_v, 3, 2) # blue (3), 2 pixels wide
h.graphList[0].append(gv) # to update vg at every fadvance()
gv.size(0, h.tstop, -80, 40)
## v along dend ("space plot")
vrvp = h.RangeVarPlot('v', dend(0), dend(1))
gvrvp = h.Graph()
gvrvp.addobject(vrvp, 3, 2) # 3, 2 == blue, 2 pixels wide
# adjust x and y axes appropriately
gvrvp.size(0, 100, -80, 40)
h.graphList[0].append(gvrvp) # update graph at every fadvance()``````
Also plot e_extracellular and vext[0] (to confirm that these variables are being driven).

Code: Select all

``````## dend.vext[0](0.5) and dend.e_extracellular(0.5) vs. t
gvext = h.Graph()
gvext.addvar('dend(0.5).e_extracellular', dend(0.5)._ref_e_extracellular, 2, 4) # 2, 4 == red, 4 pixels wide
gvext.addvar('dend(0.5).vext[0]', dend(0.5)._ref_vext[0], 3, 2) # 3, 2 == blue, 2 pixels wide
h.graphList[0].append(gvext) # to update graph at every fadvance()
gvext.size(0, h.tstop, 0, 100)
## e_extracellular and vext[0] along dendrite
eextrvp = h.RangeVarPlot('e_extracellular', dend(0), dend(1))
vextrvp = h.RangeVarPlot('dend.vext[0]', dend(0), dend(1))
gextrvp = h.Graph()
gextrvp.addobject(eextrvp, 2, 3) # red, 3 pixels wide
gextrvp.addobject(vextrvp, 3, 2) # blue, 2 pixels wide
# adjust x and y axes appropriately
gextrvp.size(0, 100, 0, 100)
h.graphList[0].append(gextrvp) # to update graph at every fadvance()``````
Finally, simulation flow control: a RunControl panel for launching and stepping through simulations

Code: Select all

``````##### simulation flow control

Three different ways you can run a simulation.

1. Click on the RunControl panel's Init & Run button. Easy, but too fast to see the graphs evolve.
2. Initialize, the manually step through a simulation. Click on RunControl/Init. Then click on Single Step to advance by dt, or on "Continue for" to advance for the number of ms specified in the adjacent numerical field.
3. Bring up a Movie Run tool by clicking on NEURON Main Menu / Tools / Movie Run. Then click on the Movie Run tool's own Init & Run button.

Now some things for you to try and think about.

1. You're not getting a spike. In fact, membrane potential isn't changing even though e_extracellular is rising to 1e6 mV (1000 V!). The problem is that the surface of the model cell remains isopotential throughout the simulation. Result: no matter how strong the stimulus is, it can't induce transmembrane current flow. Look up Gauss's law (which one? there are so many . . .).

The cell surface isopotential because the transfer resistance between the stimulus source and the surface of the cell is uniform over the surface of the cell. A quick way to fix this is to enter the following two lines at Python's >>> prompt :

Code: Select all

``````for seg in dend:
seg.xtra.rx = seg.x``````
The result emulates an extracellular field oriented parallel to dend's long axis with a steep linear gradient. Now run a simulation (with Movie Run's Init & Run) and see what happens.

Nothing that looks like a spike--because the potential gradient is way too intense: 1000 volts per 100 um length. That's a few orders of magnitude larger than what it takes to ionize dry air.

What to do? Try this one liner at the >>> prompt:
for seg in dend: seg.xtra.rx /= 10
(you'll have to hit Return twice to make that execute). Then run another simulation. In fact, follow this algorithm yourself:

Code: Select all

``````REPEAT
for seg in dend: seg.xtra.rx /= 10
(hit Return key twice)
run a simulation
UNTIL the stimulus is small enough to elicit a reasonable looking spike.``````
Now what is the value of dend(0.5).xtra.rx ?
ted
Posts: 6127
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Stimulating a Cell

One more comment: the purpose of some of these graphs is to verify that the extracellular stimulus is actually happening, and reaching the external surface of the model cell--and that it makes sense.
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

This issue has kind of sprawled for me as I looked into your suggestions - sorry if the breadth of this post is inappropriate for the forum setting, but hopefully others who are facing the same issues will find this useful.

So one thing I realized, as I was going through the read-me associated with the 'xtra' mechanism, is that I needed (okay, wanted) to implement a version of interp_zyx.hoc and calc_rxc.hoc functions in Python. This, I did as follows.

Code: Select all

``````# use interp_xyz at the SECTION level
def interp_xyz(section):
nn = xr = 0
sec = section #h.cas()

# get the data for each section
nn = int(sec.n3d()) # nn is the number of pt3d data for each section
xx = h.Vector(nn)
yy = h.Vector(nn)
zz = h.Vector(nn)
length = h.Vector(nn)
for ii in range(nn):
xx.x[ii] = sec.x3d(ii)
yy.x[ii] = sec.y3d(ii)
zz.x[ii] = sec.z3d(ii)
length.x[ii] = sec.arc3d(ii)

# h.Vector().interpolate() requires normalized section lengths
# so, length contains irregularly spaced distances of the pt3d data along the section's centroid
length.div(length.x[nn-1])

# creating the destination vector
# so, range_var contains regularly spaced normalized distances of the nodes along the section's centroid
range_var = h.Vector(definition+2)
range_var.indgen(1/(definition+1))
# range_var.sub(1/(2*definition)) # this will never creation regularly spaced distances
# range_var.x[0] = 0 # and this is redundant, indgen does it automatically
range_var[definition +1] = 1

xint = h.Vector(definition+2)
yint = h.Vector(definition+2)
zint = h.Vector(definition+2)

xint.interpolate(range_var, length, xx)
yint.interpolate(range_var, length, yy)
zint.interpolate(range_var, length, zz)

# for each node, assign the coordinate values to xtra.x, xtra.y, and xtra.z
# note that the loop excludes the first last and points
for ii in range(1, definition+1):
xr = range_var.x[ii]
section(xr).xtra.x = xint.x[ii]
section(xr).xtra.y = yint.x[ii]
section(xr).xtra.z = zint.x[ii]``````

Code: Select all

``````# use calc_rxc at the SEGMENT level
def calc_rxc(seg, e1Vec, e2Vec=[1e9,1e9,1e9], stim_type='Uni', rho=575, view=False):
# rho defaults to 575 ohms.cm, the second/ground electrode defaults to 'infinitely' far away, and the type of stimulation defaults to unipolar
# rho represents the resistivity of the bath solution, assumed to be linear and constant
# print("the segment is ", seg)

configs = ['uni', 'bi']
# this failsafe isn't working yet - come back to it later
# while (stim_type not in configs):
#     stim_type = input("Invalid electrode configuration. (Valid types are: uni, bi.) Please try again:\n")
#     print(stim_type)

# seg_idx = seg.x # used in original code to delineate ineligible segment nodes; apparently unnecessary in python?

x1 = e1Vec[0]
y1 = e1Vec[1]
z1 = e1Vec[2]

x2 = e2Vec[0]
y2 = e2Vec[1]
z2 = e2Vec[2]

def set_rx_uni(seg, X1, Y1, Z1):
r = h.sqrt( (seg.xtra.x - X1)**2 + (seg.xtra.y - Y1)**2 + (seg.xtra.z - Z1)**2 )
if r == 0:
r = seg.diam/2
seg.xtra.rx = (rho / 4 /h.PI)*(1/r)*0.01

def set_rx_bi(seg, X1,Y1,Z1,X2,Y2,Z2):
for sec in h.allsec():
for seg in sec:
r1 = h.sqrt( (seg.xtra.x - X1)**2 + (seg.xtra.y - Y1)**2 + (seg.xtra.z - Z1)**2 )
r2 = h.sqrt( (seg.xtra.x - X2)**2 + (seg.xtra.y - Y2)**2 + (seg.xtra.z - Z2)**2 )
if r1 == 0:
r1 = seg.diam/2
if r2 == 0:
r2 = seg.diam/2
seg.xtra.rx = h.abs((rho / 4 /h.PI)*(1/r1 - 1/r2)*0.01)

def put_sElec(X1, Y1, Z1):
for sec in h.allsec(sElec):
sec.pt3dclear()

def setelec(seg, xe, ye, ze):
set_rx_uni(seg, xe, ye, ze)
# put_sElec(xe, ye, ze)

# not working at the moment -- irrelevant to my current priorities
if view == True:
sElec = h.Section(name='sElec')
put_sElec(x1,y1,z1)
gElec = h.Shape(0)
gElec.view(-245.413, -250, 520.827, 520, 629, 104, 201.6, 201.28)
pElec = h.PointProcessMark(0.5)
sElec = h.PointProcessMark(0.5)
gElec.point_mark(pElec, 2)

match stim_type:
case 'Uni':
set_rx_uni(seg, x1,y1,z1)
case 'Bi':
set_rx_bi()``````
At this point, I hit my first question. Also included in the 'xtra' package is field.hoc -- what this function does exactly, I can't tell. Is this necessary to correctly stimulate a cell? Or is it ancillary/involved in other aspects of the simulation?

With these inserted into the cell creation (as below), the cell reacts *seemingly* appropriately. However, I have another question. In my first post, I used the 'hh' channel to make the cell excitable. I'm also working a variation of the Hodgkin-Huxley equations implemented in an NMODL file. In this file, I specify a membrane capacitance separate from the h.Section().cm capacitance. Is this redundant? Worse, are the two capacitances interfering with each other somehow? I noticed that the cell is taking an inordinately long time to equilibrate, ~100ms for membrane voltage to stabilize.

Code: Select all

``````class CellMake():
def __init__(self, filename, model):

self.cell = Creator(1, filename, definition) # personal model generator, basically equivalent to h.Import3d_SWC_read()
h.define_shape()

for sec in h.allsec():
# sec.insert('hh')
sec.insert('li_channels')
sec.insert('extracellular')
sec.insert('xtra')

sec.Ra = 150    # Axial resistance in Ohm * cm
sec.cm = capacitance     # Membrane capacitance in micro Farads / cm^2
interp_xyz(sec)

for seg in sec:
seg.xraxial[0] = 1e9
seg.xraxial[1] = 1e9
seg.xg[0] = 1e9
seg.xg[1] = 1e9
seg.xc[0] = 0
seg.xc[1] = 0

seg.xtra._ref_im = seg._ref_i_membrane
seg.xtra._ref_ex = seg._ref_e_extracellular
calc_rxc(seg, posE_vec)

seg.li_channels.celsius = 37
seg.li_channels.Cm = capacitance
``````
Lastly, as I said, the cell's reaction *seems* reasonable, but I'm no expert. If, in truth, it isn't a reasonable reaction to a stimulus of this magnitude (3uA/cm2, I eventually found), do you have a thought as to where else the problem could lie?
ted
Posts: 6127
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Stimulating a Cell

I needed (okay, wanted) to implement a version of interp_zyx.hoc and calc_rxc.hoc functions in Python
Have you verified that your Python code properly calculates the locations of segment centers? One way to demonstrate this is graphically. For example, treat the x3d data from a real neurite as samples from a piecewise linear function fx whose independent variable is "range" i.e. the normalized anatomical path length along the original neurite from its proximal end (end closest to the soma). The values of the dependent variable are the actual x3d values. Since fx is piecewise linear, a plot of fx vs. range will be a line that zigs and zags up and down as range increases from 0 to 1. For any value of nseg you know there will be n segments, whose indices are 0..n-1, and the center of the ith segment will correspond to a range value of (0.5 + i)/nseg. The x coordinate of the point along the actual neurite that corresponds to the ith segment's center will be the value of fx((0.5 + i)/nseg). The same analysis can be applied to the y3d and z3d data. For any value of nseg, the interpolated segment centers should lie on the plots of fx, fy, and fz vs. range.

Have you verified that your calc_rxc() generates correct extracellular potentials to emulate a uniform extracellular gradient or a gradient that would be produced by a monopolar electrode in an infinite, uniform nondispersive medium?

field.hoc calculates the potential observed at a position in space as a consequence of a model neuron's transmembrane current flow. It is useful if you are interested in extracellular recording.
In this file, I specify a membrane capacitance separate from the h.Section().cm capacitance. Is this redundant?
Yes. If you have a section, you already have a cable with longitudinal resistance and membrane capacitance. Does eliminating this redundancy take care of the problems you encountered?
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

I'm going to take these one at a time, starting with interp_xyz().. I believe the function I posted above is incorrect. I changed the range_var block when I realized I was selecting the edge of the segment, rather than the center.

Code: Select all

``````def interp_xyz(section):
nn = xr = 0
sec = section #h.cas()

# get the data for each section
nn = int(sec.n3d()) # nn is the number of pt3d data for each section
xx = h.Vector(nn)
yy = h.Vector(nn)
zz = h.Vector(nn)
length = h.Vector(nn)
for ii in range(nn):
xx.x[ii] = sec.x3d(ii)
yy.x[ii] = sec.y3d(ii)
zz.x[ii] = sec.z3d(ii)
length.x[ii] = sec.arc3d(ii)

# for i in range(nn+1):
#     print("length ", length.x[i])
# h.Vector().interpolate() requires normalized section lengths
# so, length contains irregularly spaced distances of the pt3d data along the section's centroid
length.div(length.x[nn-1])

# THIS IS THE ALTERED BLOCK - I PREVIOUSLY FAILED TO RECOGNIZE WHAT THE BLOCK IS DOING SPATIALLY
range_var = h.Vector(definition+2)
range_var.indgen(1/definition)
range_var.x[definition] = 1
range_var.sub(1/(2*definition))
range_var.x[0] = 0
range_var.x[definition+1] = 1

# this loop was used for testing, to confirm correct range_var values
# for i in range(definition+2):
#    print(range_var.x[i])

xint = h.Vector(definition+2)
yint = h.Vector(definition+2)
zint = h.Vector(definition+2)

xint.interpolate(range_var, length, xx)
yint.interpolate(range_var, length, yy)
zint.interpolate(range_var, length, zz)

# for each node, assign the coordinate values to xtra.x, xtra.y, and xtra.z
# note that the loop excludes the first last and points
for ii in range(1, definition+1):
xr = range_var.x[ii]
if xr== 0.5:
print(section(xr),xint.x[ii]) # this print statement spits out the interpolated segment center values
section(xr).xtra.x = xint.x[ii]
section(xr).xtra.y = yint.x[ii]
section(xr).xtra.z = zint.x[ii]``````
I tested its accuracy with this little script.

Code: Select all

``````sectionA = cell_n[0].axon[22] # a random neurite
count = sectionA.n3d()
for i in range(count-1):
d.append(sectionA.x3d(i))
r.append(i/(count))

for i, j in zip(r, d):
if i == 0.5:
print(i, j) # print the neurite segment and the x3d value``````
Obviously, this isn't the graphical method you recommended, but I think it accomplishes the same thing? Yet the interp_xyz() output is unchanged, still inaccurate. Is my test script wrong? Or are there more problems in interp_xyz()?
Last edited by ajpc6d on Tue Feb 14, 2023 12:25 pm, edited 1 time in total.
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

Regarding calc_rxc() -- my understanding is that this function only calculates the resistance of the extracellular medium between the electrode and the neurite segment, and that this value is based on the spatial data from interp_xyz() and is passed on to xtra for calculation of the extracellular potential. (Correct?)
I can confirm that calc_rxc() accesses the segment center xyz points derived by interp_xyz() and that the resultant segment.xtra.rx values vary appropriately with the parameters given. Assuming interp_xyz() has done its job correctly, I think it should follow that calc_rxc() is also correct; but if this reasoning isn't valid, I know no other way to verify.

And regarding redundant membrane capacitance declaration -- eliminating the redundancy did not affect the long equilibration time (~100ms).
ajpc6d
Posts: 11
Joined: Mon Oct 24, 2022 6:33 pm

### Re: Stimulating a Cell

And one last bit -- I implemented fieldrec() in this way.

Code: Select all

``````def fieldrec(v):
summa = 0.0
for sec in h.allsec():
for seg in sec:
summa += float(seg.xtra.er)
v.append(summa)

#
# later in the code
#

vrec = []
pyFcall = (fieldrec,vrec)
h.CVode().extra_scatter_gather(0,pyFcall)
# after which h.finitialize(), h.fadvance(), etc...
``````
This is slightly changed from the .hoc implementation that I adapted in that it writes to a list/vector. I think this is a necessary change in Python? Certainly comparing the plot of the fieldrec() output with, for example, membrane voltage at the soma shows good agreement