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

 Download zip file   Auto-launch 
Help downloading and running models
Accession:138379
"Coordination of neocortical oscillations has been hypothesized to underlie the "binding" essential to cognitive function. However, the mechanisms that generate neocortical oscillations in physiological frequency bands remain unknown. We hypothesized that interlaminar relations in neocortex would provide multiple intermediate loops that would play particular roles in generating oscillations, adding different dynamics to the network. We simulated networks from sensory neocortex using 9 columns of event-driven rule-based neurons wired according to anatomical data and driven with random white-noise synaptic inputs. ..."
Reference:
1 . 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 [PubMed]
Citations  Citation Browser
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 GLU cell; Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex V1 interneuron basket PV GABA cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Oscillations; Synchronization; Laminar Connectivity;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org];
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic GLU cell; Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex V1 interneuron basket PV GABA cell; GabaA; AMPA; NMDA; Gaba; Gaba; Glutamate;
/
fdemo
readme.txt
intf6_.mod
misc.mod *
nstim.mod *
stats.mod *
vecst.mod
col.hoc
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc
finish_run.hoc
grvec.hoc *
init.hoc *
labels.hoc *
local.hoc *
misc.h
mosinit.hoc
network.hoc
nload.hoc
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc *
params.hoc
python.hoc *
pywrap.hoc *
run.hoc
setup.hoc
simctrl.hoc *
spkts.hoc *
stats.hoc *
syncode.hoc *
xgetargs.hoc *
                            
// $Id: nload.hoc,v 1.189 2011/02/16 15:46:37 samn Exp $ 

strdef strvecs,strnclnq //vecs file, vq nqs file, nqin nqs file, nqfld nqs file
strdef strv
strv="09apr17.18"
sprint(tstr,"o[%d]",numcols)
declare("nqtt",tstr,"nqte",tstr,"nqtea",tstr,"nqktau",tstr,"nclnq",tstr,"vLFP",tstr,"glfp","o[1]")
declare("vspudx",new Vector(),"vspudy",new Vector())
declare("nqCO",tstr,"nqisi",tstr,"nqCTY",tstr)
if(name_declared("snq")) {
  nqsdel(snq)
  objref snq[numcols]
} else declare("snq",tstr)
if(name_declared("fnq")) {
  nqsdel(fnq)
  objref fnq[numcols]
} else declare("fnq",tstr)
if(name_declared("vq")) {
  nqsdel(vq)
  objref vq[numcols]
} else declare("vq",tstr)
declare("CDX",0)
sprint(tstr,"o[%d][10]",numcols)
declare("nqIN",tstr)
declare("skipsnq",1) // iff==1 dont get snq and isinq
declare("binsz",100) //100ms
{pflg=1 ers=0}
declare("htmin",0)//sgrdel //min time in histogram
declare("htmax",mytstop)  //max time in histogram
for i=0,numcols-1 vq[i]=lcol.o(i).cstim.vq

//* VQIsiNQS(vq,sy) - get an NQS with ISI from external inputs stored in vq
// only uses synapses type specified by sy
obfunc VQIsiNQS() { local idx,tt,jdx,sy,a localobj nqisi,vid,vd,vt,vq
  vq=$o1 sy=$2 a=allocvecs(vid,vd)
  {nqisi=new NQS("ind","visi") nqisi.odec("visi")}
  {vq.verbose=nqisi.verbose=0 vq.tog("DB")}
  vid.resize(vq.v.size) vq.getcol("ind").uniq(vid)
  for vtr(&idx,vid) if(vq.select("ind",idx,"sy",sy)>2) {
    vt=vq.getcol("time")
    vd.resize(0) vd.deriv(vt,1,1)
    nqisi.append(idx,vd)
  }
  {vq.tog("DB") vq.verbose=nqisi.verbose=1 dealloc(a)}
  return nqisi
}

//* getLV(visi) -- get LV (local variation)
func getLV () { local i,j,s,n,d localobj vec
  vec=$o1 s=0
  if(vec.size < 2) {printf("getLV ERRA: input vec too small!\n") return 0}
  for i=0,vec.size-2 {
    n = (vec.x(i)-vec.x(i+1))^2
    d = (vec.x(i)+vec.x(i+1))^2
    s += 3*n/d
  }
  return s / (vec.size-1)
}

//* addCVcol(nqisi) - add column to input NQS with CV -- input NQS has to have a visi objref col
// containing a Vector of interspike intervals
proc addCVcol () {  local i,c,m,a localobj nqi,vec
  a=allocvecs(vec) nqi=$o1 nqi.tog("DB")
  if((c=nqi.fi("cv")) == -1){ nqi.resize("cv") nqi.pad() c=nqi.m-1}
  for i=0,nqi.v.size-1 {
    vec.resize(0) vec.copy(nqi.get("visi",i).o)
    m=vec.mean()
    if(m) nqi.v[nqi.m-1].x(i) = vec.stdev / m else nqi.v[nqi.m-1].x(i)=0
  }
  dealloc(a)
}

