Extracellular stimulation of cutaneus afferents

[ksfrahm]

Hello everybody,

I am trying to develop a model of small cutaneous afferent nerves. These nerves are characterized by being myelinated inside the dermis but having an unmyelinated part inside the epidermis.
My present model consists of 6 sections, 1 unmyelinated part, 3 nodes and 2 axons with myelin. The unmyelinated and node part is modeled using HH dynamics and the rest as passive. I have tested my model using intracellular stim, and that seems to work. Now I am trying to apply an extracellular field which is the purpose of my model. However, I have got some problems understanding the extracellular mechanism. The external field is calculated using a finite element model. I have gone through the example and posts especially extracellular_stim_and_rec.zip, ic_stim.zip and relevant post like http://www.neuron.yale.edu/phpBB/viewto ... =15&t=1937

My code:

Code: Select all

load_file("nrngui.hoc")

v_init = -70 // mV
mycm = 0.1   // uF/cm2/lamella membrane//
mygm = 0.001 // S/cm2/lamella membrane
nl = 10      // myelin layers
create axon[3], node[2], unmyelinated[1]

forsec "axon" {
	nseg = 3
	L = 100   	// um
	diam = 2	 // um
	Ra = 70		// Ohm/cm
	cm = 2 		// membrane capacitance
	insert pas
	g_pas = 0.0303
	e_pas = v_init
	insert extracellular xraxial=277965 xg=mygm/(nl*2) xc=mycm/(nl*2) e_extracellular = 0
}
forsec "node"{
	nseg = 3
	L = 1
	diam = 1
	Ra = 70
	cm=2
	insert hh 
	insert extracellular xraxial=557823 xg=1e10 xc=1e9 e_extracellular = 0
	} 
forsec "unmyelinated"{
	nseg = 3
	L = 20 // um
	diam = 1
	Ra = 70
	cm=2
	insert hh
	insert extracellular xraxial=557823 xg=1e10 xc=0 e_extracellular = 0
	}
	// Connect pieces
	connect unmyelinated[0](1), axon[0](0)
	for i = 0, 1 {
		connect axon[i](1), node[i](0)
		connect node[i](1), axon[i+1](0)
	}

	finitialize(v_init)
	fcurrent() 
	
load_file("ex_stim2.hoc")
PMAX = 0 // mV
dur = 0.5
del = 1
stim_vec(PMAX, dur, del)
load_file("RC_and_v_graph.ses")
load_file("space_plot.ses")
run()

My ex_stim2.hoc file:

Code: Select all

objref tvec, pvec, veclist
objref extern_field
extern_field = new Matrix(3,6) // values define by Finite element model
extern_field.x[0][0] = -91.7
extern_field.x[1][0] = -91.7
extern_field.x[2][0] = -91.7

extern_field.x[0][1] = -91.6
extern_field.x[1][1] = -91.6
extern_field.x[2][1] = -91.5

extern_field.x[0][2] = -91.5
extern_field.x[1][2] = -91.5
extern_field.x[2][2] = -91.5

extern_field.x[0][3] = -91.5
extern_field.x[1][3] = -91.5
extern_field.x[2][3] = -91.4

extern_field.x[0][4] = -91.4
extern_field.x[1][4] = -91.4
extern_field.x[2][4] = -91.4

extern_field.x[0][5] = -91.4
extern_field.x[1][5] = -91.3
extern_field.x[1][5] = -91.3

objref sl
sl = new SectionList()
proc setstim() { localobj tmpvec
for sec_i = 0, 5 {
  veclist = new List()
  forsec sl(sec_i){
  i = 0
  for (x, 0) {  // iterate over internal nodes only, here 3
      tmpvec = pvec.c
	  tmpvec.x[i] = extern_field.x[i][sec_i]
	  i = i+1
      tmpvec.play(&e_extracellular(x), tvec)
      veclist.append(tmpvec)
    }
   }
  }
}

proc stim_vec() {
// 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 = $1 
dur = $2
del = $3
{ tvec.x[0]=0    pvec.x[0]=0 }  // field is 0 for 1 ms 
{ tvec.x[1]=del   pvec.x[1]=0 }
{ tvec.x[2]=del    pvec.x[2]=1 }  // jumps to some value
{ tvec.x[3]=del+dur  pvec.x[3]=1 }  //   for 1 ms
{ tvec.x[4]=del+dur  pvec.x[4]=0 }  // falls back to 0 ever after
setstim()
}

