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.
Sherif
Posts: 27
Joined: Thu Oct 13, 2005 4:32 pm

Extracellular Fields

Post by Sherif » Thu Oct 13, 2005 5:14 pm

Hello,

Does anybody use NEURON to simulate the effect of extracellular field potentials (e.g. homogenous and non-homogenous fields) on the behavior of neurons in anatomically detailed models?

Is there any study that is dealing with a similar problem?

I am starting a new project regarding this problem and it would be great to know different ways of using Neuron to simulate the application of different extracellular fields.

Best Regards,

costas

fluctuating (time-, space-dependent) extracellular fields

Post by costas » Tue Nov 06, 2007 2:57 pm

Hello

I am trying to apply/interface a time- and space-dependent extracellular field with a neuronal structure. For the moment, the properties of the cable are purely passive.

I am trying to define an extracellular field through Biophysics/extracellular but I haven't been successful yet in defining it as a time-, space-dependent function. Could anybody help me with this? Along the same lines, how do I interface the extracellular field calculated with another software package with NEURON?

Thanks a lot, Costas

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

Post by ted » Wed Nov 07, 2007 10:47 am

See https://www.neuron.yale.edu/phpBB2/viewtopic.php?t=168
how do I interface the extracellular field calculated with another software package with NEURON?
The most straightforward way is to use the other software to write one or more files that
contain the computed field results, and read them into Vectors with NEURON.

costas

e_extracellular

Post by costas » Wed Nov 07, 2007 9:08 pm

Dear Ted

Thanks a lot for the advice! I managed to build a simple model using the code by McIntyre.

If I am not wrong, I should be able to apply the extracellular voltage gradient parallel to the cable by setting the magnitude of the extracellular voltage through e_extracellular. When I tried to implement this in the code it didn't seem to see this variable at all, i.e. e_extracellular. Changing the extracellular properties (resistance, capacitance) didn't alter the outcome either...

For suggestions I would be thankful.

All the best, Costas

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

Post by ted » Thu Nov 08, 2007 10:56 am

My mistake. McIntyre et al. use the "activating function" approach, but the ModelDB entry
does not include their code for implementing the activating function current. If you ask them,
perhaps they will provide it.

If instead you want to drive e_extracellular, you have a choice:
1. The most general approach, which works even if the extracellular medium is nonlinear
or has reactance (capacitive impedance) that cannot be ignored. In pseudocode,

Code: Select all

for each nonzero area node {
  initialize a Vector that contains the time course of the
    extracellular field at that point
  use the Vector class's play() method to make that Vector
    drive e_extracellular at that point
}
Here's a skeleton for this in hoc:

Code: Select all

objref tvec // for the stim waveform's sample times
objref tvec = new Vector()
// here put some code that fills tvec with the sample times

objref veclist  // will hold all the stim Vectors

proc setstim() { localobj tmpvec
  veclist = new List()
  forall {
    for (x, 0) {  // iterate over internal nodes only
      // here put some code that loads tmpvec with the time course
      //   of extracellular potential at this particular location
      // might generate algorithmically,
      //   or read from a file of precalculated values
      tmpvec.play(&e_extracellular(x), tvec)
      veclist.append(tmpvec)
    }
  }
}
2. A potentially more efficient approach that can be used if the extracellular potential at
all points is separable into two components, one time-dependent and the other location-
dependent, i.e. Vextracellular(t, x, y, z) = rx(x, y, z) f(t)
Requires the extracellular medium to be linear and nonreactive--not much of a restriction,
since most papers make this assumption.
In pseudocode:

Code: Select all

initialize a single Vector that contains f(t)
use the Vector class's play() method to make that Vector
  drive a global variable that belongs to a density mechanism
code in the density mechanism scales f(t) by rx (whose value was assigned
  during model setup),
  and assigns this product to a POINTER variable
a hoc-level setpointer statement, executed once during model setup, 
  links the POINTER variable to local e_extracelluar
This is the method implemented by the code in the zip file I referred to earlier.

costas

Post by costas » Sat Nov 10, 2007 8:35 pm

Dear Ted

Thanks for the advice, I decided to follow the 1st option to apply an extracellular voltage on a 1-segment dendrite. I am attaching the code below. Although the process sees the external file ("time.vect.txt") and uploads it, when I then print e_extracellular as a function of time it constant and equal to zero.

Code: Select all