//* addLVcol(nqisi) - add column to input NQS with LV (local variability of interspike interval) --
// input NQS has to have a visi objref col containing a Vector of interspike intervals
// LV is 1 for poisson ISI, 0 for regular
proc addLVcol () {  local i,c,m,a localobj nqi,vec
  a=allocvecs(vec) nqi=$o1 nqi.tog("DB")
  if((c=nqi.fi("lv")) == -1){ nqi.resize("lv") nqi.pad() c=nqi.m-1}
  for i=0,nqi.v.size-1 {
    {vec.resize(0) vec.copy(nqi.get("visi",i).o)}
    nqi.v[nqi.m-1].x(i) = vec.lvar()
  }
  dealloc(a)
}

//* IsiNQS() - makes an NQS with ISI,CV,LV information for each cell, uses globals: CDX,col,snq!!!
obfunc IsiNQS () { local idx,tt,jdx,a localobj nqisi,vid,vd,vt,vty,vic
  if(snq[CDX]==nil)snq[CDX]=SpikeNQS(vit[CDX])
  {nqisi=new NQS("gid","id","type","ice","visi") nqisi.odec("visi")}
  snq[CDX].verbose=nqisi.verbose=0
  a=allocvecs(vid,vd)
  for idx=0,col[CDX].allcells-1 {
    if(snq[CDX].select("id",idx)>2) {
      vt=snq[CDX].getcol("t")
      vty=snq[CDX].getcol("type")
      vic=snq[CDX].getcol("ice")
      vd.resize(0) vd.deriv(vt,1,1) // ISI
      nqisi.append(col[CDX].ce.o(idx).gid,idx,vty.x(0),vic.x(0),vd)
    }
  }
  addCVcol(nqisi) // coefficient of variation
  addLVcol(nqisi) // local variation
  snq[CDX].tog("DB")
  dealloc(a)
  snq[CDX].verbose=nqisi.verbose=1
  return nqisi
}

//* SpikeNQS (time vec, id vec, cell list) returns NQS containing gid,id,type,ice,t<-spike times
obfunc SpikeNQS (){ local idx,gcol,tycol,icecol,idcol,tcol,colcol localobj vec,tvec,nq,ce,c
  {tvec=$o1 vec=$o2 ce=$o3.ce}
  nq=new NQS("col","gid","id","type","ice","t") 
  {gcol=nq.fi("gid") tycol=nq.fi("type") icecol=nq.fi("ice") idcol=nq.fi("id") tcol=nq.fi("t") colcol=nq.fi("col")}
  {nq.v[idcol].copy(vec) nq.v[tcol].copy(tvec) nq.pad()}
  for idx=0,vec.size-1 { c=ce.o(vec.x(idx))
    nq.v[colcol].x(idx)=c.col
    nq.v[gcol].x(idx)=c.gid()
    nq.v[tycol].x(idx)=c.type
  }
  {nq.v[icecol].copy(nq.v[tycol]) nq.v[icecol].apply("ice")}
  return nq
}

//* FreqNQS(snq,ce,[,duration])
//get an NQS with rate of firing of each cell
//$o1 = spikenqs, $o2 = cell list, $3 = duration (starts from 0)
obfunc FreqNQS () { local i,a,ns,tt localobj snq,fnq,vi,ce
  a=allocvecs(vi) fnq=new NQS("col","gid","id","type","ice","freq") snq=$o1 ce=$o2.ce
  if(numarg()>2) tt=$3 else tt=tstop
  snq.tog("DB") vi.indgen(0,$o2.allcells-1,1)
  for vtr(&i,vi) fnq.append(ce.o(i).col,ce.o(i).gid,i,ce.o(i).type,ice(ce.o(i).type),1e3*snq.select(-1,"id",i)/tt)
  dealloc(a) return fnq
}

//* clrMyNQs - delete nqCO[CDX], nqIN[CDX]
proc clrMyNQs () { local i
  nqsdel(nqCO[CDX])
  for i=0,9 nqsdel(nqIN[CDX][i])
  nqsdel(snq[CDX])
  nqsdel(nqisi[CDX])
  nqsdel(nqCTY[CDX])
}

//* rrounddI
func rrounddI () {
  return int($1+0.5)
}

