Motor cortex microcircuit simulation based on brain activity mapping (Chadderdon et al. 2014)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:146949
"... We developed a computational model based primarily on a unified set of brain activity mapping studies of mouse M1. The simulation consisted of 775 spiking neurons of 10 cell types with detailed population-to-population connectivity. Static analysis of connectivity with graph-theoretic tools revealed that the corticostriatal population showed strong centrality, suggesting that would provide a network hub. ... By demonstrating the effectiveness of combined static and dynamic analysis, our results show how static brain maps can be related to the results of brain activity mapping."
Reference:
1 . Chadderdon GL, Mohan A, Suter BA, Neymotin SA, Kerr CC, Francis JT, Shepherd GM, Lytton WW (2014) Motor cortex microcircuit simulation based on brain activity mapping. Neural Comput 26:1239-62 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex M1 L2/6 pyramidal intratelencephalic cell; Neocortex fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Oscillations; Laminar Connectivity;
Implementer(s): Lytton, William [billl at neurosim.downstate.edu]; Neymotin, Sam [samn at neurosim.downstate.edu]; Shepherd, Gordon MG [g-shepherd at northwestern.edu]; Chadderdon, George [gchadder3 at gmail.com]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex M1 L2/6 pyramidal intratelencephalic cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
/
src
README
infot.mod *
intf6.mod *
intfsw.mod *
matrix.mod
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
boxes.hoc *
col.hoc
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc *
gcelldata.hoc
gmgs102.nqs
grvec.hoc *
infot.hoc *
init.hoc
intfsw.hoc *
labels.hoc *
load.py
local.hoc *
main.hoc
misc.h *
miscfuncs.py
network.hoc
neuroplot.py *
nload.hoc
nqs.hoc *
nqsnet.hoc
nrnoc.hoc *
params.hoc
run.hoc
samutils.hoc *
saveoutput.hoc
saveweights.hoc
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
wdmaps2.nqs
xgetargs.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
}

Loading data, please wait...