//physical properties of the cell
create a
{access a nseg=1 L=10 diam=3
	insert pas g_pas=1 cm=.5e3*g_pas e_pas=-70
}
{insert extracellular xraxial=1e9 xg=1 xc=0}

tstop=5

// graphical interface appearance
objref g
g = new Graph()
g.size(0,5,-80,-60)
/*g.addexpr("vext(.1)")
g.addexpr("vext(.3)")*/
g.addexpr("vext(.5)")
/*g.addexpr("vext(.7)")
g.addexpr("vext(.9)")*/

// simulation control
objref m, datafile, tvec, tmpvec, veclist

proc run() {
	// Find the data file with the time-/voltage-profile
	datafile = new File()
	datafile.ropen("/Users/costas/Neuron_codes/time_vect.txt")

	// Define all the vectors that define the extracellular voltage
	m = new Matrix()
	tvec = new Vector()
	tmpvec = new Vector()
	veclist = new List()

	// Read the data file with the extracellular voltage
	// m.scanf(datafile, n_rows, n_columns)
	m.scanf(datafile,5,2)
	// The time vector is defined as the first column in the data file.
	// The numbering starts with 0.
	tvec=m.getcol(0)
	// The voltage vector for all internal nodes of each segment (the latter here is 1).
	// Note that the size of this vector should be = nseg + origin + end + 1
	// The last + 1 is due to the number of nseg that starts with 0 and not with 1
	tmpvec=m.getcol(1)

	// Define the value of the membrane voltage at the origin
	v = -60
	// e_extracellular = 1 

	// Load the temporal evolution of the extracellular voltage profile for each internal node
	// The location of the internal node is given by number x (x=0: start; x=1: end of the cable) 
	for (x, 0) { 
      	tmpvec.play(&e_extracellular(x), tvec)}
      	veclist.append(tmpvec) 
	
	t=0
	g.begin()
	g.beginline()
	
	while (t < tstop) {
		g.line(t, v(.5) + e_extracellular(.5))
		g.plot(t)
		fadvance()
	}
	g.flush()
}
// run simulation
run()

I am not sure what the problem is, any help would be greatly appreciated.
[/code]

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

Post by ted » Sun Nov 11, 2007 1:56 am

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:

Code: Select all

load_file("stim.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:

Code: Select all

load_file("plot_e_extra.ses")
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

costas

Post by costas » Tue Nov 13, 2007 5:59 pm

Dear Ted

Thank you so much for the detailed reply. I have now implemented the code in order to simulate the effect of a harmonic extracellular field on a single unbranched dendrite.

I had a question regarding this: for an unbranched cable with two von Neumann-type boundary conditions it is exactly the case that you sent me. In order to simulate a Dirichlet-type condition on the one boundary I have added a soma with e_pas=V_rest. When I simulate this simple case the numerical results from NEURON are different from the analytical expressions that I have derived (the latter is also consistent with the numerical solution of the cable equation for a harmonic extracellular field obtained using other software packages.) This is somewhat awkward and I am puzzled about the cause...

What do you think are the most prominent bug-candidates for this? Once more, I thank you for you time and effort.

All the best, Costas

costas

Post by costas » Tue Nov 13, 2007 8:44 pm

Sorry for not mentioning this in my last entry:

It seems that for the case where I apply a Dirichlet-type boundary condition (BC) on one side and a von Neumann-type BC on the other, the membrane voltage at the first internal node from the Dirichlet-type BC becomes discontinuous. This lead me to the question regarding the numerical solvers. On the other hand, it could be that I have set a parameter erroneously.

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

Post by ted » Tue Nov 13, 2007 10:26 pm

von Neumann was born too late; the credit belongs to old Carl, who wasn't a von.

To answer your latest questions I must first ask you to please be more specific about what
you mean by Dirichlet and Neumann boundary conditions in the present context, and in
particular how you tried imposing them.

costas

Post by costas » Tue Nov 13, 2007 11:52 pm

Dear Ted

Thanks for the clarification. :) Well, Dirichlet is the killed end boundary condition (V_m(x=x0)=V_0) and Neumann is the sealed end (dV_m/dx(x=x0)=0) boundary condition.