//** getnqlfpcor(nqs with columns as LFP, window size in ms)
obfunc getnqlfpcor () { local winsz,tt,et,maxt,maxi,i,j,a localobj v0,v1,nqlfp,nqo
  nqlfp=$o1 a=allocvecs(v0,v1) winsz=$2
  {maxt=nqlfp.v.size*vdt_INTF6-1e3 maxi=nqlfp.v.size-1 nqo=new NQS("c0","c1","tt","r")}
  {nqo.v.resize((maxt/winsz)*nqlfp.m*(nqlfp.m-1)/2) nqo.pad() nqo.clear()}
  for(tt=1e3;tt+winsz<=maxt;tt+=winsz) {
    {et=MINxy((tt+winsz)/vdt_INTF6,maxi) vrsz(0,v0,v1)}
    for i=0,nqlfp.m-1 { v0.copy(nqlfp.v[i],tt/vdt_INTF6,et)
      for j=i+1,nqlfp.m-1 { v1.copy(nqlfp.v[j],tt/vdt_INTF6,et)
        nqo.append(i,j,tt,v0.pcorrel(v1))
      }
    }
  }
  {dealloc(a) return nqo}
}

//** getlfppco(nqs from getnqlfpcor)
obfunc getlfppco () { local t1,t2,csz,rc,i,j,a localobj mc,v0,v1,nqc
  {a=allocvecs(v0,v1) nqc=$o1 nqc.tog("DB") nqc.verbose=0}
  {v0.resize(nqc.v.size) nqc.getcol("c0").uniq(v0)}
  {csz=(1+v0.size)*(v0.size)/2 rc=nqc.fi("r") vrsz(0,v0,v1)}
  {mc=new Matrix(nqc.v.size/csz,nqc.v.size/csz) i=j=0}
  for(t1=0;t1<nqc.v.size;t1+=csz) { v0.copy(nqc.v[rc],t1,t1+csz-1)
    {mc.x(i,i)=1 j=i+1}
    for(t2=t1+csz;t2<nqc.v.size;t2+=csz) { v1.copy(nqc.v[rc],t2,t2+csz-1)
      mc.x(i,j)=mc.x(j,i)=v0.pcorrel(v1)
      j+=1
    }
    i+=1
  }
  {nqc.verbose=1 dealloc(a) return mc}
}

//** gotolfp(nqlfp,nqspud,tmin,tmax,offy)
func gotolfp () { local a,b,c,i,offy,tmin,tmax localobj nqlfp,nqspud,vtmpx,vtmpy,vsub
  if(glfp==nil)glfp=new Graph()
  a=allocvecs(vtmpx,vtmpy,vsub)
  {nqlfp=$o1 nqspud=$o2 tmin=$3 tmax=$4 offy=$5}
  if(vLFP==nil) {
    for i=0,nqlfp.m-1 {
      vLFP[i]=new Vector()
      vLFP[i].plot(glfp,vdt_INTF6)
    }
  }
  glfp.erase
  vsub.resize(nqlfp.m)
  for i=0,nqlfp.m-1 { // draw lfps
    {vLFP[i].resize(0) vLFP[i].copy(nqlfp.v[i],tmin/vdt_INTF6,tmax/vdt_INTF6)}
    {vsub.x(i)=vLFP[i].min vLFP[i].add(offy*i-vsub.x(i)) vLFP[i].line(glfp,vdt_INTF6)}
  }
  if(nqspud!=nil) { // draw peak of bumps
    {nqspud.verbose=0 vrsz(0,vspudx,vspudy)}
    for i=0,nqlfp.m-1 if(nqspud.select("LOC","[]",tmin,tmax,"COL",i)) {
      vrsz(0,vtmpx,vtmpy)
      {vtmpy.copy(nqspud.getcol("PEAK")) vtmpy.add(i*offy-nqlfp.v[i].min)}
      {vtmpx.copy(nqspud.getcol("LOC")) vtmpx.sub(tmin)}
      {vspudx.append(vtmpx) vspudy.append(vtmpy)}
    }
    vspudy.mark(glfp,vspudx,"O",4,2,1)
    {nqspud.verbose=1 dealloc(a)}
  }
  {dealloc(a) glfp.exec_menu("View = plot")}
  return 1
}

//** make a single NQS out of nqLFP, using last full LFP in each
obfunc getnqlfp () { local i localobj nqlfp 
  nqlfp=new NQS(numcols)
  for i=0,numcols-1 nqlfp.v[i].copy(nqLFP[i].v[nqLFP[i].m-1])
  return nqlfp
}

//** getnqspud(nqs with columns as LFP)
obfunc getnqspud () { local i,off,x,a localobj vec,nqu,nqin,nqtmp
  {a=allocvecs(vec) test_do_graphs=0 nqin=$o1 printStep=vdt_INTF6}
  if(numarg()>1) off=$2 else off=1e3
  for i=0,nqin.m-1 {
    vec.resize(0) vec.copy(nqin.v[i],off/vdt_INTF6,(tstop-off)/vdt_INTF6)
    nqtmp=testupdown(vec,10,0)
    {nqtmp.resize("COL") nqtmp.pad() nqtmp.v[nqtmp.m-1].fill(i)}
    if(nqu==nil) {
      {nqu=new NQS() nqu.cp(nqtmp)}
    } else nqu.append(nqtmp)
    nqsdel(nqtmp)
  }  
  for case(&x,0,5) nqu.v[x].add(off) // add offset so have sim times
  {dealloc(a) return nqu}
}

