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 pyramidal corticothalamic L6 cell; Neocortex M1 pyramidal intratelencephalic L2-5 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 pyramidal corticothalamic L6 cell; Neocortex M1 pyramidal intratelencephalic L2-5 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 *
                            
// $Id: spkts.hoc,v 1.86 2010/07/10 02:32:11 samn Exp $

print "Loading spkts.hoc..."

// ancillary programs for handling vectors

{load_file("decvec.hoc")}
{load_file("decnqs.hoc")}


//* transfer a file into a list of strings
// usage 'f2slist(list,file)'
proc f2slist() { local i
  $o1.remove_all
  if (! tmpfile.ropen($s2)) { print "Can't open ",$s2
    return }
  while (tmpfile.gets(temp_string_)>0) {
    sfunc.head(temp_string_,"\n",temp_string_) // chop
    tmpobj = new String(temp_string_)
    $o1.append(tmpobj)
  }
}

//* spkts(attrnum[,flag,min,max]) graph spikes from the vectors

thresh = -20    // threshold for deciding which are spikes
burstlen = 1.4  // duration of spike or burst, don't accept another till after this time

proc spkts_call() {}// callback function stub
proc spkts () { local cnt, attrnum, ii, pstep, jj, num, time0, flag, min, max, tst
  revec(vec,vec1) // store times and indices respectively
  if (numarg()==0) { print "spkts(attrnum[,flag,min,max])\n\tflag 0: graph, flag 1: save vec1,vec to veclist, flag 2: save inds (vec1) and times (vec)" // infomercial
    return }
  attrnum=$1
  panobj = GRV[attrnum]
  if (attrnum==0) { cnt=printlist.count() } else { cnt = panobj.llist.count() }
  pstep = panobj.printStep
  if (numarg()>1) { flag = $2 } else { flag = 0 }
  if (numarg()>2) { min = $3 } else { min = 0 }
  if (numarg()>3) { max = $4 } else { max = cnt-1 }
  if (flag==0){
    newPlot(0,1,0,1)
    panobj.glist.append(graphItem)
  }
  for ii=min,max {
    if (attrnum==0) { 
      vrtmp.copy(printlist.object(ii).vec) 
      if (panobj.printStep==-2) tvec = printlist.object(ii).tvec
      if (panobj.printStep==-1) tvec = panobj.tvec
    } else {
      panobj.rv_readvec(ii,vrtmp)  // pick up vector from file
      if (panobj.printStep<0) tvec = panobj.tvec
    }
    if (panobj.printStep>=0) { // make a tvec
      if (!isobj(tvec,"Vector")) { print "ERR0 spkts(): tvec not a vector" return }
      tvec.resize(vrtmp.size)
      tvec.indgen(pstep)
    }
    spkts1(tvec,vrtmp,vec,vec1,ii)

    if (panobj.printStep<0) { // a tvec
      tst = vec.max
    } else {
      tst = pstep*vrtmp.size()         // calc the tst
    }
  }
  if (flag==1) { savevec(vec1) savevec(vec) }
  if (flag<=0) {
    vec1.mark(graphItem,vec,"O",panobj.line)  // graph all the times
    printf("\n")
    graphItem.size(0,tst,min,max)
    graphItem.xaxis(0)
    graphItem.label(0.1,0.9,panobj.filename)
  }
}

// spkts1(tvec,vec,timev,indexv,index)
proc spkts1 () { local ind,tm,ix
  ind=tm=allocvecs(2) tm+=1
  spkts_call()  // place to reset thresh or do other manipulations
  mso[tm].resize($o2.size)
  mso[tm].xing($o2,$o1,thresh) // times
  if (numarg()==5) {
    mso[ind].resize(mso[tm].size)  // scratch vector stores index
    mso[ind].fill($5)
    $o3.append(mso[tm])     // add the times for this to end of vec
    $o4.append(mso[ind])  // add same index for each spk to end of vec1
  } else {
    vlk(mso[tm])
  }
  dealloc(ind)
}

//* parse_spkts()
// pull the vec and vec1 files from spkts apart and put in alloc'ed vectors
func parse_spkts () { local p,q
  p=allocvecs(vec1.max+2) q=p+1
  for (ii=0;ii<=vec1.max;ii+=1) {
    mso[p].indvwhere(vec1,"==",ii)
    mso[q].index(vec,mso[p]) 
    q += 1
  }
  return p+1
}

proc line_spkts () { local ii,min,max,xmax,skip
  skip = $1
  if (numarg()==3) { min=$2 max=$3 } else {
    min = int(graphItem.size(3)+1) max = int(graphItem.size(4)) }
  xmax = graphItem.size(2)
  for (ii=min;ii<max;ii=ii+skip) {
    graphItem.beginline()
    graphItem.line(0,ii)
    graphItem.line(xmax,ii)
  }
  graphItem.xaxis()
}

burst_time=0
burst_maxfreq = 30
calc_ave = 0

