Modeling dynamic extracellular fields

[Qroid_montreal]

Hi,

This is my second project with NEURON, and my first time really generating code. I'd like to model the effects of a temporally and spatially dynamic extracellular voltage on transmembrane currents in a spatially distributed neuron. Membrane voltage "v" is internal potential minus extracellular voltage "vext". So for each section I want to read in a time-vector of vext changes ("Dv(t)"), such that at each timestep i in NEURON the ith element of Dv, dv=Dv(i) is subtracted from the v value.

My plan is to create a modified Hodgkin-Huxley mechanism ("hh2") whose only modification is an additional term in the dv/dt equation, such that

Code: Select all

dv/dt = (-Ina - Ik - Il + Istim)*1/cm - dv

Will this do what I'm expecting ie. model the sections as if situated in an external field? If so, how do I edit the hh2.mod to do it? I can't find hh.mod in my NEURON files, but from what I can tell, dv/dt isn't in there. So where would I change the voltage ODE?

I'll use ropen to read in the each DV element from an ASCII file. This is correct, right?

Many thanks.

[sgratiy]

I thought that simply inserting the existing extracellular mechanism should take care of an additional current caused by changing vext without a need to modify the hh mechanism, but I may be incorrect.

[ted]

I'd like to model the effects of a temporally and spatially dynamic extracellular voltage on transmembrane currents in a spatially distributed neuron.
. . .
My plan is to create a modified Hodgkin-Huxley mechanism ("hh2") whose only modification is an additional term in the dv/dt equation, such that

Code: Select all

dv/dt = (-Ina - Ik - Il + Istim)*1/cm - dv

Will this do what I'm expecting

Nope. And there's no reason to change any ODEs or tinker with any mod files. All you have to do is insert the extracellular mechanism into the affected sections, then force e_extracellular at each internal node to follow the desired time course. You can do this with the Vector class's play method, which is documented in the Programmer's Reference and probably illustrated in various posts here on the Forum.

[Qroid_montreal]

That's great to hear. I'll play with this tomorrow.

Just so I'm clear, what would happen if I implemented the modified hh equation as written? If I don't see the error, am I missing something on the math/conceptual or the NEURON level?

edit: Thank you!

[ted]

Just so I'm clear, what would happen if I implemented the modified hh equation as written?

Try this thought experiment: suppose you block all ionic currents in a real cell. What will happen to membrane potential? How does this differ from what your modified equation would predict?

If I don't see the error, am I missing something on the math/conceptual or the NEURON level?

Both, but let's start with "math/conceptual."

[Qroid_montreal]

Try this thought experiment: suppose you block all ionic currents in a real cell. What will happen to membrane potential? How does this differ from what your modified equation would predict?

Hmm. Excuse my density... From the original Hodgkin-Huxley equations, the membrane potential can't change without ionic currents (the RHS of dv/dt is all zero). But Vm=Vi-Ve, so what's going on in the situation when Ve is changing? Is the charge buildup on the inside wall of the membrane changing with Ve, ie. is Vi changing to negate any effect to Vm? Or does dv/dt=sum(currents)/(-cm) fail to capture it? From your question I'm guessing it's the former, but clearly I'm at a loss.

And I realise this is a forum for NEURON not biophysics help, so thank you.

[ted]

Try this thought experiment: suppose you block all ionic currents in a real cell. What will happen to membrane potential? How does this differ from what your modified equation would predict?

Hmm. Excuse my density... From the original Hodgkin-Huxley equations, the membrane potential can't change without ionic currents (the RHS of dv/dt is all zero).

So what does your equation for dv/dt predict?

[Qroid_montreal]

My equation predicts Vm can change without ionic currents.

But why is this wrong? If Vm=Vi-Ve, doesn't Vm change whenever the extracellular potential changes, regardless of currents?

[ted]

My equation predicts Vm can change without ionic currents.

But why is this wrong? If Vm=Vi-Ve, doesn't Vm change whenever the extracellular potential changes, regardless of currents?

Your equation violates fundamental physical principles.