//* initnqCTY(col id) - makes an NQS with summed spike-train activity by type, requires nqCO to be setup already
// nqCO[cid].v[celltype] is summed activity of the given celltype
proc initnqCTY () { local i,cid,ty
  {cid=$1 nqsdel(nqCTY[cid]) nqCTY[cid]=new NQS(CTYPi)}
  for i=0,nqCO[cid].m-1 {
    ty=col[cid].ce.o(i).type
    if(nqCTY[cid].v[ty].size==0) nqCTY[cid].v[ty].resize(nqCO[cid].v[i].size)
    nqCTY[cid].v[ty].add(nqCO[cid].v[i])
  }
}

//* getnqCTYnt - returns an nqs of nTE between summed counts per type (nqCTY[...]), across/within columns
obfunc getnqCTYnt () { local i,c1,c2,nt,nshuf,from,to,a localobj nq,vtmp
  a=allocvecs(vtmp) vtmp.resize(35) nshuf=30
  {nq=new NQS("col1","col2","fs","ts","from","to","ic1","ic2","nte","samec") nq.strdec("fs") nq.strdec("ts")}
  for c1=0,numcols-1 for c2=0,numcols-1 {
    for from=0,CTYPi-1 for to=0,CTYPi-1 {
      if(from==to && c1==c2) continue
      if(nqCTY[c1].v[from].size<1 || nqCTY[c2].v[to].size<1) continue
      nt = nqCTY[c1].v[from].tentropspks(nqCTY[c2].v[to],vtmp,nshuf)
      nq.append(c1,c2,CTYP.o(from).s,CTYP.o(to).s,from,to,ice(from),ice(to),nt,c1==c2)
    }
  }
  {dealloc(a) return nq}
}

//* getnTEmap(nqs from getnqCTYnt, col1, col2, type1, type2, map of types to matrix coordinates)
// returns a matrix where element x,y is nTE from row y -> column x
obfunc getnTEmap () { local ty1,ty2,c1,c2,cid,from,to localobj mc,vty1,vty2,vm,nq
  nq=$o1 c1=$2 c2=$3 vty1=$o4 vty2=$o5 vm=$o6 nq.verbose=0
  mc=new Matrix(vty1.size,vty2.size)
  for vtr(&from,vty1) for vtr(&to,vty2) if(nq.select("col1",c1,"col2",c2,"from",from,"to",to)) {    
    mc.x(vm.x(from),vm.x(to)) = MAXxy(nq.fetch("nte"),0)
  } else mc.x(vm.x(from),vm.x(to))=0
  nq.verbose=1
  return mc
}

//* colnTEmap(nq from getnqCTYnt, list for output matrices,inter flag)
// gets COLUMN nTE map, uses all cell types and columns
// returns average of list for output matrices
// iff inter == 0, then gets INTRA-column map, iff inter == 1, gets INTER-column map
obfunc colnTEmap () { local i,j,cnt,inter,a localobj mc,vty1,vty2,vm,lc,nq
  {a=allocvecs(vty1,vty2,vm) vrsz(0,vty1,vty2)}
  nq=$o1 lc=$o2 if(lc==nil){lc=new List()} inter=$3
  {vty1.append(E2,I2,I2L,E4,I4,I4L,E5R,E5B,I5,I5L,E6,I6,I6L) vty2.copy(vty1)}
  {vm.resize(CTYPi) vm.fill(0) cnt=0}
  for vtr(&i,vty1) {vm.x(i)=cnt cnt+=1} // map cell types to index in Matrix
  if(inter) {
    for i=0,numcols-1 for j=0,numcols-1 if(i!=j) lc.append(getnTEmap(nq,i,j,vty1,vty2,vm))
  } else for i=0,numcols-1 lc.append(getnTEmap(nq,i,i,vty1,vty2,vm))
  mc=avgmat(lc) // average Matrix
  {dealloc(a) return mc}
}

//* initMyNQs - init binned value NQSs using raster info from vit[CDX] and input info from vq
// nqCO[CDX] = nqs with counts/time for each cell
// nqIN[CDX][sy] = nqs with counts/time for external inputs to synapse type sy @ all the cells
proc initMyNQs () { local ct,idx,sy,a localobj vidx,vt
  a=allocvecs(vidx,vt) vq[CDX].verbose=0
  if(nqCO[CDX]==nil) nqCO[CDX]=new NQS(col[CDX].allcells) else nqCO[CDX].clear() //NQS for the binned spike-trains
  vit[CDX].vec.ihist(vit[CDX].tvec,nqCO[CDX].vl,htmin,htmax,binsz)
  for case(&sy,AM,NM,GA,AM2,NM2,GA2) if(vq[CDX].select("sy",sy)) {
    flag_stats=1 //allow vq[CDX].v[1] to be unsorted  
    if(nqIN[CDX][sy]==nil) nqIN[CDX][sy]=new NQS(vq[CDX].v[0].max+1) else nqIN[CDX][sy].clear() //NQS for the binned input events
    {vidx.copy(vq[CDX].getcol("ind")) vt.copy(vq[CDX].getcol("time"))}
    vidx.ihist(vt,nqIN[CDX][sy].vl,htmin,htmax,binsz)  
  }
  if(!skipsnq) {
    {nqsdel(snq[CDX]) snq[CDX]=SpikeNQS(vit[CDX].tvec,vit[CDX].vec,lcol.o(CDX))}
    {nqsdel(fnq[CDX]) fnq[CDX]=FreqNQS(snq[CDX],lcol.o(CDX))}
  }
  initnqCTY(CDX)
  {vq[CDX].verbose=1 dealloc(a)}
}

