Getting the 'xtra' mechanism to work without 'extracellular'

Anything that doesn't fit elsewhere.
Post Reply
pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Getting the 'xtra' mechanism to work without 'extracellular'

Post by pascal » Wed Apr 26, 2017 10:11 pm

I am trying to get a minimal model working in which I record the LFP of one cell using 'xtra', without the 'extracellular' mechanism, as described here viewtopic.php?f=2&t=3389&p=14342&hilit= ... lel#p14342. I have made the modifications to xtra.mod described in the post, and I also have beforestep_py.mod working, which I implement at the beginning of my minimal model:

Code: Select all

from neuron import h, gui

# set up recording of LFP at every time step, in conjunction with Hines's beforestep_py.mod
s = h.Section()
bscallback = h.beforestep_callback(s(.5))

lfp_rec=[]
def callback():
    vx = 0
    for sec in h.allsec():
        for seg in sec:
            vx = vx + seg.er_xtra
    #end both for loops
    lfp_rec.append(vx)
    print h.t, vx

bscallback.set_callback(callback)

h('cvode.use_fast_imem(1)') #see use_fast_imem() at https://www.neuron.yale.edu/neuron/static/new_doc/simctrl/cvode.html

soma=h.Section(name='soma') 
Bdend=h.Section(name='Bdend')

soma.diam = soma.L = 20 #microns
soma.insert('hh')
Bdend.diam = 2 #microns
Bdend.L=200 #microns
Bdend.insert('pas')

Bdend.connect(soma,0,0) 

for sec in h.allsec():
    sec.insert('xtra')

h.define_shape() #this assigns default position of cell, with soma at origin
    
#call grindaway() in order to calculate centers of all segements, and--most importantly--copy these values to the xtra mechanism of each cell, so they can be used to calculate the LFP
h('xopen("interpxyz.hoc")') #from Ted's extracellular_stim_and_rec code; see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=8&t=3649
h('grindaway()') 
#set pointers
for sec in h.allsec():
    for seg in sec:
        h.setpointer(seg._ref_i_membrane_, 'im', seg.xtra)
        
#compute transfer resistances between electrode and all segements
h('xopen("calcrxc_a.hoc")') #this is a slightly altered version from Ted's original calcrxc.hoc
h('setelec(50,0,0)') #this is what actually sets electrode location and computes transfer resistances'''
    
#record soma voltage and time
t_vec = h.Vector()
v_vec_soma = h.Vector()
v_vec_dend = h.Vector()
v_vec_soma.record(soma(0.5)._ref_v) 
v_vec_dend.record(Bdend(0.5)._ref_v)
t_vec.record(h._ref_t)

