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.
nseg

Post by nseg »

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.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

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.
nseg

Post by nseg »

Dear Ted,
thank you for your replay.
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.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

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.
nseg

Post by nseg »

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.


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

([ img ] disabled)
Last edited by nseg on Thu Oct 23, 2014 10:04 am, edited 1 time in total.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

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.
Maybe my previous message wasn't sufficiently clear about this: even if nseg is large
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.
jkm

Re: Extracellular Fields

Post by jkm »

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_)}
^
fadvance()
advance()
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 values
CATHODIC = 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, stimG
stimG = new Graph(0)


	// Make the stimulus and time vector
proc 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 it
proc 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)

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

Re: Extracellular Fields

Post by ted »

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.

Instead of debugging your code, I will use this occasion to address the task of how to modify code that works.

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
load_file("stim.hoc") // extracellular stimulus
to
// load_file("stim.hoc") // extracellular stimulus
load_file("stimbp.hoc") // biphasic stimulus

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/stati ... .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
}
mzenker

Re:

Post by mzenker »

For people reading this thread: The topic referenced by Ted in the 3rd post is now under viewtopic.php?t=168.

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

Re: Extracellular Fields

Post by ted »

Thanks for noticing the problem with the URL. Lately there has been a lot of tinkering with our server, and something broke. Old URLs that contain phpBB2 should now be working again.
mzenker

Re: Extracellular Fields

Post by mzenker »

Hi,

I have found the information in this thread very useful to set up a model with extracellular stimulation by a potential.
Now I have a simple (?) question.
In the setstim() procedure, there is a vector list veclist which is filled with the stimulus vectors. Apparently this list is just filled and apparently not used anywhere else, so I thought I could remove it. However, if I just comment out the line

Code: Select all

veclist.append(tmpvec)
the stimulation does not work any more.
Why is this? Is veclist used tacitly anywhere? I have tried to find a hint in the documentation for Vector.play() and extracellular, but without success.
What am I missing?

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

Re: Extracellular Fields

Post by ted »

In the setstim() procedure, there is a vector list veclist which is filled with the stimulus vectors.
It has been a long time since I looked at this thread. Where is the post in this thread that contains this veclist.append statement, or if it appears in a hoc file somewhere, where is that file?
mzenker

Re: Extracellular Fields

Post by mzenker »

Hi Ted,

the file is in the zip archive here: http://www.neuron.yale.edu/ftp/ted/neuron/ecstim.zip
The code is shown and described in the 7th post in this thread (I cannot make an exact URL here since it has no subject).

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

Re: Extracellular Fields

Post by ted »

Good question. The files in ecstim.zip illustrate one way to implement a model of extracellular stimulation. In that particular approach, for each internal node (i.e. each segment center) in a model, there is a corresponding Vector whose values are the sequence of extracellular potentials at that location, and that Vector's play() method is used to drive e_extracellular at that node. The number of Vectors required is equal to the number of segments in the model. This works for small models, and indeed it is the only practical way to implement extracellular stimulation if the extracellular medium is dispersive (because of the inherent distance- and frequency-dependent attenuation and phase shifts that result from signal propagation in a dispersive medium). However, it is not memory efficient. More about this below.

Regarding your particular question: in hoc (like Python and many other interpreted languages), garbage collection is controlled by reference counting. proc setstim() appends the Vectors that it creates to an instance of the List class called veclist. After setstim() exits, for each node in the model there is a corresponding element in veclist that references a Vector whose values represent samples of the (precalculated) time course of extracellular potential at that node. Indeed, those elements are the ONLY references to the extracellular potential Vectors. Get rid of veclist and you reduce the reference count of those Vectors to 0, so that they are destroyed before you can say "Norbert Wiener."

For a more memory efficient way to implement extracellular stimulation, which is suitable if the extracellular medium can be treated as non-dispersive, see
http://www.neuron.yale.edu/ftp/ted/neur ... nd_rec.zip
You'll also want to read viewtopic.php?f=28&t=168
mzenker

Re: Extracellular Fields

Post by mzenker »

Hi Ted,

thank you for clarifying this!
I will reconsider the other approach, but I will probably stay with this one, because
* it is rather simple to implement,
* it works for me so far,
* my model is rather small (~950 sections at the moment - I imagine that networks of neurons have much more sections),
* I will need to use results from an externally calculated potential distribution later,
* I am not yet sure if I will have to include frequency-dependent properties of the external medium.

Matthias
Post Reply