//** calcspkts(flag,index)
// run after spkts() so vec contains the times, vec1 contains the
// indices
proc calcspkts () { local ii,jj,flag,index,p1,p2,mn,mx
  p1 = allocvecs(2,1000) p2 = p1+1
  if (numarg()==0) {
    print "To be run after spkts(). \
Assumes times in vec, indices in vec1. \
calcspkts(flag,min,max)\nflags:\t1\tspk times\n\t2\tspk freq \
\t3\tburst times\n\t4\tburst freq\nset calc_ave to get averages for freqs"
    return
  }
  // vec contains the times, vec1 contains the indices
  flag = $1
  mn = $2
  if (numarg()==3) { mx=$3 } else { mx=mn }
  for index=mn,mx {
    mso[p2].resize(0)
    mso[p1].indvwhere(vec1,"==",index)
    mso[p1].index(vec,mso[p1])
    if (flag==1) {  
      printf("SPKS for #%d: ",index)
      for jj=0,mso[p1].size()-1 {printf("%g ",mso[p1].x[jj])}
      printf("\n")
    } else if (flag==2) {  
      printf("FREQ for #%d: ",index)
      for jj=0,mso[p1].size()-2 { 
        pushvec(mso[p2],1000./(mso[p1].x[jj+1]-mso[p1].x[jj])) }
      if (calc_ave) { print mso[p2].mean } else { vlk(mso[p2]) }
    } else if (flag==3) {  
      printf("BTIMES for #%d: ",index)
      burst_time = mso[p1].x[0]
      for jj=1,mso[p1].size()-1 {
        if (1000./(mso[p1].x[jj]-burst_time) < burst_maxfreq) {
          printf("%g ",burst_time)
          burst_time = mso[p1].x[jj]
        } 
      }
      printf("\n")
    } else if (flag==4) {  
      printf("BFREQ for #%d: ",index)
      burst_time = mso[p1].x[0]
      for jj=1,mso[p1].size()-1 {
        // should keep track of spike times in case of very long bursts
        if (1000./(mso[p1].x[jj]-burst_time) < burst_maxfreq) {
          pushvec(mso[p2],1000./(mso[p1].x[jj]-burst_time))
          burst_time = mso[p1].x[jj]
        } 
      }
      if (calc_ave) { print mso[p2].mean } else { mso[p2].printf }
    }
  }
  dealloc(p1)
}

func rvwheres () { local ii
  if ($1!=0) {
    for ii=0,panobjl.object($1).llist.count()-1 {
      if (sfunc.substr(panobjl.object($1).llist.object(ii).name,$s2)==0) {
        return ii }
    }
    errorMsg("String not found in rvwheres.")
  }
  return -2
}

supind = 0
//* spkhist assume spk times in vec 
// allows superimposing of graphs
// spkhist(bin_size)
proc spkhist () { local ii,jj,min,max,diff
  if (numarg()==0) { print "spkhist(bin_size)" return }
  if (numarg()==3) { min=$2 max=$3 } else { min=0 max=tstop }
  diff = max-min
  vrtmp.hist(vec,min,max,$1)
  vec0.resize(4*diff/$1)
  vec1.resize(4*diff/$1)
  vec0.fill(0) vec1.fill(0)
  for (ii=min;ii<int(diff/$1);ii=ii+1) {
    jj=ii*4
    vec0.x[jj+0] = ii*$1
    vec0.x[jj+1] = ii*$1
    vec0.x[jj+2] = (ii+1)*$1
    vec0.x[jj+3] = (ii+1)*$1
    vec1.x[jj+0] = 0
    vec1.x[jj+1] = vrtmp.x[ii]
    vec1.x[jj+2] = vrtmp.x[ii]
    vec1.x[jj+3] = 0
  }
  if (panobj.super==0) {
    newPlot(min,max,0,vrtmp.max)
    panobj.glist.append(graphItem)  
  } else { graphItem = panobjl.object(panobj.remote).glist.object(supind) 
    supind = supind+1 }
  vec1.line(graphItem,vec0)
  sprint(temp_string_,"Hist: %s %d",panobj.filename,$1)
  graphItem.label(0.1,0.9,temp_string_)
}

//** truncvec (vec1,margin) 
// truncate a thresholded time vector so that only one time is given for each spike
// vec1 has thresholded times, margin is duration of a spike
proc truncvec () { local a,ii,num,marg,time0 localobj vs
  marg = $2
  a=allocvecs(vs)
  num=0 time0=-1e3
  vs.resize($o1.size())
  vs.fill(-2)
  for ii=0,$o1.size()-1 {
    if ($o1.x[ii] > time0+marg) { 
      vs.x[ii] = $o1.x[ii]
      time0 = $o1.x[ii]
    }
  }
  $o1.where(vs,">",-1)
  dealloc(a)
}

//** redundkeep(vec) keeps sequential redundent entries
// destructive
proc redundkeep () { local x,ii
  $o1.sort
  x = $o1.x[0]
  for ii=1,$o1.size-1 {
    if ($o1.x[ii]!=x) { $o1.x[ii-1]=-1e20 x=$o1.x[ii] }
  }
  $o1.where($o1,">",-1e20)
}

//** after running spkall can see which cells are responsible for spks
// assumes spk times in vec, spk identities in vec1
// uses ind and vec0
proc whichspked () { local ii
  ind.indvwhere(vec,"()",$1,$2) // a range
  vec0 = vec1.ind(ind)
  ind = vec.ind(ind)
  for ii=0,ind.size()-1 { printf("%d %g\n",vec0.x[ii],ind.x[ii]) }
}

// firebtwn(ind,time,min,max) list of cells that fire between times min and max
proc firebtwn () { local ii,p1,p2,p3
  p1 = allocvecs(3) p2=p1+1 p3=p2+1
  mso[p3].indvwhere($o2,"[]",$3,$4)
  mso[p1].index($o1,mso[p3]) // indices
  mso[p2].index($o2,mso[p3]) // times
  printf("%d hits\n",mso[p3].size)
  for vtr2(&x,&y,mso[p1],mso[p2]) {
      printf("%4d::%6.2f ",x,y)
      if ((ii+1)%5==0) { print "" }
  }
  print ""
  dealloc(p1)
//  dealloc(p2) // to save the indexes
}

// elimind(ind,time,min,max) take out cells with nums between min,max
// destructive
proc elimind () { local ii,p1
  p1 = allocvecs(1)
  mso[p1].indvwhere($o1,"[]",$3,$4)
  vecelim($o1,mso[p1]) vecelim($o2,mso[p1])
  dealloc(p1)
}