When I try to apply a killed end and a sealed end boundary condition at the two ends of a single unbranched dendrite in the presence of an extracellular field. In order to do so I simply attach the dendrite on one end to a soma. Because I assume that the membrane potential at the soma remains unaffected from the extracellular field I attribute only passive properties to it (L, diam, e_pas) as well as a capacitance cm while not applying any field. Therefore, on one side the dendrite has a killed end BC with V_m(x=0)=V_rest while on the other end dV_m(x=L)/dx=0. As for the rest, the model is based on the code you sent to me on Sunday which applies for a single unbranched dendrite with two sealed ends and works fine.

Now, when I run the killed end/sealed end case the membrane voltage at the first internal node from the killed end becomes discontinuous in the sense the it starts (correctly) at x=0 from V_rest, then in the first node jumps to some value and then in the rest of the internal nodes it seems like it is trying to relax between the that value and the sealed end boundary condition at the other end. The value in this first internal node seems to be very discontinuous compared to the rest of the internal nodes. This might be a sign that something is wrong with my process-definition.

Here is the code for the soma-dendrite definition:

Code: Select all

create soma, axon
connect soma(0), axon(0)

soma {
nseg = 1
L = 50
diam = 50}

axon {
nseg = 21
L = 1000
diam = 2
}

forall {
	Ra=100
	cm=1
	insert pas
	g_pas=0.0002
	e_pas=-65
}

access axon
and here is the code that applies the extracellular field on the dendrite only:

Code: Select all

objref tvec, pvec[nseg], datafile, dummy, m
// REMEMBER: nseg has to be an ODD number

access axon

// tvec will hold the stimulus sample times
// pvec will be a 1 ms pulse of with amplitude PMAX

NUMPTS = 10
tvec = new Vector(NUMPTS)
pvec = new Vector(NUMPTS)

// Find the data file with the time-/voltage-profile
datafile = new File()
//datafile.ropen("/Users/costas/Neuron_codes/time_vect.txt")
datafile.ropen("dummy.txt")

// Read the data file with the extracellular voltage
// m.scanf(datafile, n_rows, n_columns)
m = new Matrix()
m.scanf(datafile,NUMPTS,nseg+1)

// The time vector is defined as the first column in the data file.
// The numbering starts with 0.
tvec=m.getcol(0)

// The voltage vector for all internal nodes of each segment (the latter here is 1).
// Note that the size of this vector should be = nseg + origin + end + 1
// The last + 1 is due to the number of nseg that starts with 0 and not with 1
for ii=0, nseg-1 {
	pvec[ii]=m.getcol(ii+1)}


// 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

	if (x < 0.025) {ii=0}
	if (x > 0.025 && x <= 0.080) {ii=1}
	if (x > 0.080 && x <= 0.130) {ii=2}
	if (x > 0.130 && x <= 0.180) {ii=3}
	if (x > 0.180 && x <= 0.220) {ii=4}
	if (x > 0.220 && x <= 0.280) {ii=5}
	if (x > 0.280 && x <= 0.320) {ii=6}
	if (x > 0.320 && x <= 0.370) {ii=7}
	if (x > 0.370 && x <= 0.430) {ii=8}
	if (x > 0.430 && x <= 0.470) {ii=9}
	if (x > 0.470 && x <= 0.520) {ii=10}
	if (x > 0.520 && x <= 0.560) {ii=11}
	if (x > 0.560 && x <= 0.620) {ii=12}
	if (x > 0.620 && x <= 0.660) {ii=13}
	if (x > 0.660 && x <= 0.700) {ii=14}
	if (x > 0.700 && x <= 0.750) {ii=15}
	if (x > 0.750 && x <= 0.800) {ii=16}
	if (x > 0.800 && x <= 0.850) {ii=17}
	if (x > 0.850 && x <= 0.900) {ii=18}
	if (x > 0.900 && x <= 0.950) {ii=19}
	if (x > 0.950) {ii=20}

	tmpvec = pvec[ii]

        tmpvec.play(&e_extracellular(x), tvec)
        veclist.append(tmpvec)
    }
  }
}

setstim()
Hopefully this makes the problem somewhat clearer. Maybe I haven't set a parameter correctly but I am not sure where the problem is.

Thanks a lot for your time and effort.

All the best, Costas

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

Post by ted » Thu Nov 15, 2007 1:25 am

Excellent questions, Costas. Also very clean code.

A minor programming suggestion: make more use of Lists and iterator statements.
costas wrote:Dirichlet is the killed end boundary condition (V_m(x=x0)=V_0) and Neumann is the sealed end (dV_m/dx(x=x0)=0) boundary condition.
OK. The latter you get automatically, without having to do anything.

