Modeling dynamic extracellular fields

[ted]

You're working awfully hard at this. Let me ask you a question, the answer to which might make the task easier.

I'm using Matlab to generate extracellular potential traces as a function of xyz location and time.

Is the extracellular medium purely resistive and linear?

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

Right, because each pass through the loop in proc attach_stim() asserts stim_amp.play(&e_extracellular(x).... So the same Vector (stim_amp) will be played into all e_extracellular instances. The brute force fix is to revise the contents of the loop so that it does this (in pseudocode):

Code: Select all

for each section that has extracellular {
  for each internal node of the current section {
    create a new Vector
    fill this with new data
    play it into the current e_extracellular(x)
    append it to a List
  }
}

The hoc idiom for this is

Code: Select all

objref veclist
veclist = new List() // so any top-level reference to veclist
  // that is executed before the following proc is called
  // does not complain that veclist is not a List
proc whatever() { localobj tmpvec
  forall {
    if (ismembrane("extracellular")) {
      for (x,0) { // omit the zero area nodes (i.e. at 0 and 1)
        tmpvec = new Vector()
        . . . fill tmpvec with new data . . .
        . . . play tmpvec into e_extracellular(x) . . .
        veclist.append(tmpvec)
      }
    }
  }
}

But there may be a better approach, depending on your answer to my question above.

This block of code

Code: Select all

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

attempts to position your stimulus waveform in time, right? It is far more efficient to take advantage of this fact:

On fadvance a transfer will take place if t will be (after the fadvance increment) equal or greater than the associated time of the next index.

(see the Programmer's Reference entry on the Vector class's play() method).

"What does that mean?"

Well, for example, a pulse with amplitude 1 that arises from a baseline of 0 at t = 1 ms, and falls back to 0 at t = 2 ms, is specified by a tvec, vsrc pair with these contents:

Code: Select all

index  tvec  vsrc
  0      0     0
  1      1     0
  2      1     1
  3      2     1
  4      2     0
  5      3     0

This will produce the same waveform regardless of whether dt is 1, 0.5, 0.1, or 0.025 ms (allowing, of course, for roundoff error--unless you are using adaptive integration, in which case the transitions will be at the exact times specified).

So you only need to
read your precomputed waveform into tvec and vsrc
add the "del" parameter to all elements of tvec
insert a "flat" segment at the very start and very end of your stimulus

"What's this 'flat' segment stuff?"

You want the field to be unchanged for the first del ms, right? and you also want it to stop changing after the end of your precomputed waveform, right? If there is a j such that
tvec.x[j+1] > tvec.x[j]
but
vsrc.x[j+1] == vsrc.x[j]
then the extracellular field will remain constant over the time interval tvec.x[j]..tvec.x[j+1]
And if these two points are the last two data points in your "stimulus waveform," Vector play will simply repeat the same vsrc.x value again and again for all t > tvec.x[j+1].

"Great. How do I do that?"

In semi-pseudocode:

Code: Select all

increase the size of both Vectors by 2
rotate the contents of both Vectors one element "to the right"
tvec.x[0] = 0
tvec.x[size-1] = tvec.x[size-2] + 1 // size is the size of the resized tvec and vsrc
vsrc.x[0] = vsrc.x[1] // so the segment from tvec.x[0] to tvec.x[1] is flat
vsrc.x[size-1] = vsrc.x[size-2] + 1 // ditto for the segment from tvec.x[size-2] to tvec.x[size-1]

In a separate thread I already suggested that you rethink these statements
for ii=0,del/dt + NtimeStepsStim + TimeAfterStim/dt - 1 { // LOOP through each time point
if (ii<=del/dt || ii>del/dt+NtimeStepsStim) {
and maybe some others.

edited by NTC at 12:09 PM EDT on 4/22/2012

[Qroid_montreal]

Hi Ted. As always, thanks for your thoughtful reply.

I'm coming back to this problem after a few months away. Briefly, I'm reading in data from an external file and playing it to e_extracellular in each segment. My problem is that while the whole file is being read, only the last line is being used meaningfully. I'm assuming some overwriting is going on. I'm not sure where.

I edited my code as you suggested. Since the data I'm playing is changing on each timestep, I used the "brute force fix". However the problem is still occurring. The code is below.

Code: Select all

// $Id: stim.hoc,v 1.5 2009/02/24 00:55:27 ted Exp ted $
// extracellular field file
objref f1
f1=new File()
f1.ropen("/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
		 
objref veclist
veclist = new List() // so any top-level reference to veclist
  // that is executed before the following proc is called
  // does not complain that veclist is not a List

strdef ThisSector

proc attach_stim() { localobj stim_amp stim_time
  sectcount=0 // how many sections obey ismembrane?
  forall {  // LOOP through each section
    if (ismembrane("extracellular")) { // check each section to find one that has xtra
    sectcount=sectcount+1
      for (x, 0) { // LOOP through each segment
        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)
      
        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
          } else {
            stim_amp.x[ii]=f1.scanvar()
            stim_time.x[ii]=timecount+dt
            scancount=scancount+1
          }
          timecount=timecount+dt // without timecount, would need to reference stim_time[-1] on first iteration
        }
        // print scancount
        stim_amp.play(&e_extracellular(x), stim_time)     
        veclist.append(stim_amp)     
      }
    }
  }
}

attach_stim()

Any ideas what's going wrong?

[ted]

No, but here's how you can figure it out for yourself:
1. Make a toy model with just 1 section that has nseg = 3.
2. Make a toy data set that has just a half dozen points in the stimulus waveform.
3. Embed print statements in your code that associates the stimulus with the model so that you can follow the execution as it proceeds.
4. After all data has been read and all stimulus vectors have been set up, print out all the stimulus vectors and make sure that they are correct.

After you get that working properly,

5. Make a toy model with just 2 sections that each have nseg = 1.
6. Repeat steps 2-4 with this model.

[Qroid_montreal]

I made a toy model and confirmed that the stimulus is being read correctly. Both print statements and the List of stimulus vectors show the stimulus data.

A print statement inside advance such as this:

Code: Select all

proc advance() {
	access SectionName
	print e_extracellular(.5)
	fadvance()
}

shows that only the last accessed segment has had the stimulus 'played' to it, however. Despite reading the stimulus correctly, "print e_extracellular(.5)" returns 0 for all segments except the last. "print e_extracellular(.5)" returns the correct values for the last segment. This leads me to think that subsequent uses of stimulus.play are negating previous uses. I'm confused why this would be the case. Any light you could shed would be valuable. Thanks very much.

Edit: I solved the problem by putting the read and play statements in separate (but otherwise functionally identical) loops. Not sure why this works, but it does. Thanks so much for all your help!

[ted]

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.

Just set up a range variable plot ("space plot") of e_extracellular and use Movie Run or Single step to observe the evolution over time of e_extracellular vs. distance along a path.

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.

In your proc attach_stim(), every pass through the "for loop" executes
stim_amp.play(&e_extracellular(x), stim_time)
but the objref called stim_amp always refers to the same Vector. During simulation setup, the elements in this Vector are overwritten every time you read the time course of extracellular potential for a new segment. The values that are used during a simulation will be the last values that were read. Furthermore, during a simulation these values will be played into every segment's e_extracellular. If you are going to use Vector play to drive each segment's e_extracellular with a different series of values, then you need to associate a different Vector with each segment. If the time course of each segment's e_extracellular is described by N samples, and there are M segments, you'll need M Vectors, each of which contains N elements.

I have three questions for you.

First, why bother testing for xtra when you aren't using xtra to drive e_extracellular?

Second, are you assuming that the extracellular medium is linear and nonreactive (i.e. purely resistive)?

Third, if the answer to my second question is yes, why aren't you using xtra to drive e_extracellular?