// index/time graph
// tigr(ind,vec,size,marker)
proc tigr () { local sz
  if (numarg()==0) { print "tigr(Yvec,Xvec,marker size,marker type)" 
    print "Marker types: \"o\",t,s,O,T,S,+ (circ, tri, square; CAP is filled)"
    return }
  if (numarg()>2) { sz=$3 } else { sz=6 }
  if (numarg()>3) { temp_string_=$s4 } else { temp_string_="O" }
  nvplt($o2)
  graphItem.size($o2.min,$o2.max,$o1.min,$o1.max)
  $o1.mark(graphItem,$o2,temp_string_,sz,panobj.curcol) 
}

//* p2nqs(#,panobj,nqs) -- copy an entry into an nqs
proc p2nqs () { local x,a localobj v1,q,p
  x=$1 p=$o2 q=$o3
  q.resize(0)
  a=allocvecs(v1)
  p.rv_readvec(x,v1)
  q.resize("time",p.tvec,"ind",v1)
  dealloc(a)
}

//** spkboth() determines how many cells spike in 2 time periods
proc spkboth () { local a,t1,t2,t3,t4,s1,s2,s3 localobj v1,v2,v3,o
  o=$o1 t1=$2 t2=$3 t3=$4 t4=$5
  printf("MAY NEED DEBUGGING SINCE NQS.getcol() CHANGED\n")
  a=allocvecs(v1,v2,v3,1e4)
  o.verbose=0
  o.select("time","()",t1,t2) v1.redundout(o.getcol("ind"))
  o.select("time","()",t3,t4) v2.redundout(o.getcol("ind"))
  v3.insct(v1,v2)
  s1=v1.size s2=v2.size s3=v3.size
  printf("P1: %d, P2: %d, Both: %d (%d%%, %d%%)\n",s1,s2,s3,s3/s1*100,s3/s2*100)
  o.verbose=1
  dealloc(a)
}

//returns NQS containing ID,Type,SpikeT
//doesn't check if cell is dead or alive, assumes input is valid
//$o1 = spike vitem
//or
//$o1 = spike vitem, $2 == skipI
//or
//$o1 = time vec , $o2 = id vec
obfunc SpikeNQS(){ local idx,skipI localobj vec,tvec,nq
  if(ce==nil) return nil
  skipI=0
  if(numarg()==1){
    vec = $o1.vec tvec = $o1.tvec
  } else if(numarg()==2){
    if(argtype(1)==1 && argtype(2)==1){
      tvec=$o1 vec=$o2
    } else {
      vec = $o1.vec tvec=$o1.tvec skipI=$2
    }
  } else {
    printf("SpikeNQS ERRA: invalid args!\n")
    return nil
  }
  nq = new NQS("ID","Type","SpikeT")
  if(skipI){
    for idx=0,2 { nq.v[idx].resize(vec.size) nq.v[idx].resize(0) }
    for idx=0,vec.size-1 {
      if(ice(ce.o(vec.x(idx)).type)) continue
      nq.append(vec.x(idx),ce.o(vec.x(idx)).type,tvec.x(idx))
    }
  } else {
    nq.v[0].copy(vec)
    nq.v[2].copy(tvec)
    nq.v[1].resize(vec.size)
    for idx=0,vec.size-1 nq.v[1].x(idx)=ce.o(vec.x(idx)).type
  }
  return nq
}

//returns NQS with refractory % of cell types vs time -- assumes all cells of a type
//have the same refractory period
//$o1 = nqs from SpikeNQS, $2 = dt, optional, $3=skip inhib cells, optional
obfunc refracNQ () { local ct,tt,dt,s,skipI,dotypes localobj snq,nr,vid
  snq=$o1  nr=new NQS("Type","t","r")  ct=0
  if(numarg()>1)dt=$2 else dt=0.25
  if(numarg()>2)skipI=$3 else skipI=1
  if(numarg()>3)dotypes=$4 else dotypes=0
  for(tt=0;tt<=tmax_INTF;tt+=dt){
    for ctt(&ct) if(skipI && !ice(ct)) {
      if(snq.select("Type",ct,"SpikeT","[]",tt-ce.o(ix[ct]).refrac,tt)){
        vid=snq.getcol("ID")//after select, so will use output
        s=vid.uniq
      } else {
        s=0
      }
      nr.append(ct,tt,s/numc[ct])
    }
  }
  return nr
}

//returns nqs with % of cells of each type that have activated by time=t
obfunc PActNQS () { local tinc,winsz,idx,ct,tt,tm,spks,cells,spksE,spksI,nE,nI,cellsE,cellsI\
                   localobj nqt,va,vspk,snq
  snq=$o1
  //time start,end,cell type,activated:0-1,spikes,cells:abs,cells:0-1
  nqt=new NQS("ts","te","ct","act","spks","cells","cellsn")
  if(numarg()>1)tinc=$2 else tinc=0.25
  if(numarg()>2)winsz=$3 else winsz=0.5
  va=new Vector(CTYPi+1)//keep track of % of cells of a type that have spiked
  va.fill(0)
  {vspk=new Vector(allcells) vspk.fill(0)}//keep track of which cells have spiked
  snq.verbose=0
  snq.tog("DB") tm=snq.getcol("SpikeT").max
  for(tt=0;tt<tm;tt+=tinc) {
    spksE=spksI=nE=nI=cellsE=cellsI=0
    for ctt(&ct) {
      if((spks=snq.select("Type",ct,"SpikeT","[]",tt,tt+winsz))) {
        for vtr(&idx,snq.getcol("ID"))vspk.x(idx)=1 
        if(ice(ct)) spksI+=spks else spksE+=spks
      }
      va.x(ct) = vspk.sum(ix[ct],ixe[ct]) / numc[ct]
      if(spks>0) cells=snq.out.getcol("ID").uniq else cells=0
      if(ice(ct)) {
        nI+=va.x(ct)*numc[ct]
        cellsI+=cells
      } else {
        nE+=va.x(ct)*numc[ct]
        cellsE+=cells
      }
      nqt.append(tt,tt+winsz,ct,va.x(ct),spks,cells,cells/numc[ct])
    }
    nqt.append(tt,tt+winsz,-1,nE/ecells,spksE,cellsE,cellsE/ecells)
    nqt.append(tt,tt+winsz,-2,nI/icells,spksI,cellsI,cellsI/icells)
  }
  snq.verbose=1
  return nqt
}