//* initAllMyNQs() - initMyNQs for all numcols
proc initAllMyNQs () { local cdx
  cdx=CDX
  for CDX=0,numcols-1 initMyNQs()
  CDX=cdx
}

//* clrAllMyNQs() - clrMyNQs for all numcols
proc clrAllMyNQs () { local cdx
  cdx=CDX
  for CDX=0,numcols-1 clrMyNQs()
  CDX=cdx
}
 
//* getTEKTaunq() gets transfer entropy and ktau from inputs->cells
obfunc getTEKTaunq () { local sy,loc,idx,vbt localobj nqt,vtmp
  {nqt=new NQS("id","sy","loc","te","nte","ktau") initMyNQs()}
  {vbt=verbose_infot verbose_infot=0}
  for case(&sy,AM,NM,GA,AM2,NM2,GA2) if(nqIN[CDX][sy]!=nil) {
    if(sy<AM2)loc=SOMA else loc=DEND
    for idx=0,col[CDX].allcells-1 { vtmp=normte(nqIN[CDX][sy].v[idx],nqCO[CDX].v[idx],30)
      if(idx%100==0) printf("%s:idx=%d\n",STYP.o(sy).s,idx)
      nqt.append(idx,sy,loc,vtmp.x(0),vtmp.x(2),nqIN[CDX][sy].v[idx].kcorrel(nqCO[CDX].v[idx],1))
    }
  }
  for idx=nqt.m-2,nqt.m-1 nqt.v[idx].unnan() //get rid of nans
  {verbose_infot=vbt return nqt}
}
 
// //* getTEnqAllsub -- 
// // gets nTE from all inputs of a cell type -> all outputs of a specific cell type
// obfunc getTEnqAllsub () { local a,row,ct,sy,id,poid,pls,prid,cts,tyt,\
//                        idcol,sycol,INTFidcol,pulsecol,loccol\
//                        localobj nqt,vtmp
// cts=$1
// nqin.tog("DB")  nqt=new NQS("id","type","prid","poid","pulse","sy","same","loc","te","nte")
// initMyNQs()
// {idcol=nqin.fi("id")  sycol=nqin.fi("sy")  INTFidcol=nqin.fi("INTFid")}
// {pulsecol=nqin.fi("pulse")  loccol=nqin.fi("loc")}
// {nqin.select(-1,"type",cts)  nqin.verbose=0 pls=0}
// printf("\nid:")
// for id=ix[cts],ixe[cts] {
//   ct = ce.o(id).type
//   if(id%20==0) printf("\n")
//   printf("%d ",id)
//   for vtr(&row,nqin.ind) {
//     {prid=nqin.v[INTFidcol].x(row)      poid=nqin.v[idcol].x(row)}
//     {sy=nqin.v[sycol].x(row)      loc=nqin.v[loccol].x(row)}
//     vtmp = normte(nqIN[CDX].v[prid],nqCO[CDX].v[id],30)
//     nqt.append(id,ct,prid,poid,pls,sy,poid==id,loc,vtmp.x(0),vtmp.x(2))
//   }
// }
// printf("\n")
// for id=0,nqt.m-1 nqt.v[id].unnan()
// nqin.verbose=1
// return nqt
// }
// 
// 
// //* ccntevsbinsz(prid,poid,vbins,nq)
// obfunc ccntevsbinsz () { local prid,poid localobj vbins,vtmp,nq
// prid=$1 poid=$2 vbins=$o3
// nq=$o4
// if(nq==nil) nq=new NQS("prid","poid","binsz","te","nte","prty","poty") else nq.tog("DB")
// for vtr(&binsz,vbins) {
//   print "binsz" , binsz
//   clrMyNQs()
//   initMyNQs()
//   vtmp = normte(nqCO[CDX].v[prid],nqCO[CDX].v[poid],30)
//   nq.append(prid,poid,binsz,vtmp.x(0),vtmp.x(2),ce.o(prid).type,ce.o(poid).type)
// }
// return nq
// }
// //* tyntevsbinsz(Vec of pre-type, Vec of po-type, Vec of bin sizes, nq output)
// obfunc tyntevsbinsz () { local prid,poid,prty,poty,verb localobj vbins,vtmp,nq,vprty,vpoty
// vprty=$o1 vpoty=$o2 vbins=$o3 nq=$o4
// if(numarg()>4)verb=$5 else verb=0
// if(nq==nil) nq=new NQS("prid","poid","binsz","te","nte","prty","poty") else nq.tog("DB")
// nclnq[CDX].verbose=0
// for vtr(&binsz,vbins) {
//   print "binsz " , binsz
//   clrMyNQs()
//   initMyNQs()
//   for vtr(&prty,vprty) for vtr(&poty,vpoty) {
//     if(prty==poty || div[prty][poty][0]<1) continue
//     print "prty : " , CTYP.o(prty).s, " , poty : ", CTYP.o(poty).s
//     for prid=ix[prty],ixe[prty] {
//       if(nclnq[CDX].select("prid",prid,"poty",poty)) for nclnq[CDX].qt(&poid,"poid") {
//         if(verb) print "prid " , prid , " poid " , poid
//         vtmp=normte(nqCO[CDX].v[prid],nqCO[CDX].v[poid],30)
//         nq.append(prid,poid,binsz,vtmp.x(0),vtmp.x(2),prty,poty)
//       }
//     }
//   }
// }
// nclnq[CDX].verbose=1
// return nq
// }

