/************************************************************
'ca1' model code repository
Written by Ivan Raikov, ivan.g.raikov@gmail.com and Marianne
Bezaire, marianne.bezaire@gmail.com, www.mariannebezaire.com
In the lab of Ivan Soltesz, www.ivansolteszlab.org
Routine 'interpxyz' based on code by Ted Carnevale.
Latest versions of this code are available online at:
ModelDB:
Open Source Brain: http://www.opensourcebrain.org/projects/nc_ca1
This file computes an approximation to the LFP generated by the
pyramidal cells in the network, based on the formula by Schomburg et
al., J Neurosci 2012.
The approximate LFP is calculated as the sum of current contributions
of all compartments, scaled by the distances to the recording
electrode and extracellular medium resistivity. The time resolution
of the LFP calculation may be lower than that of the simulation by
setting lfp_dt in this file.
************************************************************/
strdef s1, s2
objref lfplist, lfp_ids[numCellTypes], lfp_types[numCellTypes], lfpkmatrix[numCellTypes]
lfplist = new List()
//for celltype=0, numCellTypes-1 {
// lfpkmatrix[celltype] = new Matrix()
// lfp_ids[celltype] = new Vector()
// lfp_types[celltype] = new Vector()
//}
func npole_pos_lfp() { local avglfp, vlfp, celltype, pci, gid, i, j, ex, ey, ez, typei, MaxEDist localobj cell
// Calculate the average LFP of select pyramidal cells in the network,
// only including pyramidal cells whose somata are within MaxEDist
// microns of the (x,y,z) location passed into this function as
// arguments: $1 = x, $2 = y, $3 = z in microns
MaxEDist = $1 // microns
ex = mypoint.x[0] //$2
ey = mypoint.x[1] //$3
ez = mypoint.x[2] //$4
sumcell = 0
vlfp = 0 // Initialize the LFP variable
for celltype=0, numCellTypes-1 { // Iterate over all cell types
for i=0, lfp_ids[celltype].size()-1 { // Iterate over all cells chosen for the LFP computation
gid = lfp_ids[celltype].x[i]
cell = pc.gid2cell(gid)
forsec cell.all {
if (ismembrane("extracellular")) {
j = 0
for (x,0) {
vlfp = vlfp + (i_membrane(x) * lfpkmatrix[celltype].x[i][j])
j = j + 1
}
}
}
}
}
// Ensure all processors have completed their calculation before proceeding.
pc.barrier()
avglfp = pc.allreduce(vlfp, 1)
return avglfp
}
/*
Computes xyz coords of nodes in a model cell whose topology & geometry are defined by pt3d data.
Code by Ted Carnevale.
*/
proc interpxyz() { local ii, nn, kk, xr, nsegs localobj xx, yy, zz, xint, yint, zint, ll, range
// get the data for the section
nn = $1
nsegs = $2
xx = $o3
yy = $o4
zz = $o5
ll = $o6
xint = $o7
yint = $o8
zint = $o9
// to use Vector class's .interpolate()
// must first scale the independent variable
// i.e. normalize length along centroid
ll.div(ll.x[nn-1])
// initialize the destination "independent" vector
range = new Vector(nsegs+2)
range.indgen(1/nsegs)
range.sub(1/(2*nsegs))
range.x[0]=0
range.x[nsegs+1]=1
// length contains the normalized distances of the pt3d points
// along the centroid of the section. These are spaced at
// irregular intervals.
// range contains the normalized distances of the nodes along the
// centroid of the section. These are spaced at regular intervals.
// Ready to interpolate.
xint.interpolate(range, ll, xx)
yint.interpolate(range, ll, yy)
zint.interpolate(range, ll, zz)
}
proc setup_npole_lfp() { local i, j, k, m, n, r, h, s, avglfp, typei, gid, celltype, ex, ey, ez, sx, sy, sz, lfptype, ii, nn localobj myr, cell, dists, lst, xx, yy, zz, ll, xint, yint, zint strdef ks, s1, s2
rho = 333.0 // extracellular resistivity, [ohm cm]
fdst = 0.1 // percent of distant cells to include in the computation
ex = mypoint.x[0] //$1
ey = mypoint.x[1] //$2
ez = mypoint.x[2] //$3
//printf ("host %d: entering setup_npole_lfp\n", pc.id)
// Vector for storing longitudinal and perpendicular distances
// between recording electrode and compartments
dists = new Vector(2)
// This loop determines which cells will be used for the LFP computation and the sizes of their compartments
for celltype=0, numCellTypes-1 { // Iterate over all cell types
lfp_ids[celltype] = new Vector()
lfp_types[celltype] = new Vector()
sprint(s1,";%s;",LFPCellTypes)
sprint(s2,";%s;",cellType[celltype].cellType_string)
if (strobj.substr(s1, s2)>-1) {
m = 0
n = 0
//printf ("host %d: npole_lfp: m = %d lfp_ids[%d].size = %d\n", pc.id, m, celltype, lfp_ids[celltype].size())
for pcitr(&i, &typei, &gid, cellType[celltype].cellStartGid, cellType[celltype].cellEndGid) {
cell = pc.gid2cell(gid)
// Relative to the recording electrode position
if (sqrt((cell.x-ex)*(cell.x-ex) + (cell.y-ey)*(cell.y-ey) + (cell.z-ez)*(cell.z-ez))<MaxEDist) {
lfptype = 1 // proximal cell
} else {
myr = ranlfplist.object(i).r
myr.uniform(0, 1)
if (myr.repick() < fdst) {
lfptype = 2 // distal cell -- include only fdst fraction of total
} else {
lfptype = 0
}
}
if (lfptype > 0) {
lfp_ids[celltype].append(gid)
lfp_types[celltype].append(lfptype)
m = m+1
if (n == 0) {
forsec cell.all {
insert extracellular
n = n + nseg
}
}
}
}
if (m > 0) {
lfpkmatrix[celltype] = new Matrix(m,n)
//printf ("host %d: npole_lfp: m = %d lfpkmatrix[%d].nrow = %d lfpkmatrix[%d].ncol = %d lfp_ids[%d].size = %d\n", pc.id, m, celltype, lfpkmatrix[celltype].nrow, celltype, lfpkmatrix[celltype].ncol, celltype, lfp_ids[celltype].size())
}
}
for i=0, lfp_ids[celltype].size()-1 { // Iterates over all cells chosen for the LFP computation. If none of that type, size == 0 and loop won't run
gid = lfp_ids[celltype].x[i]
cell = pc.gid2cell(gid)
// Iterates over each compartment of the cell
j = 0
forsec cell.all {
if (ismembrane("extracellular")) {
nn = n3d()
xx = new Vector(nn)
yy = new Vector(nn)
zz = new Vector(nn)
ll = new Vector(nn)
for ii = 0,nn-1 {
xx.x[ii] = x3d(ii)
yy.x[ii] = y3d(ii)
zz.x[ii] = z3d(ii)
ll.x[ii] = arc3d(ii)
}
xint = new Vector(nseg+2)
yint = new Vector(nseg+2)
zint = new Vector(nseg+2)
interpxyz(nn,nseg,xx,yy,zz,ll,xint,yint,zint)
j = 0
for (x,0) {
sx = xint.x[j]
sy = yint.x[j]
sz = zint.x[j]
// l = L/nseg is compartment length
// r is the perpendicular distance from the electrode to a line through the compartment
// h is longitudinal distance along this line from the electrode to one end of the compartment
// s = l + h is longitudinal distance to the other end of the compartment
l = L/nseg
r = sqrt((ex-sx)*(ex-sx) + (ey-sy)*(ey-sy) + (ez-sz)*(ez-sz))
h = l/2
s = l + h
k = 0.0001 * area(x) * (rho / (4 * PI * l)) * abs(log((sqrt(h*h + r*r) - h) / (sqrt(s*s + r*r) - s)))
sprint(ks,"%g",k)
if ((strcmp(ks,"nan") == 0) || (strcmp(ks,"-nan") == 0)) {
k = 0
}
// Distal cell
if (lfp_types[celltype].x[i] == 2) {
k = (1.0/fdst)*k
}
//printf ("host %d: npole_lfp: gid = %d i = %d j = %d r = %g h = %g k = %g\n", pc.id, gid, i, j, r, h, k)
lfpkmatrix[celltype].x[i][j] = k
j = j + 1
}
}
}
}
}
}
proc sample_npole_lfp() { local avglfp, lfptype localobj vec
ex = mypoint.x[0]
ey = mypoint.x[1]
ez = mypoint.x[2]
// At t=0, calculate distances from recording electrode to all
// compartments of all pyramidal cells, calculate scaling
// coefficients for the LFP calculation, and save them in
// lfpkmatrix.
if (t == 0) {
//setup_npole_lfp(ex,ey,ez)
setup_npole_lfp()
}
// Execute this code regularly throughout the simulation
// to add datapoints to all the LFP recordings (average
// and per-cell).
// Compute LFP across a subset of cells within a certain distance from the recording electrode:
avglfp = npole_pos_lfp(MaxEDist) // x, y, z position of electrode in microns
if (pc.id == 0) {
vec = new Vector()
// For this time step, create a vector with entries of time and average LFP
vec.append(t, avglfp)
// Append the vector for this time step to the list
lfplist.append(vec.c)
}
// Add another event to the event queue, to
// execute sample_npole_lfp again, lfp_dt ms from now
cvode.event(t + lfp_dt, "sample_npole_lfp()")
}