When I run the model with the extracellular field I do not see an AP whatever stimulation parameters I set, how can that be? The extracellular field does not seem to change and the membrane potential does not change which puzzles me a bit. What am I doing wrong?

Would there be any other approaches which could more elegant and simple than a vector for each section?
My plan is to use an approach similar to the MRG model (Modeling the Excitability of Mammalian Nerve Fibers: Influence of Afterpotentials on the Recovery Cycle, McIntyre et al., 2002) to model the myelinated part of the axon. This means I will have even more sections in my model and if I apply branching axon, the extracellular mechanism could be very hard to define.

Any help will be greatly appreciated!

Steffen

[ted]

I have tested my model using intracellular stim, and that seems to work.

First it would be good to verify that your computational model's anatomical and biophysical properties are identical to those of your conceptual model. For example, do you really want your "axon" sections to have g_pas = 0.0303 S/cm2? That's enormously leaky--these sections will have a 66 microsecond time constant. What about other parameters? Is the topology of the model correct, i.e. are the sections connected to each other in the way you want? You might try
topology()
forall psection()
and use the GUI's Model View tool to explore and verify the model specification.

I have gone through the example and posts especially extracellular_stim_and_rec.zip, ic_stim.zip

What's ic_stim.zip and where did you find it?

When I run the model with the extracellular field I do not see an AP whatever stimulation parameters I set, how can that be? The extracellular field does not seem to change and the membrane potential does not change which puzzles me a bit.

Have you verified that your stimulus vectors exist and that their contents are what you expect them to be? Have you verified that they play their values into the corresponding instances of e_extracellular in the course of a simulation?

Strong suggestion:
Take a big step way back from all this complexity. Make a simple model with just one section and apply extracellular stimulation to it. After that works properly, make a model with two sections and apply extracellular stimulation to that. Get that to work, and you'll be ready to move on to more complex models. If your model is too complex to debug, and you have too many untested bits of code, you're lost.

Do you really want to throw away veclist and build a new one on every pass through this for loop?

Code: Select all

