Miscellaneous items before getting to your question--
Extracellular field cannot change the transmembrane potential of a single compartment
model. There's no return path for the current. Must have at least two compartments (in
NEURON should use odd nseg, so either the model must consist of two sections attached
to each other, both with nseg == 1, or one section with nseg == 3, 5, or any other odd number).
Locations between 0 and 1 are mapped to the nearest segment center. With nseg == 1,
this means that all rangevar(x), where 0<x<1, actually refer to rangevar(0.5)
Why Vector play() isn't working
You're not getting Vector play() to work because Vector play() depends on code that is
part of the standard run system, so you need to know how to take advantage of the
standard run system's features.
First rule: don't break the standard run system. Generally that means use the standard
run system without modification, and in particular, don't replace proc run(). Only in
special situations is it ever necessary to even customize the standard run system, and
this isn't one of them. If and when that rare occasion does arise, be sure to read
chapters 7 and 8 of The NEURON Book first.
An aside: A lot of legacy code is floating around that was written by very productive
modelers, who got started well before the standard run system reached its present
level of sophistication. These guys are still churning out lots of code that works, and it
works well--for them. They can make it work, but it tends to be idiosyncratic, and often
doesn't play well with NEURON's latest tools and features.
Two strong recommendations:
1. "Legacy programming styles" are best avoided by people who are getting started.
Otherwise a lot of time and effort are wasted in trying to duplicate already existing
functionality.
2. Use the GUI as much as possible. Especially for setting up graphs. The GUi is very
good at automatically generating a lot of "boilerplate" code that is needed to set up
graphs and properly hook them into the run control system. The tutorials at
http://www.neuron.yale.edu/neuron/docs
show how to use the GUI for this and a lot more sophisticated tasks too.
Back to your particular modeling project:
Here's how I would organize a program that implements extracellular stimulation.
1. specify anatomical and biophysical properties of model cell
2. initialize Vectors that describe the time course of the extracellular field
3. attach the Vectors to the variables that they will drive during a simulation
4. set up controls for launching simulations and displaying results
Item 4 would use the standard run system without modification. No need to use g.begin,
g.beginline, g.line, g.plot, or call fadvance. The standard run system will take care of
graph initialization and updating for you.
There is no need to execute any part of step 1, 2, or 3 on each simulation run. That's all
"setup" code, which only has to be executed once.
Here's how I did it, through a process of incremental revision and testing.
Pass 1
File 1: init.hoc
Initially it contained just
Code: Select all
load_file("nrngui.hoc")
load_file("toymodel.hoc")
FIle 2: toymodel.hoc
Code: Select all
create axon
access axon
nseg = 3
L = 100
diam = 1
insert hh
Just enough detail to get started. hh so it can spike. nseg = 3 so it can be excited by
extracellular field.
Then I used NEURON to execute init.hoc, and used the NEURON Main Menu toolbar to
set up:
--a RunControl panel to launch simulations, but mostly for convenient control over tstop,
dt, etc.
--a Voltage axis graph (that automatically plots v(.5))
--a Shape plot that I immediately turned into a space plot of v
--a Movie Run panel so I could see the space plot evolve smoothly with time
After arranging these nicely so they didn't overlap, I saved them to a file called rig.ses
Next I revised init.hoc to read
Code: Select all
load_file("nrngui.hoc")
load_file("toymodel.hoc")
load_file("rig.ses")
and verified that executing init.hoc would recreate my custom user interface (RunControl
and graphs etc.).
Pass 2
Prepare to set up the stimulus Vectors. I deferred file reading. For now, it was expedient
to merely fill some vectors with a simple rectangular pulse.
So I created a file called stim.hoc, with these contents:
Code: Select all
// create the basic stimulus time course
objref tvec, pvec
// tvec will hold the stimulus sample times
// pvec will be a 1 ms pulse of with amplitude PMAX
NUMPTS = 5
tvec = new Vector(NUMPTS)
pvec = new Vector(NUMPTS)
PMAX = 40 // found by trial and error to elicit a spike
{ tvec.x[0]=0 pvec.x[0]=0 } // field is 0 for 1 ms
{ tvec.x[1]=1 pvec.x[1]=0 }
{ tvec.x[2]=1 pvec.x[2]=PMAX } // jumps to some value
{ tvec.x[3]=2 pvec.x[3]=PMAX } // for 1 ms
{ tvec.x[4]=2 pvec.x[4]=0 } // falls back to 0 ever after
print " "
print "tvec"
tvec.printf
print " "
print "pvec"
pvec.printf
Also, I added this line to the end of init.hoc:
Now running init.hoc recreated my GUI and printed out the vectors that describe the time
course of the basic stimulus waveform.
Pass 3
Actually create the stimulus vectors and attach them to the dependent variables.
I commented out the
Code: Select all
print " "
print "tvec"
tvec.printf
print " "
print "pvec"
pvec.printf
stuff at the end of stim.hoc, and added this:
Code: Select all
// drive e_extracellular at each internal node of the model
// with a voltage that has the time course specified by pvec, tvec
forall insert extracellular
objref veclist // will hold all the stim Vectors
proc setstim() { localobj tmpvec
veclist = new List()
forall {
for (x, 0) { // iterate over internal nodes only
// specify the time course of extracellular potential at this location
// for this toy example, something very simple:
tmpvec = pvec.c
tmpvec.mul(1-x) // potential falls off with distance from 0 end
// could just as well have read values from a file
// or copied a column from a matrix into tmpvec
tmpvec.play(&e_extracellular(x), tvec)
veclist.append(tmpvec)
}
}
}
setstim()
Now running init.hoc and clicking on the RunControl panel's Init & Run elicited a spike.
Clicking on the Movie Run panel's Init & Run, I could see the spike propagate from one
end ot the axon to the other.
Pass 4
Trust, but verify. I wanted to see how e_extracellular varied along the length of the axon.
So I brought up another Shape plot, used its Plot what? menu item to bring up a variable
name browser, selected
e_extracellular
from the list, and then turned this graph into a Space plot of e_extracellular.
Now, clicking on the Movie Run panel's Init & Run button showed the linear decay of
e_extracellular with distance along the axon, during the stimulus current.
At this point, I used the Print & File Window Manager to save the e_extracellular Graph
to a file called plot_e_extracellular.ses,
and then added this line to the end of init.hoc:
Then I tested it by exiting NEURON again and running init.hoc once more.
The whole thing is in
http://www.neuron.yale.edu/ftp/ted/neuron/ecstim.zip