Apply an extracellular cathodic current

[pass1945]

Hi!
So I just started to use NEURON and I'm tring to create an extracellular cathodic current point stimulation on a passive dendrite that contains 21 compartments. Is it true that the only function that I can use for this type of stimulation is e_extracellular? Are there any functions that have already been built for extracellular current stimulation? If I can only use e_extracellular, then should I manually calculate out the extracellular voltages at each compartments, depending on the various distances r and use equation V=I/(4*pi*sigma*r), and then assign the calculated voltages to e_extracellular(x)? It's fine if I do this, but can I do these calculations in NEURON only or do I have to import the data into Matlab, calculate out the corresponding voltages, then assign the values back into NEURON? It seems like this is definitely not the easiest way.
Thanks a lot!

[ted]

Is it true that the only function that I can use for this type of stimulation is e_extracellular?

There are many ways to represent extracellular stimulation. One is to use Rattay's "activating function" approach; you'll find an example of this by McIntyre et al. in ModelDB. However, extracellular will be at the core of any approach that involves a direct representation of the potential immediately adjacent to the external surface of a cell.

import the data into Matlab, calculate out the corresponding voltages, then assign the values back into NEURON?

That would work; a similar approach can be used with COMSOL if the extracellular medium is assumed to be dispersive.

Code: Select all

Are there any functions that have already been built for extracellular current stimulation?

Take a look at
Extracellular stimulation and recording
in the Hot tips area of the NEURON Forum.

[pass1945]

Hi Ted,
So I'm trying to stimulate a cell (including a soma, an apical dendrite, and a basal dendrite) extracellularly. I had an extracellular current point source that inputs in a sine wave, and then I calculated the corresponding extracellular voltages on each end points of the segments. There are 21 segments in the apical dendrite, one in soma, and 11 in the basal dendrite. Then I used .play(&e_extracellular(x),tvec) to apply the extracellular voltages (as a function of t) in the corresponding points. Theoretically, the extracellular voltages that I calculated should be a sine wave as well, since it's just current scaled by (4*pi*sigma*r). However when I plotted the e_extracellular, it gave me a really weird looking shape (it was a like a straight line but with a sudden peak somewhere around 100ms). I don't understand why.
Here's my code for first creating a sinusoidal current:

Code: Select all

//Define extracellular resistivity and position of the electrode
x= 50 //um
//minimum distance between the position of the electrode and dendrite
//This is the distance between the edge of dendrites when the somas are touching each other in an complete simulation case
rhoa=2.5 //Ohm.m

// Define I0 which is the amplitude of the oscillatory extracellular current stimulation waveform
I0 = I_manual

// Create a time vector
tstop = 1000 //ms
dt = 0.025 //ms
num = tstop/dt +1

objref tvec
tvec = new Vector(num)
for i=0, tvec.size()-1{
tvec.x[i]=dt*i
//print tvec.x[i]
}

PI = 3.1415926
f = frequency

objref I //current
I = new Vector(num)
I = tvec.mul(2*PI*f*(180/PI))
I.apply("sin")
I.mul(I0) //in nA

And here is my code for calculating the extracellular voltages from the current, and then applying them to the corresponding points on the apical dendrite:

Code: Select all

proc setStimApical(){
for i=0, napical-1{
access apical[i]
for j= 0, apicalseg{ //For every point on apical dendrite (would be apicalseg+1 points)
position1 = j*(1/apicalseg) //This is the position of the point

r1 = sqrt((x+Lsoma/2)^2+(Lsoma/2+apicalsegL*j)^2) //Distance from the extracellular stimulation point to the point

//Calculate the corresponding extracellular voltage
ve_apical[j]=new Vector()
ve_apical[j]=I.c
ve_apical[j].div(4*PI*(1/rhoa)*r1) //multiply 10^-6 converts um (from r) to m (to cancel the m in rhoa); then divide 10^-6 to convert nA to mA

print ve_apical[j].max, " ", ve_apical[j].min


//Apply the extracellular voltage onto the corresponding point on the apical dendrite
ve_apical[j].play(&e_extracellular(position1),tvec)
}
}
}

Do you see anything wrong here?
Thank you!

[ted]

Code: Select all

tvec = new Vector(num)

at which point
tvec.indgen()
tvec.mul(dt)
gives you a vector with contents 0,dt...(num-1)*dt

BTW, you don't need to define PI--NEURON has it already.

Code: Select all

I = tvec.mul(2*PI*f*(180/PI))

Is f in Hz? NEURON's t is in ms.
What's the (180/PI) for?
tvec.mul(foo) replaces tvec's contents with foo*old_values, not at all what you intended. You probably meant

Code: Select all

I = tvec.c
I.mul(2*PI*f/1000)
// or the faster, but maybe more prone to human misreading,
// I = tvec.c.mul(2*PI*f/1000)
// which first copies tvec, then multiplies this copy by the scale factor

Code: Select all

for i=0, napical-1{
access apical[i] . . . 

Each pass of this loop pushes a new section onto the section stack, making it the new "currently accessed section"--but never pops it off. Consequently, any "naked" range variable or section variable (one that is not paired with a disambiguating section name, e.g. v, diam, L, Ra appearing by itself) will mean something different on each pass. This includes instances that appear in Graphs and print statements. When transiently changing the currently accessed section, it's much better to use "section stack syntax"
sectionname statement
For this case, the code becomes

Code: Select all

for i=0,napical-1 apical[i] {
  . . . statements that pertain to apical[i] . . .
}

Code: Select all

for j= 0, apicalseg{ //For every point on apical dendrite (would be apicalseg+1 points)
position1 = j*(1/apicalseg) //This is the position of the point
 . . . etc. . . .

You're trying to drive e_extracellular at the boundaries between compartments. That variable is not available to you. NEURON uses the central difference approximation for the 2nd spatial derivative. Consequently it solves ODEs for membrane potential at normalized distances 0.5/nseg . . . (nseg - 0.5)/nseg from the 0 end of each section. Those are the locations of a section's internal nodes, and those are where you want to drive e_extracellular.
The idiom is

Code: Select all

for (x,0) {
  statements that pertain to the internal node at location x
}

"What harm could come from trying to drive e_extracellular at 0 or 1?"

Attempts to reference a range variable at 0 or 1 end up affecting that range variable at the nearest internal node. You'd end up with two Vectors competing to drive e_extracellular at 0.5/nseg, and another two competing to drive e_extracellular at (nseg - 0.5)/nseg. The result would depend on execution sequence, and would not necessarily reflect your intent.

[pass1945]

Thank you so much Ted for the quick reply!
I fixed I and tvec, and then I used "for (x,0)" for each section. The results that came out were very similar to what I did with calculating each position, but I'm sure that with this "for (x,0)", the method was more accurate. And indeed it's much easier.
Again, thank you Ted!

[ted]

You're welcome. Did the simulation produce reasonable results or do they still look weird?