//* binarize(vec) - set all values > 0 to 1, everything else to 0
proc binarize () { local i,sz localobj vin
  vin=$o1 sz=vin.size()-1
  for i=0,sz if(vin.x(i)>0) vin.x(i)=1 else vin.x(i)=0
}

//* mknclnq - make an NQS with connectivity info
proc mknclnq () { local idx,jdx,ty1,ty2,sy,loc,ic1,ic2,lt1,gid1,gid2,c1,c2,c localobj adj,ce,nc
  if(nclnq!=nil) return
  nclnq=new NQS("prid","poid","prty","poty","sy","loc","col1","col2","gid1","gid2","w")
  for c=0,numcols-1 { // intra-column wiring for all the columns
    {ce=col[c].ce adj=AdjList(0,ce.count-1,0,ce)}
    for idx=0,adj.count-1 { ty1=ce.o(idx).type
      gid1=ce.o(idx).gid
      if(!ice(ty1)) {
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
          nclnq.append(idx,jdx,ty1,ty2,AM2,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][AM2][0])
          nclnq.append(idx,jdx,ty1,ty2,NM2,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][NM2][0])
        }
      } else {
        if(IsLTS(ty1)) sy=GA2 else sy=GA
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
          nclnq.append(idx,jdx,ty1,ty2,sy,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][sy][0])
        }
      }
    }
  }
  for ltr(nc,ncl) { // inter-column wiring for all the columns
    if(!(isojt(nc.pre,INTF6[0]) && isojt(nc.syn,INTF6[0]))) continue // only looking btwn INTF6 cells
    {idx=nc.pre.id ty1=nc.pre.type gid1=nc.pre.gid c1=nc.pre.col}
    {jdx=nc.syn.id ty2=nc.syn.type gid2=nc.syn.gid c2=nc.syn.col}
    if(!ice(ty1)) {
      nclnq.append(idx,jdx,ty1,ty2,AM2,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(AM2))
      nclnq.append(idx,jdx,ty1,ty2,NM2,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(NM2))
    } else {
      if(IsLTS(ty1)) sy=GA2 else sy=GA
      nclnq.append(idx,jdx,ty1,ty2,sy,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(sy))
    }
  }
}

//* FullAdjList(makenclnq,skipI) - make adjacency list that has intercol/intracol info
// uses global identifier (gid) of each cell (not intracol id)
obfunc FullAdjList () { local i,j,sz,skipI,mknq,a localobj adj,vt,vg1,vg2,vex
  if(nclnq==nil)mknq=1 else {
    if(numarg()>0)mknq=$1 else mknq=1
  }
  if(mknq) {nqsdel(nclnq) mknclnq()} // has all the connectivity info
  nclnq.tog("DB")
  if(numarg()>1)skipI=$2 else skipI=1 //skip inhib cells by default
  a=allocvecs(vex,vt)
  {adj=new List() vg1=nclnq.getcol("gid1") vg2=nclnq.getcol("gid2")}
  sz=MAXxy(vg1.max,vg2.max)
  for i=0,sz adj.append(new Vector())
  if(skipI) { vex.resize(sz+1) // skip inhib cells
    for i=0,sz vex.x(i)=(INTF6[i].flag("inh")==0)
    for i=0,vg1.size-1 if(vex.x(vg1.x(i)) && vex.x(vg2.x(i))) {
      adj.o(vg1.x(i)).append(vg2.x(i))
    }
  } else for i=0,vg1.size-1 adj.o(vg1.x(i)).append(vg2.x(i))
  for i=0,sz if(adj.o(i).size) { // make sure all are unique, since nclnq has connection info for each synapse type
    {vt.resize(adj.o(i).size) adj.o(i).uniq(vt)}
    {adj.o(i).resize(0) vt.sort() adj.o(i).copy(vt)}
  }
  dealloc(a)
  return adj
}