//$o1=nqs from PActNQS , gets peaks & intervals of cell activity levels,
//ct == -1 for E cells, ct == -2 for I cells
obfunc GetPeakNQ () { local idx,ct localobj vi,nqp,vc,nqpo,nqin
  nqin=$o1 nqin.tog("DB") vc=new Vector(nqin.size(-1)) nqin.getcol("ct").uniq(vc)
  nqin.verbose=0
  nqpo=new NQS("ct","ts","x","y","dx","dy")
  for vtr(&ct,vc) {
    nqin.select("ct",ct) vi=nqin.getcol("spks")
    nqp=new NQS("ct","ts","x","y") 
    for idx=1,vi.size-2 {
      if(vi.x(idx)>vi.x(idx-1) && vi.x(idx)>vi.x(idx+1)) {
        nqp.append(ct,nqin.getcol("ts").x(idx),idx,vi.x(idx))
      }
    }
    nqp.resize("dx") nqp.resize("dy")
    nqp.v[nqp.m-2]=Deriv(nqp.getcol("x"))
    nqp.v[nqp.m-1]=Deriv(nqp.getcol("y"))
    nqpo.append(nqp)
    nqsdel(nqp)
  }
  nqin.verbose=1
  return nqpo
}


// returns list containing spike times for each cell
// $o1 == raster vitem
obfunc spikelist () { local idx localobj ls,vec,tvec
  vec=$o1.vec tvec=$o1.tvec
  ls=new List()
  for idx=0,allcells-1 ls.append(new Vector())
  for idx=0,vec.size-1 ls.o(vec.x(idx)).append(tvec.x(idx))
  return ls
}

// plot a single cell's spike times
// $o1 == snq, $2 == cell id, $3 == color, $4 == size, $5==drawR
proc plotcellst () { local cid,st,clr,a,sz,drawR localobj snq,vt,vx,vy
  snq=$o1 cid=$2 clr=$3 sz=$4
  if(numarg()>4)drawR=$5 else drawR=0
  a=allocvecs(vx,vy)
  snq.select("ID",cid)
  gvmarkflag=1
  vt=snq.out.v[2]//SpikeT
  for vtr(&st,vt) vx.append(st) vy.append(cid)
  if(drawR) for vtr(&st,vt) drline(st+.05,cid,st+ce.o(cid).refrac,cid,g,1,1)
  vy.mark(g,vx,"O",sz,clr,1)
  dealloc(a)
}

//draw fancier raster
//$o1=nqs from SpikeNQS, $2==draw refractory periods, $3==skip inhib cells
proc drawrastw () { local idx,skipI,drawR,maxID,a,c,sz,drlt localobj snq,vx,vy,vtype,vc,vtu
  a=allocvecs(vx,vy) drlt=drlflush drlflush=0
  vc=new Vector(CTYP.count+1) vc.fill(0)
  snq=$o1 snq.tog("DB") maxID=snq.v[0].max//max ID
  vtype=new Vector() vtype.copy(snq.v[1])//Type
  if(numarg()>1)drawR=$2 else drawR=0
  if(numarg()>2)skipI=$3 else skipI=0
  if(numarg()>3)sz=$4 else sz=2
  vtu=new Vector(vtype.size) vtype.uniq(vtu)
  c=2
  for vtr(&idx,vtu) if(!skipI || !ice(idx)) {
    vc.x(idx)=c
    c+=1
  }
  if(skipI){
    for idx=0,maxID if(!ice(ce.o(idx).type)) {
      plotcellst(snq,idx,vc.x(ce.o(idx).type),sz,drawR)
    }
  } else {
    for idx=0,maxID {
      plotcellst(snq,idx,vc.x(ce.o(idx).type),sz,drawR)
    }
  }
  dealloc(a) drlflush=drlt  
  if(name_declared("rasterlines"))rasterlines()
  g.flush
}



// plot inhib cells in rast
// $o1 == snq, $2==color, $3==size
proc plotIrast () { local idx,clr,sz localobj xo,snq
  snq=$o1 clr=$2 sz=$3 idx=0
  for ltr(xo,ce,&idx) if(ice(xo.type)) plotcellst(snq,idx,clr,sz)
}

