gcelldata.hoc

// gcelldata.hoc -- routines for processing the empirical data used by the 
//   simulation (gcell = Gordon (Shepherd) cell).

// Globals and routines pertaining to the empirical data used from 
// (Weiler et al., 2008) and other Gordon Shepherd projects.
//
// Last update: 3/12/11 (georgec)

//* Globals
numgcs = 0           // number of Gordon cells (set later)
numgcxs = 16         // number of horizontal bins in Gordon cell conv data
xbinwid = 100        // x bin width (in microns)
numzs = 16           // number of depth bins
zbinwid = 100        // z (depth) bin width (in microns)
lambeg = 0           // the first z bin for a cell layer (set in proc)
lamend = 0           // the last z bin for a cell layer (set in proc)
sigub = -15.0        // significant upper bound (in pA) for the data values
objref wdnq          // (Weiler et al., 2008) data NQS table
objref hupnq         // hookup NQS table
objref convnq        // test cell convergence NQS table (made from hupnq)
strdef gsprepdir     // directory for Gordon Shepherd preprocessed data
gsprepdir = "."
GCELLTYPE_NORMAL = 0      // normal Gordon cell type (otherwise unlabeled)
GCELLTYPE_CXSPINAL = 1    // corticospinal layer 5 cell
GCELLTYPE_CXSTRIATAL = 2  // contralateral corticostriatal layer 5 cell

//* Routines

//** simdatainit() -- load and/or preprocess Gordon cell data
proc simdatainit () { local a,ii,tnumgcs localobj v1,v2,v3
  // Allocate temporary vectors.
  a = allocvecs(v1,v2,v3)

  // Load the NQS tables for Gordon's data.
  {sprint(tstr,"%s/gmgs102.nqs",gsprepdir)}
  rdnqss(tstr)

  // Load the appropriate preprocessed cells data table.
  wdnq = new NQS()
//  {sprint(tstr,"%s/wdmaps.nqs",gsprepdir)} // culled cells with no preprocessing
  {sprint(tstr,"%s/wdmaps2.nqs",gsprepdir)} // culled cells with interpolation applied
  wdnq.rd(tstr) // culled cells with interpolation applied

  // Set the number of Gordon cells.
  numgcs = wdnq.m

  // Make a new NQS table for only cells with neg responses over a cuttoff.
//  wdnq = cullcells(raw, sigub)

  // Create the hookup table.
  mkhuptable()

  // Scale hookup table weights (clip Gordon cell values between the negatives 
  // of the first 2 values, then scale these values to the last 2 values.
  scalehupwts(0,50,0,5)
//  scalehupwts(0,50,0,5)  // original

  // Add new columns to the cells table for (Gordon) cell number, Gordon 
  // cell type, and z bin.
  tnumgcs = cells.size
  v1.resize(tnumgcs)
  v2.resize(tnumgcs)
  v3.resize(tnumgcs)
  v1.indgen(0,tnumgcs-1,1)
  v2.fill(GCELLTYPE_NORMAL)  // all cells are normal Gordon cells for now
  for ii=0,tnumgcs-1 {
    v3.x(ii) = gcellnum2zbin(ii)
  }
  cells.resize("cellnum",v1,"gcelltyp",v2,"zbin",v3)

  // Deallocate the temporary vectors.
  dealloc(a)
}

//** gcellnum2zbin(gordon_cell_type) -- translate Gordon cell number to z 
// (depth) bin
func gcellnum2zbin () { local a,ii,zbin,ysoma localobj vv1
  a = allocvecs(vv1)
 
  // Get the micron distance of the cell from the pia.
  ysoma = cells.v[2].x($1)

  vv1.indgen(zbinwid,zbinwid * numzs,zbinwid)
  zbin = -1
  ii = numzs-1
  while (ii >= 0) {
    if (ysoma <= vv1.x(ii)) {
      zbin = ii  
    }
    ii -= 1
  }
  dealloc(a)
  return zbin
}

//** getlambin(model_cell type) - get the first and last laminar bin for cell 
// type
proc getlambin () { local lyr
  lyr=GetLyr($1)
  if (lyr == 2) {
    lambeg = 1
    lamend = 4
  } else if (lyr == 5) {
    lambeg = 5
    lamend = 10
  } else if (lyr == 6) {
    lambeg = 11
    lamend = 13
  } else {
    print "ERROR: Cannot find laminar bins for cell type #", $1
  }
}

