Extracellular Fields

Managing anatomically complex model cells with the CellBuilder. Importing morphometric data with NEURON's Import3D tool or Robert Cannon's CVAPP. Where to find detailed morphometric data.
Der Ted and Coaster.
This had been a really great topic. I downloaded the ecstim.zip and played around with it. After that I still have two questions which I hope you could explain to me.

1) Simply changing the axon length to 5 cm in the toymodel.hoc file results in a flat line (-65mV) in the voltage plot. Using smaler values shows that the voltage effect decreases with increased length. Why does the lenght of the axon has that effect?

create axon
access axon
nseg = 3
L = 50000
diam = 1
insert hh

2) Changing the stim.hoc file so that there are two stimulations does not change my voltage plot.

NUMPTS = 10
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
{ tvec.x[5]=3 pvec.x[5]=0 } // falls back to 0 ever after
{ tvec.x[6]=3 pvec.x[6]=PMAX } // jumps to some value
{ tvec.x[7]=4 pvec.x[7]=PMAX } // f for 1 ms
{ tvec.x[8]=4 pvec.x[8]=0 } // falls back to 0 ever after
{ tvec.x[9]=5 pvec.x[9]=0 } // falls back to 0 ever after

I hope you could help me understand both of the questions.
Thank you.
nseg

Posts: 3
Joined: Mon Jun 09, 2008 10:05 pm

1) Simply changing the axon length to 5 cm in the toymodel.hoc file results in a flat line (-65mV) in the voltage plot. Using smaler values shows that the voltage effect decreases with increased length. Why does the lenght of the axon has that effect?

Increasing L while leaving nseg unchanged has three principal effects that reduce the
efficacy of a stimulus in a computational model of extracellular stimulation.
--increasing the area of each compartment (increases compartment capacitance, so that
more charge must accumulate in order to depolarize the membrane)
--increasing resistance between adjacent nodes (interferes with the axial flow of current
which is necessary for membrane at opposite ends of the neurite to depolarize or
hyperpolarize, depending on orientation of the fiber and the extracellular field
--"averaging out" local maxima of the extracellular field, which are important in the near
vicinity of small electrodes

Put another way, nseg is far too small for spatial accuracy in such an anatomically
extended model. Use the d_lambda rule to figure out how long it should be. nseg may
have to be larger than that for neurites that lie within the "near field" of monopolar or
multipolar electrode configurations.

2) Changing the stim.hoc file so that there are two stimulations does not change my voltage plot.

Under what conditions? L = 50000?
Last edited by ted on Wed Jun 11, 2008 10:31 am, edited 1 time in total.
ted

Posts: 3591
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Dear Ted,

Put another way, nseg is far too small for spatial accuracy in such an anatomically
extended model. Use the d_lambda rule to figure out how long it should be.

nseg:
nseg=axon_length/(d_lambda*lambda_100)

With:
d_lambda=0.1
Ra=110
Cm=2.5
diam=0.72
axon_length=50000

lambda_100=1e5*sqrt(diam/(4*pi*100*Ra*Cm))=144.34

=> nseg=3465

I used the new value for nseg but the output did not change. After that I played around with the nseg value and used lots of different values to see if I would get any response but no. Could you tell me what went wrong?

Under what conditions? L = 50000?

I used the starting values from the example for the test. L=100 and nseg=3.
nseg

Posts: 3
Joined: Mon Jun 09, 2008 10:05 pm

Because of the way the extracellular field is represented in this particular example, there is
a fourth factor that affects the efficacy of the extracellular stimulus--a factor that depends on
L and has nothing to do with nseg.

In the original example L was 100 um, and the extracellular field was a uniform gradient with
an intensity of 40 mV/100um = 0.4 mV/um.

With the code as written, if L is changed, the extracellular potential difference from one end
to the other remains unchanged (40 mV). So increasing L 500-fold reduces the field intensity
500-fold. The field still affects membrane potential, but the effect is very small. Run a
simulation for 2 ms, then use the graphs' View = plot to see.
ted