Example: consider a spherical cell that has only membrane capacitance and "leak" channels with resting potential 0 mV. An electrical equivalent circuit model of this cell is

Code: Select all

    o Vi
    |
+---+---+
|       |
Rm      Cm
|       |
+---+---+
    |
    o Vo

Its membrane potential is Vi-Vo.

At the start of an experiment, the cell is sitting in a bath of Ringer and the bath is grounded so Vo is 0, and it has been sitting there for a long time so Vi is also 0.

Code: Select all

    o Vi
    |
+---+---+
|       |
Rm      Cm
|       |
+---+---+
    |
    o Vo
    |
   --- gnd
    -

Its membrane potential is Vi-Vo = 0.

At time t1, you insert a 3 volt battery between the Ringer bath and ground.

Code: Select all

    o Vi
    |
+---+---+
|       |
Rm      Cm
|       |
+---+---+
    |
    o Vo
    |
  ----- +3V.
   ---
    |
    o
    |
   --- gnd
    -

Result? Vo and Vi both jump up to +3V "instantaneously" (well, it takes a small fraction of a microsecond for this to happen, but on time scales that are relevant to neurons the change is instantaneous). And Vm is still 0 mV.

This stuff isn't intuitive, but the principles aren't that hard to learn. A one semester course on electrostatics would be very helpful.

Back to your original question: all you have to do to simulate the effect of an extracellular field on a NEURON model is to insert the extracellular mechanism into the model's sections, then drive each segment's e_extracellular with a signal that corresponds to the time course of extracellular potential at that point in space. This has been done, and the source code for doing it again is available for attributed re-use--see
Extracellular stimulation and recording
viewtopic.php?f=28&t=168

[Qroid_montreal]

In your example, is Vm=0 precisely because there's no current flowing across the membrane? So there's no voltage drop because V=IR=0 - is that right?

Thanks a lot. The approach in the link looks very explicit. I'm working on that now.

[Qroid_montreal]

Back to your original question: all you have to do to simulate the effect of an extracellular field on a NEURON model is to insert the extracellular mechanism into the model's sections, then drive each segment's e_extracellular with a signal that corresponds to the time course of extracellular potential at that point in space. This has been done, and the source code for doing it again is available for attributed re-use--see
Extracellular stimulation and recording
viewtopic.php?f=28&t=168

Just so I'm clear, in the example where the cell is stimulated by an extracellular electrode (initxstim.hoc), are the stimulating and recording electrode selfsame ie. at the same x,y,z location?

Very quickly, what would I need to do stimulate and record at different coordinates? Would it amount to much more than defining and calculating two transfer impedances (eg. rxstim_k and rxrec_k) for each compartment k? Could this potentially be done inside calcrx.hoc? This is more to get a handle on the code than anything else. Thanks again.

[ted]

In your example, is Vm=0 precisely because there's no current flowing across the membrane?

It starts out == 0 as its initial condition. It stays == 0 because there is no transmembrane current to charge the membrane capacitance.

[ted]

Just so I'm clear, in the example where the cell is stimulated by an extracellular electrode (initxstim.hoc), are the stimulating and recording electrode selfsame ie. at the same x,y,z location?

Yes. That's easiest to implement because it requires computing only one set of transfer resistances.

what would I need to do stimulate and record at different coordinates? Would it amount to much more than defining and calculating two transfer impedances (eg. rxstim_k and rxrec_k) for each compartment k?

Precisely.

Could this potentially be done inside calcrx.hoc?

I would suggest revising both xtra.mod and related hoc procs and funcs. Specifically

change
RANGE rx, er
to
: RANGE rx, er
RANGE rxs, rxr, er

In the INITIAL block change

Code: Select all

        ex = is*rx*(1e6)
        er = (10)*rx*im*area

to

Code: Select all

:        ex = is*rx*(1e6)
:        er = (10)*rx*im*area
        ex = is*rxs*(1e6)
        er = (10)*rxr*im*area

Also change these

Code: Select all

BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
  ex = is*rx*(1e6)
}
AFTER SOLVE { : after each solution step
  er = (10)*rx*im*area
}

to

Code: Select all

BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
:  ex = is*rx*(1e6)
  ex = is*rxs*(1e6)
}
AFTER SOLVE { : after each solution step
:  er = (10)*rx*im*area
  er = (10)*rxr*im*area
}

Then would use hoc statements to calculate the appropriate values for rxs and rxr, patterned after the hoc code that currently sets up the rx values.

[Qroid_montreal]

Thanks so much. NEURON is a fantastic tool, and a big part of that is how helpful (and patient) this forum is. Thanks again. I'll let you know how it goes.

[Qroid_montreal]

Thanks again for the help. I'm encountering this problem. I'm using Matlab to generate extracellular potential traces as a function of xyz location and time. The extracellular potential for the center of each segment is then saved in a text file, and the data from this text file is then read with Vector.scanvar() and played to e_extracellular. However, from what I can tell by looking at the space plots of e_extracellular, only the last line in the text file is being used and this is being applied to every segment identically. Also, if every value in the file is 0 save the last line I see no effect in NEURON; if I change the last line, though, I see an effect. I can confirm that it's reading the entire file, because if I shorten the file NEURON complains about an unexpected EOF. So, from what I can tell, it's reading in every line but only using the last.

The section of the code that reads the file is below. It's a modified version of "stim.hoc" from the model here http://www.neuron.yale.edu/ftp/ted/neur ... nd_rec.zip.

Code: Select all

// $Id: stim.hoc,v 1.5 2009/02/24 00:55:27 ted Exp ted $
/* Modified on April 17 2012
The stimulus is read from Ver_NEURONformat.txt and copied 
to a Vector.
For each section that has the xtra mechanism, this Vector 
is used to drive e_extracellular.
*/

// extracellular field file
objref f1
f1=new File()
f1.ropen("/Users/Mercy/Academics/Ephaptic/mfiles/Ver_NEURONformat.txt")

// read header variables
dt = f1.scanvar() // set integration dt
NtimeStepsStim = f1.scanvar() // time points in Ve (ie. stimulation)
Nlocxyzr = f1.scanvar

// when to start and stop stimulation?
del = .5 // start at del ms
TimeAfterStim = 10 // time after stimulation
tstop=del + dt*NtimeStepsStim + TimeAfterStim // stop

objref stim_amp, stim_time
stim_amp = new Vector()
stim_time = new Vector()
NtimeSteps = (del + dt*NtimeStepsStim + TimeAfterStim)/dt // total number of time points in simulation
// be careful that NtimeSteps is rounded-to-nearest-integer.
// since resize(n) floors n, use n+0.5.
// NB floor(n+0.5) is equivalent to round(n).
stim_amp.resize(NtimeSteps+0.5)
stim_time.resize(NtimeSteps+0.5)

strdef ThisSector
timecount=0

proc attach_stim() {
  forall {  // LOOP through each section
    if (ismembrane("xtra")) { // check each section to find one that has xtra
      for (x, 0) { // LOOP through each segment
        f1.scanstr(ThisSector) // sector name, from file
        timecount=0
        scancount=0 // confirm that scancount=NtimeStepsStim
        for ii=0,del/dt + NtimeStepsStim + TimeAfterStim/dt - 1 { // LOOP through each time point
          if (ii<=del/dt || ii>del/dt+NtimeStepsStim) {
            stim_amp.x[ii]=0
            stim_time.x[ii]=timecount+dt
            timecount=timecount+dt
          } else {
            stim_amp.x[ii]=f1.scanvar()
            stim_time.x[ii]=timecount+dt
            timecount=timecount+dt // without timecount, would need to reference stim_time[-1] on first iteration
            scancount=scancount+1
          }
        }
        // print scancount
        stim_amp.play(&e_extracellular(x), stim_time)          
      }
    }
  }
}

attach_stim()

Is it possible to say what I'm doing wrong? I don't know how to access the values of stim_amp for each segment as the simulation runs, so I don't know if it is indeed working and I'm plotting it incorrectly, or if it's not working. I can zip up all the hoc, ses and mod files, as well as the txt file from Matlab if that helps. Thanks.