//simple coefficient of variation of interspike interval synch measure from tiesinga03.pdf
//$o1 = spike nqs from SpikeNQS()
//$2 = interval time
//$3 = slide time
//$4 = skip inhib cells [optional] default == 1
obfunc CVPNQS(){ local idx,startt,endt,midt,N,intt,slidet,ct,skipI,CVp localobj snq,cvpnq,vs,vi,vu
  snq=$o1 intt=$2 slidet=$3 if(numarg()>3)skipI=$4 else skipI=1
  vs=new Vector(allcells*2) vi=new Vector(allcells*2)
  vs.resize(0) vi.resize(0) snq.verbose=0
  cvpnq=new NQS("Type","startt","endt","midt","sync","N","CVp","sync2")
  startt=0 endt=intt midt=intt/2
  vu=new Vector(allcells)
  for(startt=0;startt<=tmax_INTF+1-intt;startt+=slidet){
    endt=startt+intt
    if(endt>=tmax_INTF) endt=tmax_INTF
    midt=(startt+endt)/2
    if( (N=snq.select("SpikeT",">=",startt,"SpikeT","<",endt)) > 2 ){
      if(N>vu.size)vu.resize(N)
      snq.out.getcol("ID").uniq(vu) N=vu.size //# of active cells
      vs.copy(snq.out.getcol("SpikeT")) //spike times
      vs.sort //sort spike times to make ISI for all active cells
      vi.resize(0)
      for idx=0,vs.size-2 vi.append(vs.x(idx+1)-vs.x(idx))
      CVp=vi.stdev/vi.mean
      cvpnq.append(0,startt,endt,midt,(CVp-1.)/sqrt(N),N,CVp,(CVp-1.)/sqrt(vi.size))
    } else {
      cvpnq.append(0,startt,endt,midt,0,0,0,0)
    }
    for ctt(&ct) {
      if(skipI && ice(ct)) continue
      if( (N=snq.select("Type",ct,"SpikeT",">=",startt,"SpikeT","<",endt)) > 2 ){
        if(N>vu.size)vu.resize(N)
        snq.out.getcol("ID").uniq(vu) N=vu.size //# of active cells
        vs.copy(snq.out.getcol("SpikeT")) //spike times
        vs.sort //sort spike times to make ISI for all active cells
        vi.resize(0)
        for idx=0,vs.size-2 vi.append(vs.x(idx+1)-vs.x(idx))
        CVp=vi.stdev/vi.mean
        cvpnq.append(ct,startt,endt,midt,(CVp-1.)/sqrt(N),N,CVp,(CVp-1.)/sqrt(vi.size))
      } else {
        cvpnq.append(ct,startt,endt,midt,0,0,0,0)
      }
    }
  }
  snq.verbose=1
  return cvpnq
}

//return spike frequency NQS
//$o1=spike nqs from SpikeNQS()
//$2=interval [optional]
//$3=just do types with interval [optional]
//$4=skipI [optional] default==1
//$5=stop time [optional]
//$6=start time [optional]
//$7=slide time [optional]
obfunc FreqNQS(){ local idx,startt,endt,intt,ct,dotypes,starttime,stoptime,sp,slidet,skipI\
                 localobj fnq,snq
  snq=$o1 intt=$2
  if(numarg()>2) dotypes=$3 else dotypes=0
  if(numarg()>3) skipI=$4 else skipI=1
  if(numarg()>4) stoptime=$5 else stoptime=tstop
  if(numarg()>5) starttime=$6 else starttime=0
  if(numarg()>6) slidet=$7 else slidet=intt
  if(!dotypes){
    fnq=new NQS("ID","Type","Freq","StartT","EndT")
    intt=$2 startt=starttime endt=startt+intt
    for(;startt<stoptime;startt+=slidet){
      endt=startt+intt
      //check length of interval, make sure it's within time bounds of run
      if(endt >= stoptime) endt = stoptime
      if(startt>=endt)endt=startt+1
      for idx=0,allcells-1{
        if(skipI && ice(ce.o(idx).type))continue
        sp=snq.select(-1,"ID",idx,"SpikeT",">=",startt,"SpikeT","<",endt)
        fnq.append(idx,ce.o(idx).type,1e3*sp/(endt-startt),startt,endt)
      }
    }
  } else {
    fnq=new NQS("Type","Freq","StartT","EndT")
    intt=$2 startt=starttime endt=startt+intt
    for(;startt<stoptime;startt+=slidet){
      endt=startt+intt
      //check length of interval, make sure it's within time bounds of run
      if(endt >= stoptime) endt = stoptime
      if(startt>=endt)endt=startt+1
      for ctt(&ct) { 
        if(skipI && ice(ct))continue
        sp=snq.select(-1,"Type",ct,"SpikeT",">=",startt,"SpikeT","<",endt)
        fnq.append(ct,1e3*sp/(numc[ct]*(endt-startt)),startt,endt)
      }
    }
  }
  return fnq
}


func binfindtidx () { local done,val,idx,m,lo,hi,t localobj vv
  vv=$o1 t=$2
  lo=0  hi=vv.size-1  m=int(vv.size/2) done=0
  while(!done){
    if(vv.x(m)>t){
      hi=m m=int((hi+lo)/2.0)
    } else if(vv.x(m)<t){
      lo=m m=int((hi+lo)/2.0)
    }
    if(hi==m || lo==m) return m
  }
}

obfunc estconmat () { local t1,t2,idx,jdx,kdx,del,st1,st2,tdx,df localobj ls,emat,vs1,vs2,vp
  ls=$o1 del=$2  emat=new List()
  for idx=0,ls.count-1 emat.append(new Vector(ls.count))
  estconmat_vc(ls,del,emat)
  return emat
}

// geteff -- get efficiency of $1 to $2 connections
// $1 == type 1 (from)
// $2 == type 2 (to)
// $o3 == ls from spikelist
// returns vector with access by index
obfunc geteff () { local idx,jdx,kdx,ldx,ty1,ty2,tot,cnt,tt,del,fctr,cntonce\
              localobj vt,vpo,vdel,vm,vt2,vtmp,ls,vid
  ty1=$1  ty2=$2  ls=$o3 
  if(numarg()>3) fctr=$4 else fctr=2
  if(numarg()>4) cntonce=$5 else cntonce=0
  vpo=new Vector(allcells)
  vdel=new Vector(allcells)
  vm=new Vector(allcells)
  vtmp=new Vector(allcells)
  vtmp.resize(0)
  for idx=ix[ty1],ixe[ty1]{
    if(idx%100==0)printf("%d.",idx)
    if(!ls.o(idx).size) continue
    vt=ls.o(idx)
    ce.o(idx).getdvi(vpo,vdel)
    tot=0
    for jdx=0,vpo.size-1 if(ce.o(vpo.x(jdx)).type==ty2) tot+=1
    if(!tot) continue
    vtmp.resize(0)
    for jdx=0,vt.size-1 { // go thru spikes
      tt=vt.x(jdx)
      cnt=0    
      for kdx=0,vpo.size-1{ // go thru outputs
        if(ce.o(vpo.x(kdx)).type!=ty2) continue
        del=vdel.x(kdx)
        vt2=ls.o(vpo.x(kdx))
        for ldx=0,vt2.size-1{ // check spiketimes of each output
          if(vt2.x(ldx) >= tt+del && vt2.x(ldx) <= tt+fctr*del){ // within range?
            cnt += 1
            if(cntonce) break // only count once? then break
          }
        }
      }
      vtmp.append(cnt/tot)
    }
    if(vtmp.size) vm.x(idx)=vtmp.mean
  }
  printf("\n")
  return vm
}

