Modeling local field potentials (Bedard et al. 2004)

 Download zip file   Auto-launch 
Help downloading and running models
This demo simulates a model of local field potentials (LFP) with variable resistivity. This model reproduces the low-pass frequency filtering properties of extracellular potentials. The model considers inhomogeneous spatial profiles of conductivity and permittivity, which result from the multiple media (fluids, membranes, vessels, ...) composing the extracellular space around neurons. Including non-constant profiles of conductivity enables the model to display frequency filtering properties, ie slow events such as EPSPs/IPSPs are less attenuated than fast events such as action potentials. The demo simulates Fig 6 of the paper.
1 . Bédard C, Kröger H, Destexhe A (2004) Modeling extracellular field potentials and the frequency-filtering properties of extracellular space. Biophys J 86:1829-42 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Connectionist Network; Extracellular;
Brain Region(s)/Organism:
Cell Type(s):
Gap Junctions:
Simulation Environment: NEURON;
Model Concept(s): Extracellular Fields;
Implementer(s): Destexhe, Alain [Destexhe at];

	Neuron demo file to simulate local field potentials
        Model: simple EPSP/IPSP model with two synaptic currents,
	(monopolar current source) => Read data file with source current

        Method: calculate the extracellular field potentials based on
        inhomogeneous conductivity in extracellular space in radial symetry
	(Bedard, Kroger & Destexhe model)

	  - compute the Fourier transform of the source current I
          - compute the impedance Z assuming a non-homogeneous conductivity
            function (radial symetry) - MECHANISM IN MOD FILE
          - compute the extracellular voltage by inverse Fourier transform
            of the product Z(w) * I(w)

	THIS FILE: simulates fig 6 of the paper

        For details, see:

        Bedard C, Kroger H & Destexhe A.  Modeling extracellular field 
        potentials and the frequency-filtering properties of extracellular 
        space.  Biophysical Journal 86: 1829-1842, 2004.

        Written by A. Destexhe, 2002



 Fourier transform of current: 
  I(w) = INT dt exp(-iwt) I(t)
  (A-s) =   (s)  (none)   (A) 
  (nA-s) =  (s)  (none)   (nA)
		=> units of nA-s (otherwise too small numbers)
		=> time in seconds, frequency in Hz
 Fourier transform of voltage: 
  V(w) = Z(w) I(w)
  (V-s) = (Ohm) (A-s)
		=> Z(w) must be in units of Ohm (and not Ohm-s)
 Fourier transform of impedance:
  Z(w)  = 1 / (4pi sigmaR) INT dr/r^2 (sigR-iw epsR)/(sig(r)-iw eps(r))
  (Ohm) = 1 / (S/m)          (m)/(m^2)   (none)
  (Ohm) = 1 / (S/m)          (1/m)
  (Ohm) = 1 / (S/um)         (1/um)
		=> distances can be specified in microns, if sigmaR in S/um
		=> impedance in Ohm
 Inverse Fourier transform of voltage:
  Vext = INT dw exp(iwt) [Z(w) I(w)]
  (V)  =    1/s  (none)  (Ohm A s = V-s)
  (nV) =    1/s  (none)  (Ohm nA-s = nV-s)
 => all frequencies in Hz, time in seconds, current in nA
 => extracellular voltage in nV

 (see other notes labeled by "UNITS" below)

//  load and define general graphical procedures

// original code commented for below replacement
// xopen("$(NEURONHOME)/lib/hoc/stdrun.hoc")

objectvar g[20]			// max 20 graphs
ngraph = 0

