Utilizing extracellular in a for loop

The basics of how to develop, test, and use models.
Post Reply
tim
Posts: 3
Joined: Wed Dec 12, 2007 2:19 am

Utilizing extracellular in a for loop

Post by tim »

First off, here's what I'm trying to do. I'm trying to model an extracellular monopolar current stimulus and apply it to an axon with n nodes/segments (currently set at 21 for simplicity). I have MATLAB code which models the extracellular monopolar current source and calculates corresponding voltages at the membrane for each of the 21 nodes. That code has been vetted and creates a file with column one being time, and columns 2 through n+1 being nodal voltages at each time step.

The NEURON code reads in that file (as a matrix), and then is supposed to use a for loop to strip out each column 2 through n+1 in sequence as a vector and then play it via the extracellular mechanism. However, when the code is run, it only applies the extracellular voltage stimulus at the most distal node (i.e. at axon position 1.0). I (and the person I'm working with) haven't been able to figure out the problem but think it may have something to do with handling of pointers within the loop. I've been able to play vectors fine outside of the loop and the loop itself is iterating as it should.

Any thoughts on what we're missing?

Thanks very much for the help.

*****

load_file("nrngui.hoc")

create axon
access axon

// Global simulation constants
num_nodes = 21 // compartments
Dt = 0.1 // msec
tstop = 500 // msec

// Constants for active axon with HH channels
axon nseg = 100
axon diam = 1
axon L = 28000
axon Ra = 100
axon insert hh
axon insert pas
g_pas = 1/20000
e_pas = -60

// Insert extracellular mechanism. Changed nrn/src/nrnoc/options.h line #define EXTRACELLULAR 1
// Syntax:
// insert extracellular
// vext[2] -- mV
// i_membrane -- mA/cm2
// xraxial[2] -- MOhms/cm
// xg[2] -- mho/cm2
// xc[2] -- uF/cm2
// e_extracellular -- mV
axon insert extracellular

// Apply extracellular stim from MATLAB file
objref f1, v1, t1, m1
f1=new File()
num_rows = tstop/Dt+1
num_cols = num_nodes+1
m1 = new Matrix(num_rows,num_cols)
v1=new Vector(num_rows)
t1=new Vector(num_rows)
for i=1,num_nodes{
print "\nLoad compartment ", i
f1.ropen("HFAC_extra_21node_10Hz.dat")
m1.scanf(f1,num_rows,num_cols)
v1=m1.getcol(i)
// Only strip out the time vector once
if (i==1) {
t1=m1.getcol(i-1)
}
norm_step = (i-1)/(num_nodes-1)
print "Axon fraction: ", norm_step
v1.play(&e_extracellular(norm_step), t1)
}

// Current injection at proximal end of cable (to produce APs)
objectvar stim
axon stim = new IClamp(0)
stim.del = 10
stim.dur = 100
stim.amp = 1

*****

Too bad phBB won't preserve tab indents...
ted
Site Admin
Posts: 5869
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

Too bad phBB won't preserve tab indents...
Wrap your code block inside the appropriate markup.
left_square_bracket code right_square_bracket
all your code
left_square_bracket /code right_square_bracket
where
left_square_bracket means [
and
right_square_bracket means ]
Or just click and drag to select all your code, then click on the "Code" button at the top of
the "Message body" area on the "Post a message" page.

Now about your program--
it only applies the extracellular voltage stimulus at the most distal node (i.e. at axon position 1.0)
You want to drive each node with a different waveform, so you must have a different
Vector for each node. But your code tries to reuse the same Vector for each node.
Won't work.

Other comments:

The nodes at 0 and 1 are not associated with any membrane, so you can skip them.

It's risky to rely on code like
num_rows = tstop/Dt+1
because it is subject to roundoff error. For the task of reading a matrix of data from a file,
a safe and very useful strategy is to embed the matrix size as the file header, i.e.
nrow ncol datum0 datum1 . . .
Then the Matrix class's scanf can be used to read the whole thing in one pass.
tim
Posts: 3
Joined: Wed Dec 12, 2007 2:19 am

Post by tim »

ted wrote:
Too bad phBB won't preserve tab indents...
Wrap your code block inside the appropriate markup.
Thanks! I stand corrected.
it only applies the extracellular voltage stimulus at the most distal node (i.e. at axon position 1.0)
You want to drive each node with a different waveform, so you must have a different
Vector for each node. But your code tries to reuse the same Vector for each node.
Won't work.
That explains it. Our model will eventually be expanded from 21 to over 100 nodes, so my understanding is I could use list.append and list.object within a loop to create n vectors and then fill them with values from the appropriate columns of the data file. Am I on the right track with this?
Other comments:

The nodes at 0 and 1 are not associated with any membrane, so you can skip them.

It's risky to rely on code like
num_rows = tstop/Dt+1
because it is subject to roundoff error. For the task of reading a matrix of data from a file,
a safe and very useful strategy is to embed the matrix size as the file header, i.e.
nrow ncol datum0 datum1 . . .
Then the Matrix class's scanf can be used to read the whole thing in one pass.
Thanks for these suggestions. More things for the revision list. :)
ted
Site Admin
Posts: 5869
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

Is there a simple way to generate iteratively a large number of vectors (e.g. in a for loop)?
The most general and flexible way to deal with multiple items is with lists. Assuming a
Matrix m with at least nseg+1 columns,

Code: Select all

objref tvec, vvec, veclist
veclist = new List()
tvec = new Vector()
tvec = m.getcol(0)

ii = 0
for (x, 0) { // avoid zero area nodes
  ii += 1
  vvec = new Vector() // create a new Vector
    // and, after the first pass, destroy one reference to the previous Vector
  vvec = m.getcol(ii)
  vvec.play(&rangevar(x), tvec)
  veclist.append(vvec) // create a 2nd reference to the Vector
    // so the data will persist even after vvec = new Vector()
    // or objref vvec
}

objref vvec // so careless reuse of vvec won't step on the last data
Lists have many advantages over "arrays of objrefs," not least the fact that they
encourage writing code that doesn't rely on "magic numbers."
tim
Posts: 3
Joined: Wed Dec 12, 2007 2:19 am

Post by tim »

^ Implementation using a list to handle multiple vectors (in a for loop) has worked like a charm. Thanks for the help!
ted
Site Admin
Posts: 5869
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

If your work with NEURON generates publications, please be sure to cite the NEURON book
and/or one or more of the papers about NEURON, as appropriate (see
https://www.neuron.yale.edu/phpBB2/viewtopic.php?p=121 ). And just as important, please
tell us about your publications so we can include them in the Bibliography page
http://www.neuron.yale.edu/neuron/bib/usednrn.html
ny

Post by ny »

Hi, I am working on something similar - stimulating an axon extracellularly. Can you tell me why you need an IClamp "to produce APs"?
Thanks
ted
Site Admin
Posts: 5869
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

ny wrote:I am working on something similar - stimulating an axon extracellularly. Can you tell me why you need an IClamp "to produce APs"?
Unlike you, that implementer used intracellular stimulation to elicit spikes.
ny

Post by ny »

I see, so I should be able to elicit spikes just by setting e_extracellular?
ted
Site Admin
Posts: 5869
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

That can be done with NEURON. See the item
Extracellular stimulation and recording
in the Hot tips
section of The NEURON Forum.
ny

Post by ny »

I am currently doing this by setting e_extracellular via the vector play method, and using the results from a finite element model of electric field to provide the potential values. But I don't see any spikes fired by the axon, though I do observe changes in the membrane potential at the right firing frequency indicating that the method is working.
Post Reply