// get efficiency of excitatory connections between populations in a given path
// specified by $o1 , i.e. $o1 = new Vector(layer4,layer2,layer5,layer6,layer4)
// $o2 == ls , from spikelist(...)
// $3 == delay fctr for geteff
// returned as NQS
obfunc getpatheffnq () { local fctr,cntonce localobj vpath,nqp,ve,ls
  vpath=$o1
  ls=$o2
  fctr=$3
  if(numarg()>3) nqp=$o4 else nqp=new NQS("ID","from","to","e")
  if(numarg()>4) cntonce=$5 else cntonce=0
  for idx=0,vpath.size-2 {
    printf("eff from %s to %s: ",CTYP.o(vpath.x(idx)).s,CTYP.o(vpath.x(idx+1)).s)
    ve=geteff(vpath.x(idx),vpath.x(idx+1),ls,fctr,cntonce)
    for jdx=ix[vpath.x(idx)],ixe[vpath.x(idx)]{
      nqp.append(jdx,vpath.x(idx),vpath.x(idx+1),ve.x(jdx))
    }
  }
  return nqp
}

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]

References and models cited by this paper

References and models that cite this paper

Ainsworth M, Lee S, Cunningham MO, Roopun AK, Traub RD, Kopell NJ, Whittington MA (2011) Dual gamma rhythm generators control interlaminar synchrony in auditory cortex. J Neurosci 31:17040-51 [PubMed]

Anderson CT, Sheets PL, Kiritani T, Shepherd GM (2010) Sublayer-specific microcircuits of corticospinal and corticostriatal neurons in motor cortex. Nat Neurosci 13:739-44 [PubMed]

Apicella AJ, Wickersham IR, Seung HS, Shepherd GM (2012) Laminarly orthogonal excitation of fast-spiking and low-threshold-spiking interneurons in mouse motor cortex. J Neurosci 32:7021-33 [PubMed]

Bastos AM, Usrey WM, Adams RA, Mangun GR, Fries P, Friston KJ (2012) Canonical microcircuits for predictive coding. Neuron 76:695-711 [PubMed]

Beltramo R, D'Urso G, Dal Maschio M, Farisello P, Bovetti S, Clovis Y, Lassi G, Tucci V, De P (2013) Layer-specific excitatory circuits differentially control recurrent network dynamics in the neocortex. Nat Neurosci 16:227-34 [PubMed]

Bernhardt BC, Chen Z, He Y, Evans AC, Bernasconi N (2011) Graph-theoretical analysis reveals disrupted small-world organization of cortical thickness correlation networks in temporal lobe epilepsy. Cereb Cortex 21:2147-57 [PubMed]

Binzegger T, Douglas RJ, Martin KA (2004) A quantitative map of the circuit of cat primary visual cortex. J Neurosci 24:8441-53 [PubMed]

Bureau I, Shepherd GM, Svoboda K (2004) Precise development of functional and anatomical columns in the neocortex. Neuron 42:789-801 [PubMed]

Cardin JA, Carlen M, Meletis K, Knoblich U, Zhang F, Deisseroth K, Tsai LH, Moore CI (2009) Driving fast-spiking cells induces gamma rhythm and controls sensory responses. Nature 459:663-7 [PubMed]

Carnevale NT, Hines ML (2006) The NEURON Book

Chadderdon GL, Neymotin SA, Kerr CC, Lytton WW (2012) Reinforcement learning of targeted movement in a spiking neuronal model of motor cortex. PLoS One 7:e47251-57 [PubMed]

Chen D, Fetz EE (2005) Characteristic membrane potential trajectories in primate sensorimotor cortex neurons recorded in vivo. J Neurophysiol 94:2713-25 [PubMed]

Dantzker JL, Callaway EM (2000) Laminar sources of synaptic input to cortical inhibitory interneurons and pyramidal neurons. Nat Neurosci 3:701-7 [PubMed]

Dembrow NC, Chitwood RA, Johnston D (2010) Projection-specific neuromodulation of medial prefrontal cortex neurons. J Neurosci 30:16922-37 [PubMed]

Douglas RJ, Martin KA (2004) Neuronal circuits of the neocortex. Annu Rev Neurosci 27:419-51 [PubMed]

Freeman LC (1977) A set of measures of centrality based on betweenness Sociometry 40:35-41

Hattox AM, Nelson SB (2007) Layer V neurons in mouse cortex projecting to different targets have distinct physiological properties. J Neurophysiol 98:3330-40 [PubMed]

Hooks BM, Hires SA, Zhang YX, Huber D, Petreanu L, Svoboda K, Shepherd GM (2011) Laminar analysis of excitatory local circuits in vibrissal motor and sensory cortical areas. PLoS Biol 9:e1000572 [Journal] [PubMed]

   Laminar analysis of excitatory circuits in vibrissal motor and sensory cortex (Hooks et al. 2011) [Model]

Hooks BM, Mao T, Gutnisky DA, Yamawaki N, Svoboda K, Shepherd GM (2013) Organization of cortical and thalamic input to pyramidal neurons in mouse motor cortex. J Neurosci 33:748-60 [PubMed]

