calc_lfp_variable-resistivity.oc

/*----------------------------------------------------------------------------

	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
	http://cns.iaf.cnrs-gif.fr

----------------------------------------------------------------------------*/


/*-------------------------------------------------------------------------

 NOTE ON UNITS: 
 
 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")
load_file("nrngui.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()
	g[ii].size(0,tstop,$2,$3)
	g[ii].xaxis()
	g[ii].yaxis()
	g[ii].addvar($s1,1,0)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
}

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()
	g[ii].size($2,$3,$4,$5)
	g[ii].xaxis()
	g[ii].yaxis()
	g[ii].addvar($s1,1,0)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
}

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
dt=fscan()

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.xaxis()
gI.yaxis()
gI.label("Im (nA)") 

proc draw_Imem() {     // draw membrane current
  gI.beginline(1,1)
  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() {
  C.impedance(fmax,df,rext,rmax,dr,R,sigma1,sigma2,lambda,epsilon1,sigmaR)
  draw_filter()
}

objectvar g1		// draw |Z| vs frequency
g1 = new Graph()
g1.xaxis()
g1.yaxis()
g1.label("|Z|") 

fmaxdraw = 2000		// max frequency to draw

proc draw_filter() {
  g1.beginline(1,1)
  i=0
  zmax=0
  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)
        i=i+1
  }
  g1.size(0,fmaxdraw,0,zmax)
}




//----------------------------------------------------------------------------
//  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
   draw_fft()
}


objectvar g2		// create graph for Power Spectrum
g2 = new Graph()
g2.xaxis()
g2.yaxis()
g2.label("Power") 

proc draw_fft() {       // draw power spectrum of current
  g2.beginline(1,1)
  pmax=0
  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)
   draw_extracellular()
}


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

proc draw_extracellular() {     // draw extracellular V (UNITS converted to mV) 
  g3.beginline(1,1)
  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)
	xpanel()
}

draw_panel()