Posts: 3591
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Dear Ted,
thank you for your help. I am really starting to understand it much better. I hope you don't mind me aksing more questions until I solved my problem.

In the original example L was 100 um, and the extracellular field was a uniform gradient with
an intensity of 40 mV/100um = 0.4 mV/um.

The intensity is the power supplied divided by the surface area. In our example programm the extracellular stimullation is applied to all nodes.

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)
}
}
}

So I thought, that the intensity loss wouldn't be a problem because now we changed nseg proportional to the length.
I changed the PMAX to 4000 and added the axon.v(0), axon.v(0.1), axon.v(0.2), axon.v(0.3) ... and axon.v(1) to the plot. Axon.v(1) shows a peak similar to the original one, axon.v(0) is a smaller negativ peak and all the others are a straight line.

<img src='http://www.sahler.org/ruth/ecstim-example.jpg'>

The original plot with L=100, nseg=3 and PMAX=40

<img src='http://www.sahler.org/ruth/ecstim-example-100.jpg'>
([ img ] disabled)
nseg

Posts: 3
Joined: Mon Jun 09, 2008 10:05 pm

nseg wrote:The intensity is the power supplied divided by the surface area.

No. The intensity or strength of an electric field is the gradient of electrical potential and
has units of volts/length. If you have a serious interest in extracellular stimulation, you will
have to learn the elementary physics of classical electromagnetism. A one semester course
would be well worth the effort.
n our example programm the extracellular stimullation is applied to all nodes.

Has nothing to do with the fact that 40 mV/100 um is a much stronger field than
40 mV/50000 um.
So I thought, that the intensity loss wouldn't be a problem because now we changed nseg proportional to the length.

enough to give good spatial accuracy, with PMAX == 4, the extracellular field is too weak to
stimulate the axon.
I changed the PMAX to 4000 and added the axon.v(0), axon.v(0.1), axon.v(0.2),
. . .
Axon.v(1) shows a peak similar to the original one, axon.v(0) is a smaller negativ peak and all the others are a straight line.

Good way to investigate the problem. One end of the fiber is strongly depolarized, and
the other is strongly hyperpolarized, because most of the current that enters or leaves
the axon does so near either end. However, the intervening membrane is not isopotential.
If you zoom in you'll discover that, once you get away from either end, there is a gradual
shift in membrane potential along the entire length of the axon. Why? Because most of
the current that enters at one end flows along the entire length of the axon and exits at
the other end. That longitudinal current flow through cytoplasmic resistance results in a
gradual change of internal potential along the length of the axon.
ted

Posts: 3591
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Re: Extracellular Fields

Hello,

New neuron user here that has found this thread and board very helpful. I have been working with the files from extracellular_stim_and_record.zip and I am having trouble with changes I made to stim.hoc which creates the stimulus waveform. I've left all the other files in the package alone (except initxstim.hoc where I've told it to load my stimulus hoc instead). I'm trying to allow for bi-phasic stimulation for now, but NEURON keeps giving me bus errors and I've seen segmentation violation a few times. Below are two bus error messages I often see.

...MacOS/nrniv: Bus error See \$NEURONHOME/lib/help/oc.help
near line 9
{run()}
^
finitialize(-65)
init()
stdinit()
run()

and also

...MacOS/nrniv: Bus error See \$NEURONHOME/lib/help/oc.help
near line 0
{hoc_ac_ = v(hoc_ac_)+vext(hoc_ac_)}
^
step()
continuerun(10)
and others