Isomura Y, Harukuni R, Takekawa T, Aizawa H, Fukai T (2009) Microcircuitry coordination of cortical motor information in self-initiation of voluntary movements. Nat Neurosci 12:1586-93 [PubMed]

Kameda H, Hioki H, Tanaka YH, Tanaka T, Sohn J, Sonomura T, Furuta T, Fujiyama F, Kaneko T (2012) Parvalbumin-producing cortical interneurons receive inhibitory inputs on proximal portions and cortical excitatory inputs on distal dendrites. Eur J Neurosci 35:838-54 [PubMed]

Katzel D, Zemelman BV, Buetfering C, Wölfel M, Miesenböck G (2011) The columnar and laminar organization of inhibitory connections to neocortical excitatory cells. Nat Neurosci 14:100-7 [PubMed]

Kerr C, Van_albada S, Neymotin S, Chadderdon G, Robinson P, Lytton W (2012) Effects of basal ganglia on cortical computation: A hybrid network-field model Soc. Neurosci. Abstracts 301.16

Kerr CC, Neymotin SA, Chadderdon GL, Fietkiewicz CT, Francis JT, Lytton WW (2012) Electrostimulation as a prosthesis for repair of information flow in a computer model of neocortex IEEE Transactions on Neural Systems & Rehabilitation Engineering 20(2):153-60 [Journal] [PubMed]

   Prosthetic electrostimulation for information flow repair in a neocortical simulation (Kerr 2012) [Model]

Kiritani T, Wickersham IR, Seung HS, Shepherd GM (2012) Hierarchical connectivity and connection-specific dynamics in the corticospinal-corticostriatal microcircuit in mouse motor cortex. J Neurosci 32:4992-5001 [PubMed]

Lang S, Dercksen VJ, Sakmann B, Oberlaender M (2011) Simulation of signal flow in 3D reconstructions of an anatomically realistic neural network in rat vibrissal cortex. Neural Netw : [PubMed]

Lefort S, Tomm C, Floyd Sarria JC, Petersen CC (2009) The excitatory neuronal network of the C2 barrel column in mouse primary somatosensory cortex. Neuron 61:301-16 [PubMed]

Lytton W, Stewart M (2006) Rule-based firing for network simulations Neurocomputing 69:1160-1164

Lytton WW, Neymotin SA, Hines ML (2008) The virtual slice setup. J Neurosci Methods 171:309-15 [Journal] [PubMed]

   The virtual slice setup (Lytton et al. 2008) [Model]

Lytton WW, Omurtag A (2007) Tonic-clonic transitions in computer simulation. J Clin Neurophysiol 24:175-81 [PubMed]

   Tonic-clonic transitions in a seizure simulation (Lytton and Omurtag 2007) [Model]

Lytton WW, Omurtag A, Neymotin SA, Hines ML (2008) Just in time connectivity for large spiking networks Neural Comput 20(11):2745-56 [Journal] [PubMed]

   JitCon: Just in time connectivity for large spiking networks (Lytton et al. 2008) [Model]

Lytton WW, Stewart M (2005) A rule-based firing model for neural networks Int J Bioelectromagn 7:47-50

Markram H (2006) The blue brain project. Nat Rev Neurosci 7:153-60 [Journal] [PubMed]

   [241 reconstructed morphologies on NeuroMorpho.Org]

Miller MN, Okaty BW, Nelson SB (2008) Region-specific spike-frequency acceleration in layer 5 pyramidal neurons mediated by Kv1 subunits. J Neurosci 28:13716-26 [PubMed]

Morgan RJ, Soltesz I (2008) Nonrandom connectivity of the epileptic dentate gyrus predicts a major role for neuronal hubs in seizures. Proc Natl Acad Sci U S A 105:6179-84 [Journal] [PubMed]

   Dentate gyrus (Morgan et al. 2007, 2008, Santhakumar et al. 2005, Dyhrfjeld-Johnsen et al. 2007) [Model]

Morishima M, Kawaguchi Y (2006) Recurrent connection patterns of corticostriatal pyramidal cells in frontal cortex. J Neurosci 26:4394-405 [PubMed]

Neymotin S, Kerr C, Francis J, Lytton W (2011) Training oscillatory dynamics with spike-timing-dependent plasticity in a computer model of neocortex Signal Processing in Medicine and Biology Symposium (SPMB), IEEE :1-6

Neymotin SA, Jacobs KM, Fenton AA, Lytton WW (2011) Synaptic information transfer in computer models of neocortical columns. J Comput Neurosci. 30(1):69-84 [Journal] [PubMed]

   Synaptic information transfer in computer models of neocortical columns (Neymotin et al. 2010) [Model]

Neymotin SA, Lazarewicz MT, Sherif M, Contreras D, Finkel LH, Lytton WW (2011) Ketamine disrupts theta modulation of gamma in a computer model of hippocampus Journal of Neuroscience 31(32):11733-11743 [Journal] [PubMed]

   Ketamine disrupts theta modulation of gamma in a computer model of hippocampus (Neymotin et al 2011) [Model]

Neymotin SA, Lee H, Park E, Fenton AA, Lytton WW (2011) Emergence of physiological oscillation frequencies in a computer model of neocortex. Front Comput Neurosci 5:19-75 [Journal] [PubMed]

   Emergence of physiological oscillation frequencies in neocortex simulations (Neymotin et al. 2011) [Model]

Oviedo HV, Bureau I, Svoboda K, Zador AM (2010) The functional asymmetry of auditory cortex is reflected in the organization of local cortical circuits. Nat Neurosci 13:1413-20 [PubMed]

Packer AM, Yuste R (2011) Dense, unspecific connectivity of neocortical parvalbumin-positive interneurons: a canonical microcircuit for inhibition? J Neurosci 31:13260-71 [PubMed]

   [37 reconstructed morphologies on NeuroMorpho.Org]