//* getsy(from,to) - get synapse type , either returns AM2 or GA,GA2 , never returns NMDA
func getsy () { local from,to
  from=$1 to=$2
  if(ice(from)) {
    if(IsLTS(from)) return GA2 else return GA
  } else return AM2
}
 
//* getvpow(nqf,minf,maxf)
//get the power in a range of frequencies from nqf, the matspecgram spectrogram nqs from nqfld
obfunc getvpow () { local a,minf,maxf,i localobj vtmp,vout,nqf,vf
  nqf=$o1 minf=$2 maxf=$3
  a=allocvecs(vtmp,vf) vout=new Vector()
  {vtmp=nqf.get("pow",0).o print vtmp vout.resize(vtmp.size) vout.fill(0)}
  {vf.resize(nqf.v.size)  nqf.getcol("f").uniq(vf)}
  for i=0,vf.size-1 {
    if(vf.x(i)>=minf && vf.x(i)<=maxf) {
      vtmp=nqf.get("pow",i).o
      vout.add(vtmp)
    }
    if(vf.x(i)>maxf) break
  }
  dealloc(a)
  return vout
}

//* getpfcormat(nq1[,nq2]) -- nq1,nq2 are from matspecgram, getpfcormat gets power fluctuation correlation matrix
// can get nq1,nq2 from nqfld+matspecgram like this: nq1=matspecgram(nqLFP.v,sampr,100,0,1)
// nq2 is optional
obfunc getpfcormat () { local r1,r2,sz,a localobj mc,v1,v2,nq1,nq2
  a=allocvecs(v1,v2) 
  if(numarg()==1){nq1=nq2=$o1} else {nq1=$o1 nq2=$o2}
  {sz=MAXxy(nq1.v.size,nq2.v.size) mc=new Matrix(sz,sz)}
  if(nq1==nq2) {
    for r1=0,nq1.v.size-1 { v1.copy(nq1.get("pow",r1).o)
      mc.x(r1,r1)=1
      for r2=r1+1,nq2.v.size-1 { v2.copy(nq2.get("pow",r2).o)
        mc.x(r1,r2) = mc.x(r2,r1) = v1.pcorrel(v2)
      }
    }
  } else {
    for r1=0,nq1.v.size-1 { v1.copy(nq1.get("pow",r1).o)
      for r2=0,nq2.v.size-1 { v2.copy(nq2.get("pow",r2).o)
        mc.x(r1,r2) = v1.pcorrel(v2)
      }
    }
  }
  dealloc(a)
  return mc
}

//* getsyncarea(nqf,mc,minf,maxf,[th]) -- nqf is from matspecgram, mc is from getpfcormat
// gets the 'synchronization area' -- the ratio of pixels in matrix with r > th / that area of the matrix
// frequencies/area used are in range [minf,maxf]
func getsyncarea () { local r1,c1,sa,cnt,minf,maxf,th localobj nqf,mc
  nqf=$o1 mc=$o2 minf=$3 maxf=$4 sa=cnt=0
  if(numarg()>4)th=$5 else th=0.2
  for r1=0,mc.nrow-1 if(nqf.v[0].x(r1)>=minf && nqf.v[0].x(r1)<=maxf) {
    for c1=0,r1-1 if(nqf.v[0].x(c1)>=minf && nqf.v[0].x(c1)<=maxf) {
      if(mc.x(r1,c1)>th) sa+=1
      cnt+=1
    }
  }
  if(cnt>0) return sa/cnt else return 0
}
 
//* FileExists(path) - tests whether file exists by opening in 'r' mode
func FileExists(){ localobj f
  f = new File()
  f.ropen($s1)
  if(f.isopen()){
    f.close()
    return 1
  }
  return 0
}