//** cullcells() -- make a new NQS table for cells only above neg. cutoff
obfunc cullcells () { local a,cutoff,cullcount localobj iq,oq,v1
  iq = $o1  // input nq (with 102 cells)
  cutoff = $2   // cutoff value for minimum

  oq = new NQS() // output nq
  a = allocvecs(v1)  // alloc vector for keeping culled cell ID nums

  // Loop over all column vectors of the NQS table...
  for ii=0,iq.m-1 {
    // If the minimum data value is below the cutoff value (e.g. -20 pA)...
    if (iq.v[ii].min < cutoff) { 
      oq.resize(iq.s[ii].s,iq.v[ii])
      v1.append(ii)
    }
  }

  // Copy the cull cell ID numbers into the .x vector of the table
  oq.x.copy(v1)

  // Deallocate vectors.
  dealloc(a)

  return oq
}

//** chkgapcontig() -- check for occurences of places in data where 
//  points were removed because of closeness to soma
proc chkgapcontig () { local a,ii,jj localobj v1,v2
  // Allocate vectors
  a = allocvecs(v1,v2)

  // Loop over all of the cells in the table...
  for ii=0,wdnq.m-1 {
//    print "Cell #", ii

    // Loop over the 16 rows of the raw data matrix...
    for jj=0,15 {
//      printf("Row #%d: ", jj) 
      v1.resize(16)
      v2.resize(16)
      v1.mrow(wdnq.v[ii],jj,16)  // put the jjth row in v1
      v2.indvwhere(v1,"==",1e-9) // pull out indices where number = 1e-9
      if (v2.size > 0) {
        iscontig = v2.ismono(3)    // set true if indices are consecutive
        if (iscontig) { 
//          print "contiguous" 
        } else { 
//          print "NONCONTIGUOUS" 
//          if (v1.min() < -10.0) {
            printf("Cell #%d, Row #%d: ",ii,jj)
            vlk(v1) 
//          }
        }
      } else { 
//        print "no gap"
      }
    }
//    print "------"
  }

  // Dealloc vectors.
  dealloc(a)
}

//** lininterp() -- linear interpolation function, returning a vector
obfunc lininterp () { local val1,val2,veclen,intstep localobj vint
  val1 = $1
  val2 = $2
  veclen = $3

  vint = new Vector(veclen)

  // Calculate a (positive-valued) even step.
  intstep = abs(val2 - val1) / (veclen - 1)  

  // If 2nd value higher, do normal indgen.
  if (val2 > val1) {
    vint.indgen(val1,val2,intstep)
  // If 1st value higher, flip value, do indgen, and flip result.
  } else if (val1 > val2) {
    vint.indgen(val2,val1,intstep)
    vint.reverse()
  // If values equal, just fill with the first value.
  } else {
    vint.fill(val1)
  }

  return vint
}

//** interpgaps() -- interpolate through the "infinity" gaps in the data
proc interpgaps () { local a,ii,jj,kk,val,infbin,srb,erb,intstp localobj v1,v2,v3

  // Allocate vectors
  a = allocvecs(v1,v2,v3)

  // Loop over all of the cells in the table...
  for ii=0,wdnq.m-1 {
//    print "Cell #", ii

    // Loop over the 16 rows of the raw data matrix...
    for jj=0,15 {
//      printf("Row #%d: ", jj) 
      v1.resize(16)
      v2.resize(16)
      v1.mrow(wdnq.v[ii],jj,16)  // put the jjth row in v1
      v2.indvwhere(v1,"==",1e-9) // pull out indices where number = 1e-9
      // Only worry about rows where we have 1e-9...
      if (v2.size > 0) {
        srb = -1
        erb = -1
        for kk=0,v2.size-1 {
          // Get the index of the inf (1e-9) bin.
          infbin = v2.x(kk)

          // Deal with the edge case where bin 0 is 1e-9.
          // Note: if we had runs of more than 1, we'd need to deal with that.
          if (infbin == 0) {
            v1.x(0) = v1.x(1)  // Set to the next bin
          // Deal with the edge case where bin 15 is 1e-9.
          // Note: if we had runs of more than 1, we'd need to deal with that.
          } else if (infbin == 15) {
            v1.x(15) = v1.x(14)  // Set to the previous bin
          // If we have not started a run of Infs yet, start a new run
          } else if (srb == -1) {
            srb = infbin - 1
            erb = infbin + 1
          // If we have started a run and the new bin contig with the run
          } else if (infbin == erb) {
            erb += 1
          // If the next bin is not contiguous with the latest run
          } else {
            // Linearly interpolate for the run.
//            printf("Cell #%d, Row #%d (%d to %d): ",ii,jj,srb,erb)
//            printf("interp between %f and %f\n",v1.x(srb),v1.x(erb))
//            vlk(v2)
            v3 = lininterp(v1.x(srb),v1.x(erb),erb-srb+1)
            v1.copy(v3,srb)  // copy linear interpolation vector back to row
//            vlk(v1)

            // Start a new run.
            srb = infbin - 1
            erb = infbin + 1
          }

          // If we're at the end of the list of inf bins, and have started a 
          // run, output.
          if ((kk == v2.size-1) && (srb != -1)) {
            // Linearly interpolate for the run.
//            printf("Cell #%d, Row #%d (%d to %d): ",ii,jj,srb,erb)
//            printf("interp between %f and %f\n",v1.x(srb),v1.x(erb))
//            vlk(v2)
            v3 = lininterp(v1.x(srb),v1.x(erb),erb-srb+1)
            v1.copy(v3,srb)  // copy linear interpolation vector back to row
//            vlk(v1)
          }
        }

        // Put row v1 back in the NQS data (saving changes).
        v1.msetrow(wdnq.v[ii],jj,16)
      }
    }
//    print "------"
  }

  // Dealloc vectors.
  dealloc(a)
}