Assuming that nseg is at least 3, the former you can get merely by setting the equilibrium
potential of a linear ion channel (pas mechanism would do) at that location to 0, and
increasing its conductance density to a very large number, e.g. 1e6 siemens/cm2 would
probably be far more than enough. This would make the potential at the last internal node
of the "crushed" end compartment identical to the potential of the extracellular bath at that
point. This is exactly what you want, because that node is associated with membrane.
extracellular must be present, unless you are assuming that the killed end is lying in a
grounded pool of Ringer solution that is isolated from the main bath, as by a vaseline seal.

Attaching a "soma" to the 0 or 1 end might at first seem like a way to emulate a crushed
end, but don't forget to make conductance density in the soma extremely large. Also,
extracellular must be attached to the soma, and soma e_extracellular(0.5) must be driven
by the same voltage as drives the adjacent internal node of axon. Otherwise there will
be a huge capacitive transient if e_extracelluar along the axon changes abruptly--
remember that axon's internal nodes are capacitively coupled to axon's e_extracellular,
so they must follow e_extracelluar at short times, but the soma's membrane capacitance,
one side of which is grounded, prevents its internal node from following such quick
changes.

However, even if you fix that problem, anoher remains which is best seen by taking
advantage of the principle of superposition. By which I mean, set e_pas to 0 throughout
the model so that it is more easy to observe divergence of internal potentials from rest.
Then you will note that, even though the potential at the soma node remains 0 during the
field stimulus, the proximal portion of axon--starting with its first internal node, and
continuing for some distance, diverges from 0 in the opposite direction from which the
sealed end diverges. A killed end boundary condition would require that potential at the
first internal node of axon be 0.

The problem arises not because soma is an inefficient "ground" but because soma is
separated from the first internal node of axon by half the axial resistance of one
segment of axon. Axial current flowing through this resistance is sufficient to allow the
potential at the 1st internal node of axon to escape from 0. This can be overcome by
simulating axon with a much finer spatial grid.


http://www.neuron.yale.edu/ftp/ted/neur ... ed_end.zip contains files that
demonstrate all this. Unzip it, then use NEURON to load init_costas.hoc
Click on the Movie Run panel's Init & Run.
Then at the oc> prompt execute
soma g_pas*=1e6
and run another simulation. Note that soma v now stays at 0, but the proximal end of axon
diverges from 0.

Exit NEURON, then edit costas_cell.hoc to change
nseg = 21
to
nseg = 21*5
and save that file.
Then start NEURON again and load init_costas.hoc
Now the axial voltage drop between soma and the 0 end of axon is almost obliterated,
but soma and the proximal end of axon still diverge from 0 . . . because soma g_pas is
too small. So increase soma g_pas by 1e6 and now you see that the simulation is much
closer to killed end boundary condition.

Of course, a very similar result could have been obtained without using any soma
section at all, and with a much smaller number of segments in the entire model, just by
making g_pas in the first segment of axon very large--
axon g_pas(0.01)*=1e6

costas

Post by costas » Fri Nov 16, 2007 6:30 pm

Dear Ted

Thank you so much for the detailed answer and the codes. I implemented your suggestions, i.e., I defined extracellular over the soma too and made the spatial discretization much finer, both for the soma as well as for the axon. The solutions now are really smooth both for sealed end and killed end boundary conditions.

There is still a problem though regarding the resulting membrane voltage distribution. That is, I have derived the analytical solution for the induced membrane potential in the presence of a harmonic extracellular field for various boundary conditions (stationary case). The results from the numerical simulations (after all the transients die out) are different from these analytical expressions and I am not sure why.

Would it be possible to Email these analytical solutions, their derivation and the results from the NEURON code for you to take a look at. I am sorry to have to ask you (I know that I am not the only one with questions about NEURON) but I really cannot find the mistake....

Please let me know if this is OK with you.

All the best, Costas

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

Post by ted » Sat Nov 17, 2007 9:02 am

Sure, go ahead
ted dot carnevale at yale dot edu

costas

Post by costas » Sun Nov 18, 2007 2:12 am

Dear Ted

Thanks for the reply. I think I found the problem... I guess I was applying an erroneous sealed end boundary condition. I'll recalculate the whole thing to see if it now reproduces the numerical results (I'd be surprised if it didn't). I'll keep you posted. Until then, thanks so much for your time & effort.

All the best, Costas

Post Reply