//* loadmydat() - load in data from sim - this func depends on strv being set
// $s1 == vecs file, [$s2 == nclnq file]
func loadmydat () { local cdx,idx,a localobj co,vi,str,vtmp,svq
  {a=allocvecs(vtmp) strvecs=$s1 str=new String2() svq=new String()}
  clrAllMyNQs() // free histogram NQSs
  printf("reading in %s vecs file\n",strvecs)
  gvnew(strvecs) // load the Vectors
  for ltr(co,lcol,&cdx) {
    {vit[cdx].vec.resize(0) vit[cdx].tvec.resize(0)}
    {sprint(str.s,"%s_SPKS",co.name) sprint(str.t,"%s_LFP",co.name)}
    for ltr(vi,panobj.printlist,&idx) {
      if(!strcmp(vi.name,str.s)) {
        print "reading in spikes for ",co.name, " from vitem ", vi.name
        panobj.rv_readvec(idx,vit[cdx].tvec,vit[cdx].vec) //read in raster
      }
    }
    sprint(svq.s,"data/%s_%s_vq.nqs",strv,co.name)
    if(FileExists(svq.s)) {
      print "reading in vq for ", co.name, " from ", svq.s
      vq[cdx].rd(svq.s)
    } else print "loadmydat WARNING: file ", svq.s , " doesn't exist."
    sprint(str.s,"data/%s_%s_LFP.nqs",strv,col[cdx].name)
    if(FileExists(str.s)) {
      {nqsdel(nqLFP[cdx]) nqLFP[cdx]=new NQS(str.s)}
    } else print "loadmydat WARNING: file ", str.s , " doesn't exist."
  }
  if(numarg()>1) {
    {strnclnq=$s2 printf("reading in %s nclnq nqs file\n",strnclnq)}
    if(nclnq==nil) nclnq=new NQS(strnclnq) else nclnq.rd(strnclnq)
  }
  {dealloc(a) return 1}
}
 
//* pravgrates([duration,allcolumns,strdef output]) - print average rates of cell types after run
//print average firing rates from vit[...], assumes sgrdur,numc same as when ran sim
//when allcolumns==0 (default) uses vit[CDX] (CDX is global index for COLUMN)
proc pravgrates () { local idx,ct,dur,cdx,sidx,eidx,all,a localobj vspk,vnc,str
  if(numarg()>0)dur=$1 else dur=tstop
  if(numarg()>1)all=$2 else all=0   
  {a=allocvecs(vspk,vnc) vrsz(CTYPi+1,vspk,vnc) idx=ct=0 sidx=eidx=CDX str=new String2()}
  if(all){sidx=0 eidx=numcols-1}
  for cdx=sidx,eidx {
    for vtr(&idx,vit[cdx].vec) vspk.x(col[cdx].ce.o(idx).type)+=1
    for ct=0,CTYPi-1 vnc.x(ct)+=col[cdx].numc[ct]
  }
  for ct=0,CTYPi-1 if(vnc.x(ct)) {
    sprint(str.t,"avg rate for %s: %gHz\n",CTYP.o(ct).s,(1e3*vspk.x(ct)/dur)/vnc.x(ct))
    strcat(str.s,str.t)
  }
  print str.s
  if(numarg()>2) $s3=str.s
}
//* prctcellsfired - get % of cells that fired after a run
func prctcellsfired () { local a,i,pr,ct localobj vf
  a=allocvecs(vf) 
  if(numarg()>0) {
    ct=$1
    if(ct<0 || ct>=CTYPi) {
      print "invalid ct = ", ct
      return 0
    }
    vrsz(numc[ct],vf) vf.fill(0)
    for vtr(&i,vit[CDX].vec) if(col[CDX].ce.o(i).type==ct) vf.x(i-ix[ct])=1
    pr=vf.count(1)/numc[ct]
  } else {
    vrsz(col[CDX].allcells,vf) vf.fill(0)
    for vtr(&i,vit[CDX].vec) vf.x(i)=1  
    pr=vf.count(1)/col[CDX].allcells
  }
  dealloc(a)
  return pr
}
 
//* minrunsv() minimal run/save of data , just save vecs, vq, nclnq
// $1 - optional , specifices whether to save nclnq
proc minrunsv () { local svc,cdx,svq localobj st,co
  {st=new String() time("run()")}
  printf("\nran sim for %g\n",tstop)
  sprint(st.s,"data/%s_.vecs",strv)
  {panobj.pvplist(st.s,strv) printf("saved vecs\n")}
  for cdx=0,numcols-1 { // save NQS with LFPs
    sprint(tstr,"data/%s_%s_LFP.nqs",strv,col[cdx].name)
    {nqLFP[cdx].tog("DB") nqLFP[cdx].sv(tstr)}
  }
  if(numarg()>0) svc=$1 else svc=1
  if(numarg()>1) svq=$2 else svq=1
  if(svc) {
    sprint(st.s,"data/%s_nclnq.nqs",strv)
    {mknclnq() nclnq.tog("DB") nclnq.sv(st.s)}
    printf("saved connectivity\n")
  }
  if(svq) for ltr(co,lcol,&cdx) {
    sprint(st.s,"data/%s_%s_vq.nqs",strv,co.name)
    {vq[cdx].tog("DB") vq[cdx].sv(st.s) printf("saved inputs for %s\n",co.name)}
  }
}

//$s1=version string, loads the data created from minrunsv - should set strv first
proc loadminrundat () { localobj st
  {st=new String2() sprint(st.s,"data/%s_.vecs",$s1) sprint(st.t,"data/%s_nclnq.nqs",$s1)}
  if(FileExists(st.t)) loadmydat(st.s,st.t) else loadmydat(st.s)
}