for sec_i = 0, 5 {
  veclist = new List()
  . . .

Because that's what this code will do. Or is your real aim to build just one list that holds all the Vectors that will be used to stimulate your model? in which case the
veclist = new List()
statement should be executed _before_ entering the for loop.

This
forsec sl(sec_i)
is not how to iterate over the contents of a SectionList. Have you read the Programmer's Reference documentation of the SectionList class? Have you created a SectionList, appended a few sections to it, then iterated over the list to verify that you know how to do it? A useful test is to print out the name of each section that has been appended to a SectionList.

Would there be any other approaches which could more elegant and simple than a vector for each section?

Yes. If the potential at each point j is the product of a stimulus current i(t) and a transfer resistance r(j), you could use the approach illustrated in extracellular_stim_and_rec.zip

if I apply branching axon, the extracellular mechanism could be very hard to define.

What would make this more difficult than dealing with an unbranched axon? Conceptually it's the same task. The beauty of digital computers is ease with which they handle repetitious tasks.

[Qroid_montreal]

Have you verified that they play their values into the corresponding instances of e_extracellular in the course of a simulation?

What are the general ways to do this? If I'm playing a vector to e_extracellular, how would I access e_extracellular(t) for each section and each time point?

[ted]

e_extracellular is a range variable. You can access the value of any range variable in any segment that has the variable. Whether you print or plot it is up to you. For this particular variable, you might find a space plot (plot of variable vs. distance along a path, which evolves in the course of a simulation) particularly helpful.

[Qroid_montreal]

Say I want to store e_extracellular(0.5) for each section at each time point in another variable for later analysis, how would I do this? I'm confused how to access and store variables (range or otherwise) as the simulation runs. As always, many thanks.

[ksfrahm]

Hi ted

Thanks for the great replies to my post. I am trying to do all the steps you mentioned.

e_extracellular is a range variable. You can access the value of any range variable in any segment that has the variable. Whether you print or plot it is up to you. For this particular variable, you might find a space plot (plot of variable vs. distance along a path, which evolves in the course of a simulation) particularly helpful.

To Qroid_montreal's post you mention to use the space plot, I actually have been trying to use this based on the extracellular_stim_and_rec.zip example, but I am not sure I am doing it correctly. I am also plotting the v and v + vext. How do I define where 0 is placed on the x-axis is?
my code is

Code: Select all

objectvar save_window_, rvp_
objectvar scene_vector_[1]
origin_val = 25
{
save_window_ = new Graph(0)
save_window_.size(-376.823,328.433,-80,40)
scene_vector_[0] = save_window_
{save_window_.view(-300, -200, 350, 250, 300, 370, 300.48, 200.32)}
flush_list.append(save_window_)
save_window_.save_name("flush_list.")
objectvar rvp_
rvp_ = new RangeVarPlot("v")
unmyelinated[0] rvp_.begin(1)
axon[2] rvp_.end(1)
rvp_.origin(origin_val)
save_window_.addobject(rvp_, 2, 9, 0.640256, 0.904792)

objectvar rvp_
rvp_ = new RangeVarPlot("v($1)+vext($1)")
unmyelinated[0] rvp_.begin(1)
axon[2] rvp_.end(1)
rvp_.origin(origin_val)
save_window_.addobject(rvp_, 3, 7, 0.64984, 0.880831)

objectvar rvp_
rvp_ = new RangeVarPlot("e_extracellular")
unmyelinated[0] rvp_.begin(1)
axon[2] rvp_.end(1)
rvp_.origin(origin_val)
save_window_.addobject(rvp_, 4, 7, 0.64984, 0.880831)
}

What's ic_stim.zip and where did you find it?

I apologize I meant the ecstim.zip

First it would be good to verify that your computational model's anatomical and biophysical properties are identical to those of your conceptual model. For example, do you really want your "axon" sections to have g_pas = 0.0303 S/cm2? That's enormously leaky--these sections will have a 66 microsecond time constant. What about other parameters? Is the topology of the model correct, i.e. are the sections connected to each other in the way you want? You might try
topology()
forall psection()
and use the GUI's Model View tool to explore and verify the model specification.

You are correct my g_pas was an old value I found in another example. I have changed it to 0.0001 S/cm2, which the standard value used for similar approaches in mammalian fibers. I used the topology() and forall psection() and the connections seem to be OK, however, axon[2](0) was connected to node[1](1) instead of node[1](0), so after re-arranging a bit I got

Code: Select all

|---|       axon[2](0-1)
       `--|       node[1](0-1)
            `--|       axon[1](0-1)
                 `--|       node[0](0-1)
                      `--|       axon[0](0-1)
                           `--|       unmyelinated[0](0-1)

Instead of

Code: Select all

|---|       axon[2](0-1)
 `--|       node[1](1-0)
     `--|       axon[1](1-0)
         `--|       node[0](1-0)
             `--|       axon[0](1-0)
                 `--|       unmyelinated[0](1-0)

I have changed my stimulation code into the following based on your suggestions.

Code: Select all

objref sl
sl = new SectionList()
sl.wholetree()
veclist = new List()
proc setstim() { localobj tmpvec
  sec_i = 0
  forsec sl{
  //print secname()
  i = 0
  for (x, 0) {  // iterate over internal nodes only, here 3
      tmpvec = pvec.c
	  tmpvec.x[2] = extern_field.x[i][sec_i]
	  tmpvec.x[3] = extern_field.x[i][sec_i]
	  tmpvec.play(&e_extracellular(x), tvec)
	  veclist.append(tmpvec)
	  i = i+1
        }
	sec_i = sec_i+1
	}
  }

After changing the code above, my space plot does show that the extracellular potential changes, to approx - 91 mV as it should, but the 'v' values does not indicate the initiation of an AP. Nor does the voltage graph .

Thanks again for any further help
Steffen

[ted]

How do I define where 0 is placed on the x-axis is?

The x coordinate of any point in a space plot is the path distance from the origin of the model cell, i.e. from the 0 node of the root section (the section that has no parent).

I have changed my stimulation code

In order to provide further assistance at this point, I'll have to be able to reproduce and observe your simulation results myself. If you zip up just the hoc, ses, mod, and whatever data files (time course of extracellular potentials?) are necessary and email them to me
ted dot carnevale at yale dot edu
I'll tell you what I discover.

[ted]

Summary of email to ksfrahm: The simulation code appears to be working properly. Membrane potential at the nodes depolarizes, then hyperpolarizes, and the sodium conductance's m gating variable responds properly. A larger stimulus may help, or it may be possible to get an anode break spike by changing the polarity of the stimulus.