#run simulation
h.tstop = 200 #ms
h.run()
And here is calcrxc_a.hoc (modified version of Ted's calcrxc.hoc, stripping out all the GUI-related stuff):

Code: Select all

rho = 351 
XE = -50  // um
YE = 0
ZE = 0

proc setrx() {  // now expects xyc coords as arguments
  forall {
    if (ismembrane("xtra")) {
      for (x,0) {
        r = sqrt((x_xtra(x) - $1)^2 + (y_xtra(x) - $2)^2 + (z_xtra(x) - $3)^2)
        rx_xtra(x) = (rho / 4 / PI)*(1/r)*0.01
      }
    }
  }
}

proc setelec() {
	xe = $1
	ye = $2
	ze = $3
	setrx(xe, ye, ze)
}
My problem is that attempting to run this simulation gives the error:
NEURON: division by zero, near line 1, setelec(50,0,0)

I checked whether the electrode was placed directly on top of one of the segment centers, and that is not the problem. I am actually able to get the program to run and give a coherent result if I insert the 'extracellular' mechanism into both sections, and also use Ted's setpointers() function (rather than using my own Python-based approach for setting pointers). But I really want to be able to get this to work without using 'extracellular', both for the performance speedup and so that I feel like I actually know what is going on. Thanks for the help!

ted
Site Admin
Posts: 5742
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by ted » Thu Apr 27, 2017 9:03 am

Your first task is to discover the point at which the error occurs. Either trace execution step by step, or insert print statements at multiple points in the code to indicate execution progress.

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Thu Apr 27, 2017 3:33 pm

Hmmm...I see what you're getting at now. Upon further investigation, the simulation was not actually choking within setelec, but rather somewhere within h.run(). The error message associated with this is simply "Segmentation fault (core dumped)". So yes, it would seem I need to somehow debug the actual execution of the simulation, as you say.

However, I've never before had to figure out exactly where in h.run() a simulation chokes. How do I "trace execution step by step"? Or, alternatively, how do I print statements within h.run()? Do I need to add some print statements to some file within the NEURON source code?

ted
Site Admin
Posts: 5742
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by ted » Fri Apr 28, 2017 12:23 am

Who said it happens in h.run? What makes you sure it's not happening during model setup, or that model setup isn't laying the ground for it to happen?

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Fri Apr 28, 2017 8:58 am

Model setup is laying the groundwork for this to happen, because print statements show the program gets to h.run and then chokes with Segmentation fault (core dumped). I'm having a hard time debugging from here, because usually when NEURON chokes within h.run, it gives an error message that indicates how I screwed up the model setup. Looking up the documentation for h.run, I see that it's just a wrapper for h.stdinit and h.continuerun, and that h.stdinit is a wrapper that includes h.init, which is a wrapper that includes h.finitialize. Using print statements I've been able to deduce that the segmentation fault occurs in h.finitialize.

Another big clue is that I can get the simulation to run without a segmentation fault by making two changes:
1. Insert 'extracellular' in both sections.
2. Use your setpointers() function (from extracellular_stim_and_rec) instead of

Code: Select all

for sec in h.allsec():
    for seg in sec:
        h.setpointer(seg._ref_i_membrane_, 'im', seg.xtra)
I need to do BOTH of these things for the simulation to run properly--just one or the other is not sufficient. Of course, at that point I'm back to relying on the extracellular mechanism to record the LFP, which is what I'm trying to avoid.

The fact that the segmentation fault occurs in h.finitialize makes me wonder if perhaps there is a problem with my modified xtra.mod file (which I modified according to your instructions here: viewtopic.php?f=2&t=3389&p=14342&hilit= ... lel#p14342), but I have triple-checked this and am pretty sure I did exactly what you said. Or perhaps there is something wrong with the line h.setpointer(seg._ref_i_membrane_, 'im', seg.xtra) (though it appears to my eyes to be in line with the syntax outlined here: viewtopic.php?f=2&t=1685&p=14075&hilit= ... ter#p14075).

Thanks for the help.

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Fri May 05, 2017 9:16 am

All right, I figured out what was causing the segmentation fault in h.finitialize: I wasn't linking the ex pointer within xtra.mod. This post was helpful: viewtopic.php?f=16&t=3403&p=14429&hilit ... ize#p14429

xtra.mod defines two pointers, im and ex. im is used for extracellular recording, and ex is used for extracellular stimulation. From my understanding of this post (viewtopic.php?f=2&t=3389&p=14342&hilit= ... lel#p14342), im no longer requires the 'extracellular' mechanism, but ex still does. This explains why I could only get my simulation to run when I inserted 'extracellular' and used your original setpointers.hoc function: ex still needed them.

But at this point I care only about extracellular recording, and not at all about stimulation. So in order to do away with my need for the 'extracellular' mechanism (and all its computational overhead), I believe I need to comment out all the parts of xtra.mod involving ex. When I do this, re-compile xtra.mod, and run the simulation, it runs without a glitch. Is this the correct approach to recording extracellular potentials without the 'extracellular' mechanism? Here is what my xtra.mod looks like now:

Code: Select all

NEURON{
    SUFFIX xtra
    RANGE rx, er
    RANGE x, y, z
:   GLOBAL is
:   POINTER im, ex  : get rid of ex, because it requires the 'extracellular' mechanism, which we do not want to use
    POINTER im
}

PARAMETER{
    rx =1 (megohm) : mV/nA
    x = 0 (1) :spatial coords
    y = 0 (1)
    z = 0 (1)
}

ASSIGNED {
    v (millivolts)
:   is (milliamp)   : get rid of is, becasue it is only used for extracellular stimulation, which we are not interested in
:   ex (millivolts)  : again, get rid of ex because we don't want to use the extracellular mechanism
:   im (milliamp/cm2)  : see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3389&p=14342&hilit=extracellular+recording+parallel#p14342
    im (nA)
    er (microvolts)
:   area (micron2)  : see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3389&p=14342&hilit=extracellular+recording+parallel#p14342
}

INITIAL {
:   ex = is*rx*(1e6)  : getting rid of ex, as described above
:   er = (10)*rx*im*area  : see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3389&p=14342&hilit=extracellular+recording+parallel#p14342
    er = rx*im
}

:BEFORE BREAKPOINT {
:    ex = is*rx*(1e6)  : getting rid of ex, as described above
:}

AFTER SOLVE {
:   er = (10)*rx*im*area  : see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3389&p=14342&hilit=extracellular+recording+parallel#p14342
    er = rx*im
}
Although this gets the simulation to run, it is unfortunately not the whole story, because the LFP it produces is simply flat-lined, regardless of what the cell is actually doing. I added a synaptic event to the simulation, and it gave the following result for the intracellular voltage (top figure) and the extracellular potential (bottom figure):

Image

Any idea what the problem might be? Am I doing something wrong in my revisions to xtra.mod?

ted
Site Admin
Posts: 5742
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by ted » Fri May 05, 2017 5:08 pm

Your NMODL revisions look good. You checked with modlunit to make sure that units are compatible?

The variable of interest is i_membrane_ (notice that terminal underscore character). This is a range variable (each segment will have its own i_membrane_). But i_membrane_ is calculated only if you call the CVode class's use_fast_imem method--read about this in the Programmer's Reference
https://www.neuron.yale.edu/neuron/stat ... _fast_imem
Sorry that's so buried in the documentation.

Also, I know the documentation says that i_membrane_ will exist at even the 0 and 1 nodes, but its value at those locations should be 0--those nodes are not associated with membrane, so they shouldn't produce transmembrane current. But that's of no consequence since you're going to use the
for (x,0)
syntax to iterate only over internal nodes when you set up the pointers.

Note that an instance of the CVode class is created automatically if you are starting NEURON from the command line via nrngui--that is
nrngui name_of_your_program.hoc
or
nrngui -python progname.py
or if your program contains a
load_file("nrngui.hoc")
statement, or even a
load_file("noload.hoc")
or
load_file("stdgui.hoc")
statement. The name of that instance of the CVode class will be
CVode[0]
so all you have to do is execute the statement
CVode[0].use_fast_imem(1)
(do that after the load_file("nrngui.hoc") etc. statement).

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Sat May 06, 2017 6:15 pm

Thanks, Ted. Good point on using modlunit. When I did so, I found that I needed to use 'nanoamps' instead of 'nA'. Also, I needed to change er = rx * im to er = (1000) * rx * im, which makes sense since er is in microvolts, rx is in megohms, and im is in nanoamps. (Although, weirdly, modlunit continued to tell me to make this change even after I made it.)

And thanks for the tip with use_fast_imem. I knew about this before, but apparently coded it slightly incorrectly (I had been using h('cvode.use_fast_imem(1)') ). Unfortunately, even with your code, I still get that the LFP is flat-lined throughout the entire simulation (including a synaptic event) . And even though my model is simple, it has two different sections, so that's not the problem. Even more strange is when I set nseg=5 for both the soma and the dendrite, simulation results are not consistent for the LFP. The voltage traces are always the same, but the value the LFP is flat-lined at changes from simulation to simulation (with the value sometimes being as large as 10^228...!), as shown in the two figures below, taken from successive simulations:

Image

Image

Is it possible there is a bug associated with use_fast_imem or i_membrane_ ?

Here is the code that produced the above figures...

Code: Select all

from __future__ import division
from neuron import h, gui
import matplotlib.pyplot as plt 

h('CVode[0].use_fast_imem(1)') #see use_fast_imem() at https://www.neuron.yale.edu/neuron/static/new_doc/simctrl/cvode.html

soma=h.Section(name='soma') 
Bdend=h.Section(name='Bdend')

soma.diam = soma.L = 20 #microns
soma.insert('hh')
soma.nseg=5
Bdend.diam = 2 #microns
Bdend.L=200 #microns
Bdend.insert('pas')
Bdend.nseg=5

Bdend.connect(soma,0,0) 

for sec in h.allsec():
    sec.insert('xtra')
    
'''now define spatial location of the cell'''
h.define_shape() #this assigns default position of cell, with soma at origin
#check to make sure h.define_shape() did what we thought it should do
for sec in h.allsec():
    print "section=",sec.name() #or use print h.secname()
    for i in xrange(int(h.n3d())): print i, h.x3d(i), h.y3d(i), h.z3d(i), h.diam3d(i)
    
#call grindaway() in order to calculate centers of all segements, and--most importantly--copy these values to the xtra mechanism of each cell, so they can be used to calculate the LFP
h('xopen("interpxyz.hoc")') #from Ted's extracellular_stim_and_rec code; see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=8&t=3649
h('grindaway()') 
#set pointers
for sec in h.allsec():
    if h.ismembrane('xtra'):
        for seg in sec:
            h.setpointer(seg._ref_i_membrane_, 'im', seg.xtra)
    
#this computes transfer resistances between electrode and all segements;
#note this is a slightly altered version from Ted's original calcrxc.hoc, since we don't need the GUI code
h('xopen("calcrxc_a.hoc")')
h('setelec(50,0,0)') #this is what actually sets electrode location and computes transfer resistances'''

# add synapse
syn = h.Exp2Syn(Bdend(0.5)) #attach synapse to designated section(x)
syn.e = 0 #set equilibrium potential to make synapse excitatory

#add a netcon to the synapse https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3596
#netstim might also be helpful: http://neuron.yale.edu/neuron/static/docs/neuronpython/ballandstick3.html#make-a-netstim
h("objref nil")
nil=h.nil
nc=h.NetCon(nil,syn)
nc.weight[0]=0.002
# this results in the event being created at some point within finitialize(), in particular AFTER all the events have been cleared (see p. 185 of TNB)
fih = h.FInitializeHandler('nrnpython("nc.event(100)")')

#record soma voltage and time
t_vec = h.Vector()
v_vec_soma = h.Vector()
v_vec_dend = h.Vector()
v_vec_soma.record(soma(0.5)._ref_v) 
v_vec_dend.record(Bdend(0.5)._ref_v)
t_vec.record(h._ref_t)

# set up recording of LFP at every time step, in conjunction with Hines's beforestep_py.mod
s = h.Section()
bscallback = h.beforestep_callback(s(.5))

lfp_rec=[]
def callback():
    vx = 0
    for sec in h.allsec():
        if h.ismembrane('xtra'):
            for seg in sec:
                #print h.t, seg.er_xtra
                vx = vx + seg.er_xtra
            #end for seg in sec
        #end if h.ismembrane('xtra')
    #end for sec in h.allsec()
    lfp_rec.append(vx)
    #print h.t, vx

bscallback.set_callback(callback)

#run simulation
h.tstop = 200 #ms
h.run()

#plot voltage versus time
plt.figure(figsize=(8,4))
plt.plot(t_vec,v_vec_soma,label='soma')
plt.plot(t_vec,v_vec_dend,'r',label='dendrite')
plt.xlabel('Time (ms)')
plt.ylabel('mV')
plt.legend(loc='upper left')
plt.show()

#plot LFP versus time
plt.figure(figsize=(8,4))
plt.plot(t_vec,lfp_rec)
plt.xlabel('Time (ms)')
plt.ylabel('LFP')
plt.show()

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Sun May 07, 2017 12:10 pm

Hmmm...I solved the reproducibility problem by commenting out from neuron import gui and replacing it with h('load_file("nrngui.hoc")'). Apparently these two statements are not identical?

In any case, I still have the flat-line problem. Here is the result I consistently get:

Image

pascal
Posts: 89
Joined: Thu Apr 26, 2012 11:51 am

Re: Getting the 'xtra' mechanism to work without 'extracellu

Post by pascal » Fri Jun 09, 2017 7:11 pm

All right, finally figured out the solution. Much thanks to Ted and Robert for their help with this. For some reason beforestep_py.mod messes with i_membrane_'s location in memory, and this problem can be avoided by instead using cvode's extra_scatter_gather method to call the function which computes the LFP at every time step. Here is the Python code:

Code: Select all

from __future__ import division
from neuron import h#, gui
import matplotlib.pyplot as plt 
import numpy as np

h.load_file("stdgui.hoc") #for some reason, need this instead of 'from neuron import gui' to get the simulation to be reproducible, and for the LFP not to be flat-lined
h.cvode.use_fast_imem(1) #see use_fast_imem() at https://www.neuron.yale.edu/neuron/static/new_doc/simctrl/cvode.html

soma=h.Section(name='soma') 
Bdend=h.Section(name='Bdend')

soma.diam = soma.L = 20 #microns
soma.nseg=5
soma.insert('hh')
Bdend.diam = 2 #microns
Bdend.L=200 #microns
Bdend.nseg=5
Bdend.insert('pas')

Bdend.connect(soma,0,0) 

for sec in h.allsec():
    sec.insert('xtra')
    
h.define_shape() #this assigns default position of cell, with soma at origin
    
#call grindaway() in order to calculate centers of all segements, and--most importantly--copy these values to the xtra mechanism of each cell, so they can be used to calculate the LFP
h.load_file('interpxyz.hoc') #from Ted's extracellular_stim_and_rec code; see https://www.neuron.yale.edu/phpBB/viewtopic.php?f=8&t=3649
h.grindaway()
#set pointers
for sec in h.allsec():
    if h.ismembrane('xtra'):
        for seg in sec:
            h.setpointer(seg._ref_i_membrane_, 'im', seg.xtra)
    
#this computes transfer resistances between electrode and all segements;
#note this is a slightly altered version from Ted's original calcrxc.hoc, since we don't need the GUI code
h.load_file("calcrxc_a.hoc")
h.setelec(50,0,0) #this is what actually sets electrode location and computes transfer resistances

# add synapse
syn = h.Exp2Syn(Bdend(0.5)) #attach synapse to designated section(x)
syn.e = 0 #set equilibrium potential to make synapse excitatory

#add a netcon to the synapse https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3596
nc=h.NetCon(None,syn)
nc.weight[0]=0.002
# this results in the event being created at some point within finitialize(), in particular AFTER all the events have been cleared (see p. 185 of The Neuron Book)
fih = h.FInitializeHandler((nc.event,100)) #passes a tuple to FInitializeHandler: a function (in this case nc.event), and an argument to the function (in this case 100)

#record soma voltage and time
t_vec = h.Vector()
v_vec_soma = h.Vector()
v_vec_dend = h.Vector()
v_vec_soma.record(soma(0.5)._ref_v) 
v_vec_dend.record(Bdend(0.5)._ref_v)
t_vec.record(h._ref_t)

imem_soma_vec=h.Vector()
imem_soma_vec.record(soma(0.5)._ref_i_membrane_)

imem_dend_vec=h.Vector()
imem_dend_vec.record(Bdend(0.5)._ref_i_membrane_)

# set up recording of LFP at every time step
lfp_rec=[]
def callback():
    vx = 0
    for sec in h.allsec():
        if h.ismembrane('xtra'):
            for seg in sec:
                vx = vx + seg.er_xtra
    lfp_rec.append(vx)

h.cvode.extra_scatter_gather(0,callback) #instead of using mod file mechanism for callback, use extra_scatter_gather

#run simulation
h.tstop = 200 #ms
h.run()

#plot voltage versus time
plt.figure(figsize=(8,4))
plt.plot(t_vec,v_vec_soma,label='soma')
plt.plot(t_vec,v_vec_dend,'r',label='dendrite')
plt.xlabel('Time (ms)')
plt.ylabel('mV')
plt.legend(loc='upper left')

t_vec_np=t_vec.as_numpy() #allows you to treat hoc Vector as nunpy array, and it doesn't take any more memory! Both t_vec and t_vec_np points to the same slots in memory
#plot LFP versus time
plt.figure(figsize=(8,4))
plt.plot(t_vec_np[0:-1],lfp_rec,marker='.') #throw out the last time point in t_vec, since for some reason using extra_scatter_gather doesn't get the last time point for the LFP (I think it calls the callback function before advance is called)
plt.xlabel('Time (ms)')
plt.ylabel('LFP')

#plot i_membrane_ versus time
plt.figure(figsize=(8,4))
plt.plot(t_vec,imem_soma_vec,label='soma')
plt.plot(t_vec,imem_dend_vec,label='dend')
plt.legend
plt.xlabel('Time (ms)')
plt.ylabel('i_membrane_')
plt.show()
And here are the ancillary files:

Code: Select all

//calcrxc_a.hoc, adapted from Ted's calcrxc.hoc
rho = 351    // for living mammalian cells, change to brain tissue or Ringer's value  351 from brain resistivity paper

// assume monopolar stimulation and recording
// electrode coordinates:
// for this test, default location is 50 microns horizontal from the cell's 0,0,0

XE = -50  // um
YE = 0
ZE = 0

proc setrx() {  // now expects xyc coords as arguments
  forall {
    if (ismembrane("xtra")) {
// avoid nodes at 0 and 1 ends, so as not to override values at internal nodes
      for (x,0) {
//        r = sqrt((x_xtra(x) - xe)^2 + (y_xtra(x) - ye)^2 + (z_xtra(x) - ze)^2)
        r = sqrt((x_xtra(x) - $1)^2 + (y_xtra(x) - $2)^2 + (z_xtra(x) - $3)^2)
        // 0.01 converts rho's cm to um and ohm to megohm
        rx_xtra(x) = (rho / 4 / PI)*(1/r)*0.01
	//print secname(), "  \t", rx_xtra(x)
      }
    }
  }
}

proc setelec() {
	xe = $1
	ye = $2
	ze = $3
	setrx(xe, ye, ze)
}

Code: Select all

// $Id: interpxyz.hoc,v 1.2 2005/09/10 23:02:15 ted Exp $
/* Computes xyz coords of nodes in a model cell 
   whose topology & geometry are defined by pt3d data.
   Expects sections to already exist, and that the xtra mechanism has been inserted
 */


// original data, irregularly spaced
objref xx, yy, zz, length
// interpolated data, spaced at regular intervals
objref xint, yint, zint, range

proc grindaway() { local ii, nn, kk, xr
	forall {
	  if (ismembrane("xtra")) {
		// get the data for the section
		nn = n3d()
		xx = new Vector(nn)
		yy = new Vector(nn)
		zz = new Vector(nn)
		length = new Vector(nn) 

		for ii = 0,nn-1 {
			xx.x[ii] = x3d(ii)
			yy.x[ii] = y3d(ii)
			zz.x[ii] = z3d(ii)
			length.x[ii] = arc3d(ii)
		}

		// to use Vector class's .interpolate() 
		// must first scale the independent variable
		// i.e. normalize length along centroid
		length.div(length.x[nn-1])

		// initialize the destination "independent" vector
		range = new Vector(nseg+2)
		range.indgen(1/nseg)
		range.sub(1/(2*nseg))
		range.x[0]=0
		range.x[nseg+1]=1

		// length contains the normalized distances of the pt3d points 
		// along the centroid of the section.  These are spaced at 
		// irregular intervals.
		// range contains the normalized distances of the nodes along the 
		// centroid of the section.  These are spaced at regular intervals.
		// Ready to interpolate.

		xint = new Vector(nseg+2)
		yint = new Vector(nseg+2)
		zint = new Vector(nseg+2)
		xint.interpolate(range, length, xx)
		yint.interpolate(range, length, yy)
		zint.interpolate(range, length, zz)

		// for each node, assign the xyz values to x_xtra, y_xtra, z_xtra
//		for ii = 0, nseg+1 {
// don't bother computing coords of the 0 and 1 ends
// also avoid writing coords of the 1 end into the last internal node's coords
		for ii = 1, nseg {
			xr = range.x[ii]
			x_xtra(xr) = xint.x[ii] //x_xtra is a range variable, and xr is a normalized distance, x_xtra(xr) is the value of x_xtra in the compartment whose center is nearest the noramlized distance xr
			y_xtra(xr) = yint.x[ii]
			z_xtra(xr) = zint.x[ii]
		}
	  }
	}
}

Code: Select all

: $Id: modified version of Ted's xtra.mod, which allows for extracellular recording without the extracellular mechanism $

COMMENT
Now uses i_membrane_ (see CVode class's use_fast_imem documentation)
See  https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3389&p=14342&hilit=extracellular+recording+parallel#p14342
ENDCOMMENT

NEURON {
	SUFFIX xtra
	RANGE rx, er
	RANGE x, y, z
:	GLOBAL is
:	POINTER im, ex
	POINTER im
}

PARAMETER {
	: default transfer resistance between stim electrodes and axon
	rx = 1 (megohm) : mV/nA
	x = 0 (1) : spatial coords
	y = 0 (1)
	z = 0 (1)
}

ASSIGNED {
	v (millivolts)
:	is (milliamp)
:	ex (millivolts)
:	im (milliamp/cm2)
	im (nanoamp)
	er (microvolts)
:	area (micron2)
}

INITIAL {
:	ex = is*rx*(1e6)
:	er = (10)*rx*im*area
	er = (1000)*rx*im
: this demonstrates that area is known
: UNITSOFF
: printf("area = %f\n", area)
: UNITSON
}

: Use BREAKPOINT for NEURON 5.4 and earlier
: BREAKPOINT {
:	SOLVE f METHOD cvode_t
: }
:
: PROCEDURE f() {
:	: 1 mA * 1 megohm is 1000 volts
:	: but ex is in mV
:	ex = is*rx*(1e6)
:	er = (10)*rx*im*area
: }

: With NEURON 5.5 and later, abandon the BREAKPOINT block and PROCEDURE f(),
: and instead use BEFORE BREAKPOINT and AFTER BREAKPOINT

:BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
:  ex = is*rx*(1e6)
:}
AFTER SOLVE { : after each solution step
:  er = (10)*rx*im*area
   er = (1000)*rx*im
}
I implemented this approach in a large-scale network simulation, and it gave approximately a 25% speed-up over recording the LFP using 'xtra' in combination with 'extracellular'.

Post Reply