Calculating LFP from multicompartment neuron

[vladimirov]

Hello,
Could anyone give a hint how to efficiently calculate LFP of a multicompartment cell? I know the basic theory, the line source model - take all transmembrane currents and sum them with weights inversely proportional to the distance from electrode. How to implement this in NEURON in a neat way?
Specifically,
a) How to sum up all currents in a cable segment in a computationally efficient and esthetically acceptable way? Should I take each channel name, say NaF, and retrieve it's i_NaF at every time step? Or there is a better way to know a lump transmembrane current from different channels in a segment?
b) After I know total current in each segment at each time step, how can I record it during the simulation? Should I add customary proc advance(), as suggested by Ted Carnevale in another post? : http://www.neuron.yale.edu/phpBB/viewtopic.php?f=12&t=2269&hilit=local+field
Any feedback is greatly appreciated!
Thanks!

[ted]

Could anyone give a hint how to efficiently calculate LFP of a multicompartment cell?

The extracellular mechanism reports total transmembrane current--you can read about it in the Programmer's Reference. Of course there are many steps between that and a complete implementation. For an approach implemented entirely with hoc see Extracellular stimulation and recording in the Hot tips area of The NEURON Forum. The "line source" approach was used by
Gold, C., Henze, D., Koch, C. and Buzsaki, G..
On the origin of the extracellular action potential waveform: a modeling study.
Journal of Neurophysiology 95:3113-3128, 2006 (source code available from ModelDB http://modeldb.yale.edu/).
You may also want to check out the URLs listed in the Forum post you cited, to see if others have made their code available.

How to sum up all currents in a cable segment in a computationally efficient and esthetically acceptable way?

Calculate the weighted sum of membrane currents by iterating over all segments of all sections (forall for (x,0) does this), at each point applying the appropriate weight (actually a transfer resistance) and adding the result to the sum.

After I know total current in each segment at each time step, how can I record it during the simulation?

Use the Vector class's record method. If you really want all that data, right after model setup do the following (expressed here in a mix of pseudocode and hoc):

Code: Select all

objref dvecs
dvecs = new List() // will hold all of the Vectors to which data will be recorded
forall for (x,0) {
  tmpvec = new Vector()
  tmpvec.record(&i_membrane(x)) // records i_membrane of currently accessed segment
    // of currently accessed section
  dvecs.append(tmpvec)
}
objref tmpvec // prevent careless programmer damage to last element in the List

At this point, every simulation run will cause a fresh recording of membrane current throughout the model; specifically, dvec.o(i) will be the Vector that holds the recording of the time course of i_membrane at the ith internal node of the model cell. After a run, just open a file and write all the data to that file (or do whatever else you want with it).

Should I add customary proc advance(), as suggested by Ted Carnevale in another post?

?? Not at all necessary if you just want to record all membrane currents. The thread you cited was about calculating the field in the course of simulation.

[vladimirov]

Many thanks, Ted. Your response is solving the question in full, as always!

[vladimirov]

Ted, I am going through your example code Extracellular stimulation and recording. When I insert xtra mechanism into a minimal hh model, such as

Code: Select all

create soma
soma {
    insert hh
    insert pas
    insert xtra
}

access soma



I get the following error, when I press Run:

Code: Select all

nrniv: Segmentation violation
near line 20
{run()}
        ^
finitialize(-65        )
init(      )
stdinit(    )
run(  )


What do I do wrong? Thank you.

The xtra.mod file:

Code: Select all

NEURON {
   SUFFIX xtra
   RANGE rx, er
   RANGE x, y, z
   GLOBAL is
   POINTER im, ex
}

PARAMETER {
   : default transfer resistance between stim electrodes and axon
   rx = 1 (megohm) : mV/nA
   x = 0 (1) : spatial coords
   y = 0 (1)
   z = 0 (1)
}

ASSIGNED {
   v (millivolts)
   is (milliamp)
   ex (millivolts)
   im (milliamp/cm2)
   er (microvolts)
   area (micron2)
}

INITIAL {
   ex = is*rx*(1e6)
   er = (10)*rx*im*area
: this demonstrates that area is known
: UNITSOFF
: printf("area = %f\n", area)
: UNITSON
}

: Use BREAKPOINT for NEURON 5.4 and earlier
: BREAKPOINT {
:   SOLVE f METHOD cvode_t
: }
:
: PROCEDURE f() {
:   : 1 mA * 1 megohm is 1000 volts
:   : but ex is in mV
:   ex = is*rx*(1e6)
:   er = (10)*rx*im*area
: }

: With NEURON 5.5 and later, abandon the BREAKPOINT block and PROCEDURE f(),
: and instead use BEFORE BREAKPOINT and AFTER BREAKPOINT

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


[ted]

Have you compiled the mod file?
Have you linked the pointer variables im and ex to hoc variables with setpointer statements?

[vladimirov]

Got it, I did not insert the pointers. Now it seems to work. Thanks!

[ehagen]

We've developed a Python module for this; LFPy: [url]compneuro.umb.no/LFPy[/url]

Best,
Espen