proc addgraph() { local ii	// define subroutine to add a new graph
				// addgraph("variable", minvalue, maxvalue)
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()

proc makegraph() { local ii	// define subroutine to add a new graph
				// makeraph("variable", xmin,xmax,ymin,ymax)
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()

nrnmainmenu()			// create main menu
nrncontrolmenu()		// crate control menu

// Read in data points of source current

ropen("curr.dat")	// source current of EPSP/spike model

npoints = fscan()	// nb points should be an exponent of 2

objectvar Curr		// create vectors of data points
Curr = new Vector(npoints)

for i=0,npoints-1 { 	// read source currents
 Curr.set(i, fscan() )

//  create graphs

tmaxdraw = 10           // set time max to draw

objectvar gI		// create graph for membrane current
gI = new Graph()
gI.label("Im (nA)") 

proc draw_Imem() {     // draw membrane current
  for i=0, npoints-1 {
    if(i*dt <= tmaxdraw) {
        gI.line( i*dt, Curr.get(i))
  gI.size(0, tmaxdraw, Curr.min(), Curr.max())

//  Calculate and store the impedance

// UNITS: compute the impedance in Ohm

npt = npoints/2			// nb of points in frequency 
df = 1000/(npoints*dt)		// freq interval in Hz
fmax = df*npt			// max frequency in Hz

R = 105				// radius of source (UNITS: um)
rext = 500			// distance of the recording electrode (um)
rmax = 200*R			// max distance for integration (um)
dr = rmax/1000                  // integration step (um)

lambda = 500			// space cst. of conductivity decay (UNITS=um)
sigma1 = 1			// high-amplitude of conductivity (normalized)
sigma2 = 0.2			// low-amplitude of conductivity (normalized)
				// (this is the asymptotic value, macroscopic)
sigma2 = 2e-9			// low-amplitude of conductivity (normalized)
				// (this is the minimal value, membranes)
epsilon1 = 6e-11  		// value of permittivity (constant, normalized)

sigmaR = 1.56			// absolute value of conductivity at the source
                                // UNITS: S/m
sigmaR = sigmaR * 1e-6          // converted to UNITS of S/um so that distances
				// can be specified in um

epsilon1 = 0.0001  // heuristic

create soma		// create fake compartment for impedance
access soma
objref C
C = new ImpedanceFM(0.5)	// create object that will contain the C code
soma C.loc(0.5)		// locate the object in the compartment

double Zr[npt+1],Zi[npt+1]
setpointer C.vec1, Zr[0]	// assign pointers to vectors
setpointer C.vec2, Zi[0]

// Calculate filter:
//  loop over frequencies -> calculate impedances -> store in vectors
proc calc_filter() {

objectvar g1		// draw |Z| vs frequency
g1 = new Graph()

fmaxdraw = 2000		// max frequency to draw

proc draw_filter() {
  for(f=df; f<=fmaxdraw; f=f+df) {			// draw impedance
        ReZ = Zr[i]
        ImZ = Zi[i]
        z = sqrt(ReZ*ReZ+ImZ*ImZ)               // norm of impedance
	if(z>zmax) zmax=z
	g1.line(f, z)

//  Procedures to calculate FFT (from Vector objects)

// UNITS: current in nA => FFT (Ir,Ii) will be in nA-s

objref Ir, Ii, Vr, Vi, Vext		// create vectors
Ir = new Vector(npt+1)		// Fourier transform of the current
Ii = new Vector(npt+1)
Vr = new Vector(npt+1)		// Fourier transform of extracellular voltage
Vi = new Vector(npt+1)
Vext = new Vector(npoints)      // Extracellular voltage

objref RFT			// Vector for FFT
RFT = new Vector(npoints)

// fft function from Vector in ivoc:
// The "fft" function is taken from the "realft" procedure from Numerical
// Recipes (Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. 
// Numerical Recipes. The Art of Scientific Computing. Cambridge University
// Press, 1986).  It inputs a real data vector of N points (N must be exponent
// of 2) and returns the half of its Fourier transform (the other half being
// symetric), as complex numbers.  The first and second points returned are
// the first and last real-valued points of the FFT.  The other points are
// returned as couples of real/complex values.  The frequency interval is 
// delta-f = 1/(N*dt).  To yield the correct amplitudes of the FT, all real
// (cosine) and complex (sine) transforms must be divided by (N/2), except the 
// first and second points returned which must be divided by N.
// Inverse Fourier transform: use fft with second argument of "-1".

//  Procedure to calculate FFT of source current
proc runfft() {
   draw_Imem()				// draw source current
   RFT.fft(Curr)			// calculate FFT of current
   for(i=0; i<npt; i=i+1) {     	// decompose real/complex components
        Ir.set(i, RFT.get(2*i)/npt )	// store into Ir, Ii
        Ii.set(i, RFT.get(2*i+1)/npt )
   Ir.set(npt, RFT.get(1)/npt)		// last real component is 2nd pt of FFT
   Ii.set(0, 0)				// first complex component set to zero
   Ii.set(npt, 0)			// last complex component set to zero

objectvar g2		// create graph for Power Spectrum
g2 = new Graph()

proc draw_fft() {       // draw power spectrum of current
  for i=0,npt {
    if(i*df<=fmaxdraw) {
        pow = Ir.get(i)^2 + Ii.get(i)^2
	if(pow > pmax) pmax=pow
        g2.line(i*df, pow)
  g2.size(0, fmaxdraw, 0, pmax)

//  Procedure to calculate extracellular voltage at a given distance
//  (assumes that Ir, Ii contains the Re/Im values of current; UNITS nA-s)
//  (assumes that Zr, Zi contains the Re/Im values of impedance; UNITS Ohm)
proc calc_extracellular() {
   for i=0, npt {		// compute the product (Zr,Zi)*(Ir,Ii)
	a = Zr[i]*Ir.get(i) - Zi[i]*Ii.get(i)
        b = Zr[i]*Ii.get(i) + Zi[i]*Ir.get(i)
        Vr.set(i, a)		// store in Vr, Vi (w-freq component of Vext)
        Vi.set(i, b)            // (UNITS: nV-s)
   for i=0, npt-1 {		// store into RFT vector for inverse FFT
        RFT.set(2*i, Vr.get(i) )
        RFT.set(2*i+1, Vi.get(i) )
   RFT.set(1, Vr.get(npt) )     // 2nd point is last real component

   RFT.set(0, 0 )		// remove dc (nonzero dc is due to num error)

   Vext.fft(RFT,-1)             // compute inverse FFT and store into Vext
                                // Vext is the extracellular voltage (UNITS nV)

objectvar g3		// create graph for Vext
g3 = new Graph()
g3.label("Vext (mV)") 

proc draw_extracellular() {     // draw extracellular V (UNITS converted to mV) 
  for i=0, npoints-1 {
    if(i*dt <= tmaxdraw) {
        g3.line( i*dt, Vext.get(i) * 1e-6 )
  g3.size(0, tmaxdraw, Vext.min()*1e-6, Vext.max()*1e-6)

//  Create command panel

proc draw_panel() {
	xpanel("Calculate extracellular field potentials")
	xpvalue("Position of electrode (um)",&rext)
	xpvalue("Max distance for integration",&rmax)
	xpvalue("Integration step",&dr)
	xpvalue("Spatial decay of conductivity",&lambda)
	xpvalue("Minimal conductivity (normalized)",&sigma2)
	xpvalue("Conductivity at the source (absolute)",&sigmaR)
	xpvalue("Value of permittivity (normalized)",&epsilon1)
	xbutton("=> 1. Calculate impedances","calc_filter()")
	xbutton(" => Redraw impedances","draw_filter()")
	xbutton("=> 2. Calculate current FFT","runfft()")
	xbutton(" => Redraw FFT","draw_fft()")
	xbutton("=> 3. Calculate extracellular potential","calc_extracellular()")
        xbutton(" => Redraw extracellular potential","draw_extracellular()")
	xpvalue("Max frequency to draw",&fmaxdraw)
	xpvalue("Max time to draw",&tmaxdraw)