//** mkhuptable() -- make a hookup NQS table from wdnq.
proc mkhuptable () { local a,ii,jj,kk,cz localobj vv1
  // Allocate scratch vectors.
  a = allocvecs(vv1)

  // Make the NQS table.
  hupnq = new NQS("srcz","destcell","destxshift","destz","wt")

  // Loop over all of the depth bins...
  for ii=0,numzs-1 {
    // Loop over all "Gordon cells"...
    for jj=0,wdnq.m-1 {
      // Grab row ii, the row corresponding to the depth bin we're at.
      vv1.mrow(wdnq.v[jj],ii,numgcxs)   // numgcxs = 16 now

      // Loop over all x bins in the row...
      for kk=0,vv1.size-1 {
        // If the value is less than the significant upper bound...
        if (vv1.x(kk) < sigub) {
          // Record the entry (source z bin, destination cell #, 
          // destination x shift, destination z bin (with respect to 
          // 7.5), weight).
          cz = gcellnum2zbin(wdnq.x.x(jj))  // get the cell's z bin
          hupnq.append(ii, wdnq.x.x(jj), 7.5-kk, cz, vv1.x(kk))
        }
      } // for kk     
    } // for jj
  } // for ii

  // Deallocate scratch vectors.
  dealloc(a)
}

//** scalehupwts() -- scale the hookup table weights so they are converted 
//  from -pA values to sensible weight values
proc scalehupwts () { local a,ii,gvmin,gvmax,wvmin,wvmax localobj vv1
  gvmin = $1
  gvmax = $2
  wvmin = $3
  wvmax = $4

  // Allocate scratch vectors.
  a = allocvecs(vv1)

  // Read the old Gordon cell weights.
  vv1 = hupnq.getcol("wt")

  // Make the weights positive.
  vv1.mul(-1)
  
  // Clip the Gordon cell weights between gvmin and gvmax.
  for ii=0,vv1.size-1 {
    if (vv1.x(ii) < gvmin) vv1.x(ii) = gvmin
    if (vv1.x(ii) > gvmax) vv1.x(ii) = gvmax
  }

  // Scale between wvmin and wvmax.
  vv1.scale(wvmin,wvmax)

  // Set the new weights.
//  hupnq.setcol("wt",vv1)

  // Deallocate scratch vectors.
  dealloc(a)
}

//** huptable2cellconv() -- test function for showing that huptable has 
// all significant convergence info
obfunc huptable2cellconv () { local a,ii,jj localobj v1
  // Allocate vectors.
  a = allocvecs(v1)

  // Create the new empty convergence table.
  convnq = new NQS()
  convnq.cp(wdnq)  // copy from original table
  for ii=0,convnq.m-1 {  // zero out all data
    convnq.v[ii].fill(0)
  }

  // Loop over all "Gordon" cells...
  for ii=0,numgcs-1 {
    // Select just the rows of the hupnq corresponding to the desired cell.
    hupnq.select("destcell",wdnq.x.x(ii))

    // Loop over all found rows for the desired cell...
    for jj=0,hupnq.size-1 {
      v1 = hupnq.getrow(jj)
      
      // Set the convergence matrix for cell ii.
      convnq.v[ii].mset(v1.x(0),int(v1.x(2)-7.5),numgcxs,v1.x(4))
    } // for jj
  } // for ii

  // Toggle back to the full table for hupnq.
  hupnq.tog

  // Deallocate scratch vectors.
  dealloc(a)

  return convnq
}