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

 Download zip file   Auto-launch 
Help downloading and running models
This model is an extension of a model (<a href="">138379</a>) recently published in Frontiers in Computational Neuroscience. This model consists of 4700 event-driven, rule-based neurons, wired according to anatomical data, and driven by both white-noise synaptic inputs and a sensory signal recorded from a rat thalamus. Its purpose is to explore the effects of cortical damage, along with the repair of this damage via a neuroprosthesis.
1 . 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 [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 V1 pyramidal intratelencephalic L2-5 cell; Neocortex V1 interneuron basket PV cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell;
Channel(s): I Chloride; I Sodium; I Potassium;
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA; Gaba;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Deep brain stimulation; Information transfer; Brain Rhythms;
Implementer(s): Lytton, William [billl at]; Neymotin, Sam [samn at]; Kerr, Cliff [cliffk at];
Search NeuronDB for information about:  Neocortex V1 pyramidal corticothalamic L6 cell; Neocortex V1 pyramidal intratelencephalic L2-5 cell; Neocortex V1 interneuron basket PV cell; GabaA; AMPA; NMDA; Gaba; I Chloride; I Sodium; I Potassium; Gaba; Glutamate;
infot.mod *
intf6_.mod *
intfsw.mod *
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
decmat.hoc *
decnqs.hoc *
default.hoc *
drline.hoc *
infot.hoc *
local.hoc *
misc.h *
ratlfp.dat *
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
xgetargs.hoc *
// $Id: nload.hoc,v 1.182 2010/09/20 00:45:08 samn Exp $ 

print "Loading nload.hoc..."

strdef strvecs,strnclnq //vecs file, vq nqs file, nqin nqs file, nqfld nqs file
strdef strv
{declare("vspudx",new Vector(),"vspudy",new Vector())}
if(name_declared("snq")) {
  objref snq[numcols]
} else declare("snq",tstr)
if(name_declared("fnq")) {
  objref fnq[numcols]
} else declare("fnq",tstr)
if(name_declared("vq")) {
  objref vq[numcols]
} else declare("vq",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

//* Deriv() return derivative of input vector
//output vector has same size as input
//$o1 = input vector
//vp = output vector
obfunc Deriv(){ localobj vp
  vp=new Vector($o1.size)
  if($o1.size < 2){
    printf("Deriv ERRAA: input vec size < 2!\n")
    return vp
  return vp

//* IsiNQS()
obfunc IsiNQS () { local idx,tt,jdx,a localobj vtmp,nqisi,vid
  nqisi=new NQS("id","isi")
  for idx=0,col[CDX].allcells-1 {
    vtmp = Deriv(snq[CDX].getcol("SpikeT"))
  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") 
  {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[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*,"id",i)/tt)
  dealloc(a) return fnq

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

//* 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)
  {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"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)
  {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()
  {nqlfp=$o1 nqspud=$o2 tmin=$3 tmax=$4 offy=$5}
  if(vLFP==nil) {
    for i=0,nqlfp.m-1 {
      vLFP[i]=new Vector()
  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("LOC","[]",tmin,tmax,"COL",i)) {
      {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)}
    {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.resize("COL") nqtmp.pad() nqtmp.v[nqtmp.m-1].fill(i)}
    if(nqu==nil) {
      {nqu=new NQS() nqu.cp(nqtmp)}
    } else nqu.append(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 {
    if(nqCTY[cid].v[ty].size==0) nqCTY[cid].v[ty].resize(nqCO[cid].v[i].size)

//* 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)
  {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("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
  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
  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"))}
  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))}
  {vq[CDX].verbose=1 dealloc(a)}

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

//* clrAllMyNQs() - clrMyNQs for all numcols
proc clrAllMyNQs () { local cdx
  for CDX=0,numcols-1 clrMyNQs()
//* 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)
  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()
// {"id")"sy")"INTFid")}
// {"pulse")"loc")}
// {,"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
      if(!ice(ty1)) {
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
      } else {
        if(IsLTS(ty1)) sy=GA2 else sy=GA
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
  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
    { ty1=nc.pre.type gid1=nc.pre.gid c1=nc.pre.col}
    { ty2=nc.syn.type gid2=nc.syn.gid c2=nc.syn.col}
    if(!ice(ty1)) {
    } else {
      if(IsLTS(ty1)) sy=GA2 else sy=GA

//* 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
  if(numarg()>1)skipI=$2 else skipI=1 //skip inhib cells by default
  {adj=new List() vg1=nclnq.getcol("gid1") vg2=nclnq.getcol("gid2")}
  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))) {
  } 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)}
  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

//* GraphNQ - get an NQS with graph properties (path length,clustering-coeff,centrality,etc.)
obfunc GraphNQ () { local i,j,c,ediv,econv,iconv,idiv,intraediv,intraeconv,intraidiv,intraiconv,\
                   localobj gnq,vexl,vexcc,vcent,cel,vc,vd,viec,vied,viid,viic,xo,adjf,vwexl,\
  {adjf=adj=FullAdjList() vexl=GetNetEXL() vexcc=GetNetEXCC()}
  {vcent.resize(adj.count) GetCentrality_intfsw(adj,vcent)}  
  gnq=new NQS("gid","id","ty","ice","col")
  gnq.resize("exl","excc","cent","ediv","econv","idiv","iconv") //pathlen,clustering,centrality,total div/conv
  gnq.resize("interEdiv","interEconv","interIdiv","interIconv") //intercolumn div/conv
  gnq.resize("intraEdiv","intraEconv","intraIdiv","intraIconv") //intracolumn div/conv
  gnq.resize("wediv","weconv","widiv","wiconv") //weighted div/conv
  gnq.resize("winterEdiv","winterEconv","winterIdiv","winterIconv") //weighted intercolumn div/conv
  gnq.resize("wintraEdiv","wintraEconv","wintraIdiv","wintraIconv") //weighted intracolumn div/conv
  {vrsz(CTYPi,vc,vd) vrsz(vcent.size,viec,vied,viid,viic,vwiec,vwied,vwiid,vwiic) vtmp.resize(adj.count)}
  for ltr(xo,ncl) if(isojt(xo.pre,INTF6[0]) && isojt(xo.syn,INTF6[0])) { // intercol div/conv
    if(ice(xo.pre.type)) {
    } else {
    if(ice(xo.syn.type)) {
    } else {
  for i=0,vcent.size-1 { cel=INTF6[i]
    {cel.getdvi(vtmp) vd.fill(0)}
    for j=0,vtmp.size-1 vd.x(col[cel.col].ce.o(vtmp.x(j)).type)+=1
    {cel.getconv(1.2,vc) ediv=econv=iconv=idiv=wediv=weconv=wiconv=widiv=0}
    for j=0,CTYPi-1 {
      if(ice(j)) {
        idiv+=vd.x(j) iconv+=vc.x(j)
        widiv += vd.x(j)*col[cel.col].wmat[cel.type][j][getsy(cel.type,j)]
        wiconv += vc.x(j)*col[cel.col].wmat[j][cel.type][getsy(j,cel.type)]
      } else {
        econv+=vc.x(j) ediv+=vd.x(j)
        wediv += vd.x(j)*col[cel.col].wmat[cel.type][j][getsy(cel.type,j)]
        weconv += vc.x(j)*col[cel.col].wmat[j][cel.type][getsy(j,cel.type)]
    {intraediv=ediv intraeconv=econv intraidiv=idiv intraiconv=iconv}
    {wintraediv=wediv wintraeconv=weconv wintraidiv=widiv wintraiconv=wiconv}
    {idiv+=viid.x(i) iconv+=viic.x(i) ediv+=vied.x(i) econv+=viec.x(i)} // add in intercol div/conv
    {widiv+=vwiid.x(i) wiconv+=vwiic.x(i) wediv+=vwied.x(i) weconv+=vwiec.x(i)}// add in intercol weighted div/conv
  {gnq.resize("intraexl","intraexcc","intracent") gnq.pad()} //intracolumn pathlength,clustering-coeff,centrality
  for c=0,numcols-1 { 
    {vexl=GetNetEXL() vexcc=GetNetEXCC()}
    {vcent.resize(adj.count) vcent.fill(0) GetCentrality_intfsw(adj,vcent)}
    for i=0,col[c].allcells-1 { gdx=col[c].ce.o(i).gid
  {adj=adjf dealloc(a) return gnq}
//* getnTEvecs(cell id, synapse, location, internal inputs vec, external inputs vec, vpl inputs vec,[list of internal input vecs])
// fills $o4,$o5,$o6 with corresponding inputs. returns 1 if VPL inputs were found, 0 otherwise
// last optional arg will be a list of internal input vecs
// can be used as a helper function in allntevsbinsz
// +++++++ NB: vvpl is not currently used since no VPL inputs implemented for this sim yet!! ++++++++
func getnTEvecs () { local cidx,sy,loc,hasvpl,prid,asl,a localobj vint,vext,vvpl,vprid,lint
  {cidx=$1 sy=$2 loc=$3 vint=$o4 vext=$o5 vvpl=$o6 a=allocvecs(vprid)}
  if(numarg()>6) {lint=$o7 asl=1 lint.remove_all()} else asl=0
  if(vext.size!=nqCO[CDX].v[cidx].size) vrsz(nqCO[CDX].v[cidx].size,vext,vint,vvpl)
  {vext.fill(0) vint.fill(0) vvpl.fill(0) hasvpl=0}
  // if(("id",cidx,"syn",sy,"loc",loc,"pulse",-1))>0) {          
  //   intfid=nqin.fetch("INTFid") // VPL input from an INTF
  //   vvpl.add(nqIN[CDX].v[intfid])
  // }
  vext.add(nqIN[CDX][sy].v[cidx]) // external inputs [but not VPL]  
  vprid.resize("poid",cidx,"sy",sy)) //select internal inputs to this cell
  if(vprid.size>0) { 
    nclnq.getcol("prid").uniq(vprid) // IDs to use
    if(asl) {
      for vtr(&prid,vprid) {
        lint.append(new Vector())
    } else {
      for vtr(&prid,vprid) vint.add(nqCO[CDX].v[prid])      
  return hasvpl
//* allntevsbinsz(Vec of bin sizes, nq output, verbose)
// columns of output NQS:
//  id = id of cell, type = type of cell, ice = whether inhibitory, binsz = bin sizes used in ms,
//  intnte = nTE from cells within network at the synapse, extnte = nTE from external inputs at the
//  synapse, sumnte = nTE from internal + external inputs at the synapse, sy = synapse type,
//  loc = location of synapse, vplnte = nTE from external VPL inputs at the synapse, ocnte = nTE from
//  inputs from other COLUMNs
obfunc allntevsbinsz () { local i,hasvpl,hasint,hasoc,sy,loc,prid,poid,prty,poty,verb,cidx,a,intfid,\
                         localobj vbins,vtmp,nq,vprty,vpoty,vint,vext,vvpl,vprid,vc,vsy,vloc,ce,voc
  vbins=$o1 nq=$o2
  if(numarg()>2)verb=$3 else verb=0
  if(nq==nil) nq=new NQS("gid","col","id","type","ice","binsz","intnte","extnte","sumnte","sy","loc","vplnte","ocnte") else nq.tog("DB")
  {mknclnq() nclnq.verbose=hasvpl=vplnte=0}
  {vsy.append(AM2,NM2,GA2,GA) vloc.append(DEND,DEND,DEND,SOMA)} // synapse/location pairs
  for vtr(&binsz,vbins) { // go thru different bin sizes
    print "binsz " , binsz    
    {clrAllMyNQs()    initAllMyNQs()} // initialize NQS with histograms
    for CDX=0,numcols-1 { print "COL ", CDX // go thru all COLUMNs
      {vrsz(nqCO[CDX].v.size,vext,vint,vvpl,voc) ce=col[CDX].ce}
      for cidx=0,col[CDX].allcells-1 {     //go thru all cells
        {ty=ce.o(cidx).type ic=ice(ty)} // type of cell & whether inhibitory
        if(cidx%100==0) print cidx
        for i=0,vsy.size-1 { // go thru all synapse/location pairs
          {sy=vsy.x(i) loc=vloc.x(i)} // synapse/location pair
          {vext.fill(0) vint.fill(0) vvpl.fill(0) voc.fill(0)} // initialize to 0
          //vplnte = -2 // invalid value
          //// nTE @ AMPA on DEND from VPL inputs    
          //if(("id",cidx,"syn",sy,"loc",loc,"pulse",-1))>0) {          
          //  intfid=nqin.fetch("INTFid") // VPL input from an INTF
          //  vvpl.add(nqIN[CDX].v[intfid])
          //  vtmp=normte(vvpl,nqCO[CDX].v[cidx],30) // calculate the nTE
          //  vplnte=vtmp.x(2) // external nTE from VPL        

          vext.add(nqIN[CDX][sy].v[cidx]) // external inputs [but not VPL]
          vtmp=normte(vext,nqCO[CDX].v[cidx],30) // calculate the nTE
          extnte=vtmp.x(2) // external nTE

          vint.fill(0) // sum of internal INTRA-COLUMN inputs to this cell synapse
          vprid.resize("poid",cidx,"sy",sy,"col1",CDX,"col2",CDX)) //select intra-COLUMN inputs to this cell
          if((hasint=vprid.size)>0) { 
            nclnq.getcol("prid").uniq(vprid) // IDs to use
            for vtr(&prid,vprid) vint.add(nqCO[CDX].v[prid])      
            intnte=vtmp.x(2) // internal nTE (cells within COLUMN)
            vext.add(vint) // external + internal inputs

          voc.fill(0) // sum of INTER-COLUMN inputs to this cell synapse
          vprid.resize("gid2",ce.o(cidx).gid,"sy",sy,"col1","!=",CDX)) //select inter-COLUMN inputs to this cell
          if((hasoc=vprid.size)>0) {
            nclnq.getcol("gid1").uniq(vprid) // GIDs to use
            for vtr(&prid,vprid) voc.add(nqCO[INTF6[prid].col].v[INTF6[prid].id])
            ocnte=vtmp.x(2) // internal nTE (cells within COLUMN)
            vext.add(voc) // external + internal + from other COLUMNs inputs

          //if(hasvpl) vext.add(vvpl) // add VPL inputs, if available
          if(hasoc || hasint) {
            vtmp=normte(vext,nqCO[CDX].v[cidx],30) // nTE of external+internal+other COLUMN inputs to cell's spikes
          } else sumnte=extnte
          nq.append(ce.o(cidx).gid,CDX,cidx,ty,ic,binsz,intnte,extnte,sumnte,sy,loc,vplnte,ocnte) // save it
  {nclnq.verbose=1 dealloc(a)}
  return nq

//* getpopcktaunq() gets correlations between pairs of cells - with sliding time window - used
// for population coordination matrix
// $1=winsz,$2=wininc
obfunc getpopcktaunq () { local a,i,j,kt,ty1,ty2,ic1,winsz,wininc,tt,mx localobj nqcc,v1,v2,vc
  {winsz=$1 wininc=$2}
  nqcc=new NQS("id1","id2","ty1","ty2","ic1","ic2","ktau","t")
  if(nqCO[CDX]==nil) initMyNQs()
  for(tt=0;tt<nqCO[CDX].v.size;tt+=wininc) {
    for i=0,nqCO[CDX].m-1 {
      if(i%100==0) printf("i=%d\n",i)
      for j=i+1,nqCO[CDX].m-1 {
  return nqcc
//* getpopcmat(nqs from getpopcktaunq, skipice) - get population coordination matrix
//$o1=nqs from getpopcktaunq
obfunc getpopcmat () { local a,tt,idx,jdx,skipice localobj mc,nq,vt,ls
  nq=$o1 nq.tog("DB") if(numarg()>1)skipice=$2 else skipice=0
  ls=new List()
  if(skipice) {
    for vtr(&tt,vt) {"t",tt,"ic1",0,"ic2",0)
      ls.append(new Vector())
  } else {
    for vtr(&tt,vt) {"t",tt)
      ls.append(new Vector())
  mc=new Matrix(ls.count,ls.count)
  for idx=0,ls.count-1 for jdx=idx+1,ls.count-1 mc.x(idx,jdx)=mc.x(jdx,idx)=ls.o(idx).pcorrel(ls.o(jdx))
  for idx=0,ls.count-1 mc.x(idx,idx)=1
  return mc
//* getcellcelltaunq([vc1,vc2]) gets correlations between pairs of cells
// 2 optional args are vectors of pairs of COLUMN IDs for which to get correlations
obfunc getcellcelltaunq () { local a,i,j,ii,kt,ty1,ty2,ic1,rel,c1,c2,gid1,gid2 localobj nqcc,vpv,ce1,ce2,vc1,vc2,nqtmp
  a=allocvecs(vpv,vc1,vc2) vpv.resize(1)
  nqcc=new NQS("id1","id2","ty1","ty2","ic1","ic2","ktau","pval","rel","col1","col2","gid1","gid2")
  {clrAllMyNQs()    initAllMyNQs()}
  if(numarg()>=2) {
    vc1.copy($o1) vc2.copy($o2)
  } else {
    for c1=0,numcols-1 for c2=0,numcols-1 {
      if(c2<c1) continue
      {vc1.append(c1) vc2.append(c2)}
  for ii=0,vc1.size-1 { 
    {c1=vc1.x(ii) ce1=col[c1].ce}
    {c2=vc2.x(ii) ce2=col[c2].ce}
    print "c",c1," c",c2
    nqtmp=new NQS("id1","id2","ty1","ty2","ic1","ic2","ktau","pval","rel","col1","col2","gid1","gid2")
    if(c1==c2) { // INTRA-COLUMN
      for i=0,nqCO[c1].m-1 {
        if(i%100==0) printf("i=%d\n",i)
        {ty1=ce1.o(i).type ic1=ice(ty1) gid1=ce1.o(i).gid}
        for j=i+1,nqCO[c2].m-1 {
          rel=0 //independent
          if(vpv.x(0) < 0.05) {
            if(kt>0) {
              rel=1 //positive
            } else if(kt<0) {
              rel=-1 //negative
    } else { // INTER-COLUMN
      for i=0,nqCO[c1].m-1 {
        if(i%100==0) printf("i=%d\n",i)
        {ty1=ce1.o(i).type ic1=ice(ty1) gid1=ce1.o(i).gid}
        for j=0,nqCO[c2].m-1 {
          rel=0 //independent
          if(vpv.x(0) < 0.05) {
            if(kt>0) {
              rel=1 //positive
            } else if(kt<0) {
              rel=-1 //negative
    {nqtmp.v["ktau")].unnan() nqcc.append(nqtmp) nqsdel(nqtmp)}
  {dealloc(a) return nqcc}
//* getcellcelltenq() gets transfer entropies between pairs
// of cells (in both directions)
// iff monosynapticflag == 1 only gets for cell pairs that are connected
obfunc getcellcelltenq () { local a,i,j,k,kt,ty1,ty2,ic1,spks1,spks2,c1,c2,gid1,gid2,nte,nshuf,fas\
                           localobj nqte,vtmp,ce1,ce2,vs0,vs1,vte,vnte,vfrom,vto
  if(numarg()>0)fas=$1 else fas=0
  nqte=new NQS("id1","id2","ty1","ty2","ic1","ic2","spks1","spks2","col1","col2","gid1","gid2","te","nte")
  {clrAllMyNQs() initAllMyNQs()}
  if(fas) {
    a=allocvecs(vtmp,vs0,vs1,vte,vnte,vfrom,vto) vrsz(4,vtmp)
    for c1=0,numcols-1 { ce1=col[c1].ce
      for i=0,nqCO[c1].m-1 vs0.x(i)=nqCO[c1].v[i].sum()
      for c2=0,numcols-1 { ce2=col[c2].ce
        print "C", c1, "-> C",c2
        vrsz(nqCO[c2].m,vs1) // save spike counts
        for i=0,nqCO[c2].m-1 vs1.x(i)=nqCO[c2].v[i].sum()
        for i=0,nqCO[c1].m-1 { // get from/to indices
          for j=0,nqCO[c2].m-1 { if(c1==c2 && j==i)continue
            {vfrom.append(i) vto.append(j)}
        vfrom.ntel2(vto,nqCO[c1].vl,nqCO[c2].vl,vnte,vte,nshuf) // calculate transfer entropy
        k=0 //index into vte,vnte
        for i=0,nqCO[c1].m-1 { // save pairwise nTE values and other info into NQS
          if(i%100==0) printf("i=%d\n",i)
          {ty1=ce1.o(i).type ic1=ice(ty1) spks1=vs0.x(i) gid1=ce1.o(i).gid}
          for j=0,nqCO[c2].m-1 {
            if(c1==c2 && j==i)continue
            {ty2=ce2.o(j).type spks2=vs1.x(j)}
  } else {
    a=allocvecs(vtmp,vs0,vs1) vrsz(4,vtmp)
    for c1=0,numcols-1 { ce1=col[c1].ce
      for i=0,nqCO[c1].m-1 vs0.x(i)=nqCO[c1].v[i].sum()
      for c2=0,numcols-1 { ce2=col[c2].ce
        print "C", c1, "-> C",c2
        vrsz(nqCO[c2].m,vs1) // save spike counts
        for i=0,nqCO[c2].m-1 vs1.x(i)=nqCO[c2].v[i].sum()
        for i=0,nqCO[c1].m-1 { // get pairwise nTE values
          if(i%100==0) printf("i=%d\n",i)
          {ty1=ce1.o(i).type ic1=ice(ty1) spks1=vs0.x(i) gid1=ce1.o(i).gid}
          for j=0,nqCO[c2].m-1 {
            if(c1==c2 && j==i)continue
            {ty2=ce2.o(j).type spks2=vs1.x(j)}
  for"te"),"nte") nqte.v[i].unnan()
  return nqte
// //* getincelltenqt() - get nqs of input->cell normalized transfer entropy vs time
// //getincelltenqt(bmin,bmax,binc,winsz)
// // bmin is min bin, bmax is max bin, binc is # of bins to move over by, winsz = # of bins to examine
// obfunc getincelltenqt () { local a,sy,id,intfid,idx,winsz,bmin,bmax,binc,bb,ty,pls,sh2\
//                         localobj nqt,vtmp,v1,v2,visi
// bmin=$1 bmax=$2 binc=$3 winsz=$4
// if(numarg()>4 && nqisi[CDX]!=nil)sh2=$5 else sh2=0
// sy=id=intfid=0
// nqin.tog("DB")
// nqt=new NQS("id","type","syn","t","nte","INTFid","pulse")
// initMyNQs()
// a=allocvecs(v1,v2,visi)
// vrsz(nqt.v[0].size,v1,v2)  vrsz(0,v1,v2)
// idx=0
// for nqin.qt(&id,"id",&sy,"sy",&intfid,"INTFid",&ty,"type",&pls,"pulse") {
//   if(sh2) {
//     nqisi[CDX].select("id",id)
//     visi.copy(nqisi[CDX].getcol("isi"))
//   }
//   for(bb=bmin;bb+winsz<=bmax;bb+=binc){
//     v1.copy(nqIN[CDX].v[intfid],bb,bb+winsz)
//     v2.copy(nqCO[CDX].v[id],bb,bb+winsz)
//     if(sh2) vtmp=normte(v1,v2,30,visi,binsz,htmin) else vtmp=normte(v1,v2,30)
//     vtmp.unnan()
//     nqt.append(id,ty,sy,htmin+bb*binsz,vtmp.x(2),intfid,pls)
//   }
//   if(idx%100==0) printf("idx=%d\n",idx)
//   idx+=1
// }
// dealloc(a)
// return nqt
// }

//* drawEIktau() $o1=nqc,[$2=pulse ID optional]
proc drawEIktau () { local x,pdx,nn localobj nqc,str
  if(numarg()>1)pdx=$2 else pdx=-1
  str=new String()
  if(pflg==0) g.label(0,0,"counts") else g.label(0,0,"probability")
  if(pdx>-1)"syn",AM,"pulse",pdx) else"syn",AM)
  sprint(str.s,"E input/output ktau, avg=%g",nqc.getcol("ktau").mean)
  if(pdx>-1)"syn",GA,"pulse",pdx) else"syn",GA)
  if(nn>0) {
    sprint(str.s,"I input/output ktau, avg=%g",nqc.getcol("ktau").mean)
//* drawEIcmp() -- draws histo comparing AM vs GA input->output relationship on column $s2
//$o1=nqs, $s2=column to perform comparison on,opt arg 3=pulse id, opt 4th arg=cell-type
proc drawEIcmp () { local x,pdx,nn localobj nq,str
  nq=$o1  str=new String()  g.color(1)
  if(pflg==0) g.label(0,0,"counts") else g.label(0,0,"probability")
  {ers=0   clr=3  g.color(clr)}
  if(numarg()>2)pdx=$3 else pdx=-1
  if(numarg()>3)"syn",AM,"type",$4) else"syn",AM)
  if(numarg()>3) {
    sprint(str.s,"%s:E input/output %s, avg=%g",CTYP.o($4).s,$s2,nq.getcol($s2).mean)
  } else {
    sprint(str.s,"E input/output %s, avg=%g",$s2,nq.getcol($s2).mean)
  if(numarg()>3)"syn",GA,"type",$4,"loc",SOMA) else"syn",GA,"loc",SOMA)
  {clr=2 g.color(clr)}
  if(numarg()>3) {
    sprint(str.s,"%s:SOMA I input/output %s, avg=%g",CTYP.o($4).s,$s2,nq.getcol($s2).mean)
  } else {
    sprint(str.s,"SOMA I input/output %s, avg=%g",$s2,nq.getcol($s2).mean)
  if(numarg()>3)"syn",GA,"type",$4,"loc",DEND) else"syn",GA,"loc",DEND)
  {clr=4 g.color(clr)}
  if(numarg()>3) {
    sprint(str.s,"%s:DEND I input/output %s, avg=%g",CTYP.o($4).s,$s2,nq.getcol($s2).mean)
  } else {
    sprint(str.s,"DEND I input/output %s, avg=%g",$s2,nq.getcol($s2).mean)

//* 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) {
    if(vf.x(i)>maxf) break
  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
  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)
      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)
  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
  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()
    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", sprint(str.t,"%s_LFP",}
    for ltr(vi,panobj.printlist,&idx) {
      if(!strcmp(,str.s)) {
        print "reading in spikes for ",, " from vitem ",
        panobj.rv_readvec(idx,vit[cdx].tvec,vit[cdx].vec) //read in raster
    if(FileExists(svq.s)) {
      print "reading in vq for ",, " from ", svq.s
    } else print "loadmydat WARNING: file ", svq.s , " doesn't exist."
    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))
  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
  if(numarg()>0) {
    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
  } else {
    vrsz(col[CDX].allcells,vf) vf.fill(0)
    for vtr(&i,vit[CDX].vec) vf.x(i)=1  
  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)
  {panobj.pvplist(st.s,strv) printf("saved vecs\n")}
  for cdx=0,numcols-1 { // save NQS with LFPs
    {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) {
    {mknclnq() nclnq.tog("DB")}
    printf("saved connectivity\n")
  if(svq) for ltr(co,lcol,&cdx) {
    {vq[cdx].tog("DB") vq[cdx].sv(st.s) printf("saved inputs for %s\n",}

//$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)

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[PubMed]

References and models cited by this paper

References and models that cite this paper

Adesnik H, Scanziani M (2010) Lateral competition for cortical space by layer-specific horizontal circuits. Nature 464:1155-60 [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]

Carnevale NT, Hines ML (2006) The NEURON Book

Cui J, Xu L, Bressler SL, Ding M, Liang H (2008) BSMART: a Matlab-C toolbox for analysis of multichannel neural time series. Neural Netw 21:1094-104 [PubMed]

Francis JT, Xu S, Chapin JK (2008) Proprioceptive and cutaneous representations in the rat ventral posterolateral thalamus. J Neurophysiol 99:2291-304 [PubMed]

Gisiger T, Boukadoum M (2011) Mechanisms Gating the Flow of Information in the Cortex: What They Might Look Like and What Their Uses may be. Front Comput Neurosci 5:1-304 [PubMed]

Hines ML, Carnevale NT (2001) NEURON: a tool for neuroscientists. Neuroscientist 7:123-35 [Journal] [PubMed]

   Spatial gridding and temporal accuracy in NEURON (Hines and Carnevale 2001) [Model]

Kamiński M, Ding M, Truccolo WA, Bressler SL (2001) Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance. Biol Cybern 85:145-57 [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]

Lizier JT, Heinzle J, Horstmann A, Haynes JD, Prokopenko M (2011) Multivariate information-theoretic measures reveal directed information structure and task relevant changes in fMRI connectivity. J Comput Neurosci 30:85-107

Lloyd KG, Davidson L, Hornykiewicz O (1975) The neurochemistry of Parkinson's disease: effect of L-dopa therapy. J Pharmacol Exp Ther 195:453-64 [PubMed]

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

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

Meyer JS, Obara K, Muramatsu K (1993) Diaschisis. Neurol Res 15:362-6 [PubMed]

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]

Quilodran R, Gariel MA, Markov NT, Falchier A, Vezoli J, Sallet J, Anderson JC, Dehay C, Doug (2008) Strong loops in the neocortex Society for Neuroscience Abstracts 853.4

Rasche D, Rinaldi PC, Young RF, Tronnier VM (2006) Deep brain stimulation for the treatment of various chronic pain syndromes. Neurosurg Focus 21:E8

Richman JS, Moorman JR (2000) Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol Heart Circ Physiol 278:H2039-49 [PubMed]

Schiff ND, Giacino JT, Kalmar K, Victor JD, Baker K, Gerber M, Fritz B, Eisenberg B, Biondi T (2007) Behavioural improvements with thalamic stimulation after severe traumatic brain injury. Nature 448:600-3 [PubMed]

Schroeder CE, Mehta AD, Foxe JJ (2001) Determinants and mechanisms of attentional modulation of neural processing. Front Biosci 6:D672-84

Shipp S (2005) The importance of being agranular: a comparative account of visual and motor cortex. Philos Trans R Soc Lond B Biol Sci 360:797-814 [PubMed]

Stefani A, Lozano AM, Peppe A, Stanzione P, Galati S, Tropepi D, Pierantozzi M, Brusa L, Scar (2007) Bilateral deep brain stimulation of the pedunculopontine and subthalamic nuclei in severe Parkinson's disease. Brain 130:1596-607 [PubMed]

Stoerig P, Cowey A (1997) Blindsight in man and monkey. Brain 120 ( Pt 3):535-59 [PubMed]

Traub RD, Contreras D, Cunningham MO, Murray H, Lebeau FE, Roopun A, Bibbig A, et al (2005) A single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles and epileptogenic bursts J Neurophysiol 93(4):2194-232 [Journal] [PubMed]

   A single column thalamocortical network model (Traub et al 2005) [Model]
   Collection of simulated data from a thalamocortical network model (Glabska, Chintaluri, Wojcik 2017) [Model]

Van Essen DC, Anderson CH, Felleman DJ (1992) Information processing in the primate visual system: an integrated systems perspective. Science 255:419-23 [PubMed]

Von_monakow C (1914) Die Lokalisation im Grosshirn und der Abbau der Funktion durch kortikale Herde

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 [Journal] [PubMed]

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

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]

Dura-Bernal S, Zhou X, Neymotin SA, Przekwas A, Francis JT, Lytton WW (2015) Cortical Spiking Network Interfaced with Virtual Musculoskeletal Arm and Robotic Arm. Front Neurorobot 9:13 [Journal] [PubMed]

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

Kerr CC, Van Albada SJ, Neymotin SA, Chadderdon GL, Robinson PA, Lytton WW (2013) Cortical information flow in Parkinson's disease: a composite network-field model. Front Comput Neurosci 7:39:1-14 [Journal] [PubMed]

   Composite spiking network/neural field model of Parkinsons (Kerr et al 2013) [Model]

Neymotin SA, Chadderdon GL, Kerr CC, Francis JT, Lytton WW (2013) Reinforcement learning of 2-joint virtual arm reaching in a computer model of sensorimotor cortex Neural Computation 25(12):3263-93 [Journal] [PubMed]

   Sensorimotor cortex reinforcement learning of 2-joint virtual arm reaching (Neymotin et al. 2013) [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, 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]

Rowan MS, Neymotin SA, Lytton WW (2014) Electrostimulation to reduce synaptic scaling driven progression of Alzheimer's disease. Front Comput Neurosci 8:39 [Journal] [PubMed]

   Electrostimulation to reduce synaptic scaling driven progression of Alzheimers (Rowan et al. 2014) [Model]

(38 refs)