My stimulus hoc defaults the stimulus to be just like the default stimulus from the original zip package, and my hoc runs and duplicates the results from the original hoc file (as far as I have noticed). Once I start changing the stimulus though, things will either start looking very erratic in the graphs and soon I get an error message from the terminal and NEURON is done working for that session (if that doesn't happen as soon as I start changing the stimulus and hit Init and Run).

Below is the code I've come up with for bi-phasic stimulation. Any ideas as to why I'm crashing NEURON? Thanks in advance for your help and time.
-jkm

Code: Select all
`/* * Current for extracellular stimulation */   // Default valuesCATHODIC = 0               // 1-yes 0-no : is it a cathodic pulse first?BIPHASIC = 0               // 1-yes 0-no : is it biphasic?DELAY = 1                  // delay before application of stimulus (ms)AMPLITUDE = .05               // stimulus amplitude (mA)DURATION = 1               // pulse durations (ms)INTERPULSE = 0               // inter pulse width for biphasic stimulus (ms)objref stim, time, stimGstimG = new Graph(0)   // Make the stimulus and time vectorproc stimWaveform() {   stim = new Vector()   time = new Vector()   j = 0               // What time is it?   Cathodic = \$1   Biphasic = \$2   Delay = \$3   Amplitude = \$4   Duration = \$5   InterPulse = \$6      if (Delay>0) {       for i = 0, Delay {          stim.append(0)         time.append(j)         if ( i != Delay) { j += 1 }         }       }      if (Biphasic) {               if (Cathodic) {                        for i = 0, Duration {             stim.append(Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }                  if (InterPulse>0) {             for i = 0, InterPulse {                stim.append(0)                time.append(j)               if ( i != InterPulse ) { j += 1 }            }          }                  for i = 0, Duration {             stim.append(-Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }               } else {      // anodic first                  for i = 0, Duration {             stim.append(-Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }                  if (InterPulse>0) {             for i = 0, InterPulse {                stim.append(0)                time.append(j)               if ( i != InterPulse ) { j += 1 }            }          }                  for i = 0, Duration {             stim.append(Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }      }         } else {      // monophasic            if (Cathodic) {                        for i = 0, Duration {             stim.append(Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }               } else {                  for i = 0, Duration {             stim.append(-Amplitude)             time.append(j)            if ( i != Duration ) { j += 1 }         }               }   }      stim.append(0)   time.append(j)   }   // is_xtra is GLOBAL, so can play amplitude into itproc attachStim() {   forall {      if (ismembrane("xtra")) {         stim.play(&is_xtra, time)      }   }}proc graphStim() {   stimG.unmap()   stimG.view(0, -AMPLITUDE*1.3, stim.size()*1.3, AMPLITUDE*2.6, 1060, 650, 300, 200)   stim.plot(stimG, time)}proc setStim() {      Cath = \$1   Biph = \$2   Del = \$3   Amp = \$4   Dur = \$5   IntPul = \$6      stimWaveform(Cath, Biph, Del, Amp, Dur, IntPul)   attachStim()   graphStim()}proc go() {   setStim(CATHODIC, BIPHASIC, DELAY, AMPLITUDE, DURATION, INTERPULSE)}go()xpanel("Extracellular Stimulus", 0)   xcheckbox("Cathodic First", &CATHODIC, "go()")   xcheckbox("Biphasic", &BIPHASIC, "go()")   xvalue("Latency (ms)", "DELAY", 1, "go()", 0, 1)   xvalue("Amplitude (mA)", "AMPLITUDE", 1, "go()", 0, 1)   xvalue("Pulse-Width (ms)", "DURATION", 1, "go()", 0, 1)   xvalue("Inter-Pulse-Width (ms)", "INTERPULSE", 1, "go()", 0, 1)xpanel(1052,384)`
jkm

Re: Extracellular Fields

I'm glad you have chosen NEURON for your work, and that you have found the NEURON Forum to be helpful.

That's good use of a custom xpanel to display and control stimulus parameters. Also nice control of the position (and scaling) of the graph.

Before modifying code that works, it is best to identify what needs to be changed and what should be left unaltered. In this particular case, the stuff that needs to be changed is in the file stim.hoc, but not everything in that file needs to be changed. The part that should be left unaltered is the part that attaches the stimulus vector to is_xtra--i.e. everything from
ATTACHED__ = 0
to the end of proc attach_stim(). As per the comment in the original version of attach_stim(), is_xtra is GLOBAL, so it is only necessary to execute the Vector play statement one time.

This raises the side issue of whether it is a good idea to start introducing new variable names for stuff that already exists, e.g. the stimlus Vector stim_amp and its associated time Vector stim_time, procedure names, whatever else. It is usually a bad idea to start renaming things unless there is a compelling reason to do so. For one thing, it is possible that a variable declared and used in one file may be reused somewhere else. For another, it encourages confusion on the part of the programmer. Nontrivial programs often have more names than a Russian novel (even if programs doesn't force one to deal with unfamiliar diminutives, and variables are never called by their patronymics). Better to avoid arbitrary expansion of the name space.

"Well, I didn't want any of my new stuff to conflict with stuff that was in the original stim.hoc."

Fine, but the way to avoid such a conflict is to
1. copy the original initxstim.hoc to initbp.hoc
2. in initbp.hoc change the line
to

Then you can either
(1) copy the original stim.hoc to stimbp.hoc and make your changes in stimbp.hoc,
or, since you have already put some time & effort into your own code, you might want to
(2) rename your own stimulus code file as stimbp.hoc and work on its contents.

Returning to what parts of stim.hoc need to be changed--

The code you present can be easily modified for execution as a standalone program. Just comment out proc attachStim(), and also the line in proc setStim() that calls attachStim(). Then use NEURON to execute it, and play with the latency, pulse width, and interpulse interval.

And you will discover that noninteger values for any of these parameters results in ugly and unanticipated distortions of the stimulus. For example, specify a biphasic stimulus with latency 0.01 ms, pulse width 0.01 ms, and interpulse interval 0.01 ms, and what you get instead is a triangle wave that starts 0 ms, peaks at 1 ms, has an opposite peak at 3 ms, and ends at 4 ms.

proc stimWaveform() is far more complex than it has to be, and it expresses an algorithm that is doomed to fail grotesquely.

It's much easier to steal the proc stim_waveform() from the original stim.hoc and modify it. To see what is necessary, sketch a biphasic waveform and identify each point at which the waveform makes a sudden change of direction. Then place a dot at the origin and imagine another dot at time = 1 ms after the end of the biphasic waveform, stim = 0.* These are the 10 points whose coordinates you must append to the time and stimulus vectors. For a monophasic stimulus, the first 2 and the last 6 points all have stim values of 0. For a biphasic stimulus, the first 2, 5th and 6th, and 9th and 10th have stim values of 0.

You could modify stim_waveform() so that it constructs a complete biphasic waveform, then if BIPHASIC is 0, it sets the 7th and 8th stim values to 0. Or you could write more clever code that appends just enough points to each vector as needed for a monophasic or a biphasic waveform.

*--The last point at time = 1 ms after the end of the biphasic waveform, stim = 0, is necessary because of how NEURON now (version 6.2 and later) deals with interpolated Vector play--see http://www.neuron.yale.edu/neuron/static/ ... .html#play
So you need to give the waveform a flat tail or else NEURON will extrapolate an infinitely large third phase of the stimulus waveform. Talk about grotesque failures.

And to avoid the same, I must now modify the distributed code for stim_waveform() for the same reason--it will become
Code: Select all
`proc stim_waveform() {   // this uses interpolated play  // index    0  1    2    3        4        5  // stim vec 0, 0,   1,   1,       0        0  // time vec 0, DEL, DEL, DEL+DUR, DEL+DUR, DEL+DUR+1  //  really  0, \$1,  \$1,  \$1+\$2,   \$1+\$2,   \$1+\$2+1  // first the stim vector  stim_amp.resize(6)  stim_amp.fill(0)  stim_amp.x[2]=1  stim_amp.x[3]=1  stim_amp.mul(\$3)  // now the time vector  stim_time.resize(6)  stim_time.x[1]=\$1  stim_time.x[2]=\$1  stim_time.x[3]=\$1+\$2  stim_time.x[4]=\$1+\$2  stim_time.x[4]=\$1+\$2+1}`
ted

Posts: 3591
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Previous

Who is online

Users browsing this forum: No registered users and 1 guest