Roopun AK, Kramer MA, Carracedo LM, Kaiser M, Davies CH, Traub RD, Kopell NJ, Whittington MA (2008) Period concatenation underlies interactions between gamma and beta rhythms in neocortex. Front Cell Neurosci 2:1-11 [PubMed]

Sakata S, Harris KD (2009) Laminar structure of spontaneous and sensory-evoked population activity in auditory cortex. Neuron 64:404-18 [PubMed]

Schubert D, Staiger JF, Cho N, Kotter R, Zilles K, Luhmann HJ (2001) Layer-specific intracolumnar and transcolumnar functional connectivity of layer V pyramidal cells in rat barrel cortex. J Neurosci 21:3580-92 [PubMed]

Song W, Kerr CC, Lytton WW, Francis JT (2013) Cortical plasticity induced by spike-triggered microstimulation in primate somatosensory cortex. PLoS One 8:e57453-18 [PubMed]

Stepanyants A, Chklovskii DB (2005) Neurogeometry and potential synaptic connectivity. Trends Neurosci 28:387-94 [PubMed]

Stepanyants A, Hirsch JA, Martinez LM, Kisvarday ZF, Ferecskó AS, Chklovskii DB (2008) Local potential connectivity in cat primary visual cortex. Cereb Cortex 18:13-28 [PubMed]

   [10 reconstructed morphologies on NeuroMorpho.Org]

Suter BA, Migliore M, Shepherd GM (2013) Intrinsic electrophysiology of mouse corticospinal neurons: a class-specific triad of spike-related properties. Cereb Cortex 23:1965-77 [PubMed]

Tanaka YH, Tanaka YR, Fujiyama F, Furuta T, Yanagawa Y, Kaneko T (2011) Local connections of layer 5 GABAergic interneurons to corticospinal neurons. Front Neural Circuits 5:12

Thomson AM, Bannister AP (2003) Interlaminar connections in the neocortex. Cereb Cortex 13:5-14 [PubMed]

Tiesinga P, Sejnowski TJ (2009) Cortical enlightenment: are attentional gamma oscillations driven by ING or PING? Neuron 63:727-32 [PubMed]

van Diessen E, Hanemaaijer JI, Otte WM, Zelmann R, Jacobs J, Jansen FE, Dubeau F, Stam CJ, Go (2013) Are high frequency oscillations associated with altered network topology in partial epilepsy? Neuroimage 82:564-73 [PubMed]

Vierling-Claassen D, Cardin JA, Moore CI, Jones SR (2010) Computational modeling of distinct neocortical oscillations driven by cell-type selective optogenetic drive: separable resonant circuits controlled by low-threshold spiking and fast-spiking interneurons. Front Hum Neurosci 4:198 [Journal] [PubMed]

   Engaging distinct oscillatory neocortical circuits (Vierling-Claassen et al. 2010) [Model]

Wang XJ (2010) Neurophysiological and computational principles of cortical rhythms in cognition. Physiol Rev 90:1195-268 [PubMed]

Wang Y, Markram H, Goodman PH, Berger TK, Ma J, Goldman-Rakic PS (2006) Heterogeneity in the pyramidal network of the medial prefrontal cortex. Nat Neurosci 9:534-42 [PubMed]

   [204 reconstructed morphologies on NeuroMorpho.Org]

Weiler N, Wood L, Yu J, Solla SA, Shepherd GM (2008) Top-down laminar organization of the excitatory network in motor cortex. Nat Neurosci 11:360-6 [Journal] [PubMed]

   Laminar connectivity matrix simulation (Weiler et al 2008) [Model]

Wilke C, Worrell G, He B (2011) Graph analysis of epileptogenic networks in human partial epilepsy. Epilepsia 52:84-93 [PubMed]

Wilson HR, Cowan JD (1972) Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J 12:1-24 [Journal] [PubMed]

   Excitatory and inhibitory interactions in populations of model neurons (Wilson and Cowan 1972) [Model]

Yu J, Anderson CT, Kiritani T, Sheets PL, Wokosin DL, Wood L, Shepherd GM (2008) Local-Circuit Phenotypes of Layer 5 Neurons in Motor-Frontal Cortex of YFP-H Mice. Front Neural Circuits 2:6-405 [PubMed]

Dura-Bernal S, Li K, Neymotin SA, Francis JT, Principe JC, Lytton WW (2016) Restoring behavior via inverse neurocontroller in a lesioned cortical spiking model driving a virtual arm. Front. Neurosci. Neuroprosthetics 10:28 [Journal]

   Cortical model with reinforcement learning drives realistic virtual arm (Dura-Bernal et al 2015) [Model]

Dura-Bernal S, Neymotin SA, Kerr CC, Sivagnanam S, Majumdar A, Francis JT, Lytton WW (2017) Evolutionary algorithm optimization of biological learning parameters in a biomimetic neuroprosthesis. IBM Journal of Research and Development (Computational Neuroscience special issue) 61(2/3):6:1-6:14 [Journal]

   Motor system model with reinforcement learning drives virtual arm (Dura-Bernal et al 2017) [Model]

Neymotin SA, Dura-Bernal S, Lakatos P, Sanger TD, Lytton WW (2016) Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex. Front Pharmacol 7:157 [Journal] [PubMed]

   Multitarget pharmacology for Dystonia in M1 (Neymotin et al 2016) [Model]

Neymotin SA, McDougal RA, Bulanova AS, Zeki M, Lakatos P, Terman D, Hines ML, Lytton WW (2016) Calcium regulation of HCN channels supports persistent activity in a multiscale model of neocortex Neuroscience 316:344-366 [Journal] [PubMed]

   Ca+/HCN channel-dependent persistent activity in multiscale model of neocortex (Neymotin et al 2016) [Model]

(64 refs)