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

 Download zip file   Auto-launch 
Help downloading and running models
Accession:154096
"... As cells die and synapses lose their drive, remaining cells suffer an initial decrease in activity. Neuronal homeostatic synaptic scaling then provides a feedback mechanism to restore activity. ... The scaling mechanism increases the firing rates of remaining cells in the network to compensate for decreases in network activity. However, this effect can itself become a pathology, ... Here, we present a mechanistic explanation of how directed brain stimulation might be expected to slow AD progression based on computational simulations in a 470-neuron biomimetic model of a neocortical column. ... "
Reference:
1 . Rowan MS, Neymotin SA, Lytton WW (2014) Electrostimulation to reduce synaptic scaling driven progression of Alzheimer's disease. Front Comput Neurosci 8:39 [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; 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; Python;
Model Concept(s): Long-term Synaptic Plasticity; Aging/Alzheimer`s; Deep brain stimulation; Homeostasis;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [samn at neurosim.downstate.edu]; Rowan, Mark [m.s.rowan at cs.bham.ac.uk];
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; Glutamate;
/
RowanEtAl2014
batchscripts
mod
README
alz.hoc
alzinfo.m
autotune.hoc *
basestdp.hoc *
batch.hoc *
batch2.hoc *
batchcommon
checkirreg.hoc *
clusterrun.sh
col.dot *
col.hoc *
comppowspec.hoc *
condisconcellfig.hoc *
condisconpowfig.hoc *
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
e2hubsdisconpow.hoc *
e2incconpow.hoc *
filtutils.hoc *
flexinput.hoc
geom.hoc *
graphplug.hoc *
grvec.hoc *
infot.hoc *
init.hoc *
labels.hoc *
load.hoc *
local.hoc *
makepopspikenq.hoc *
matfftpowplug.hoc *
matpmtmplug.hoc *
matpmtmsubpopplug.hoc *
matspecplug.hoc *
mosinit.hoc
network.hoc *
nload.hoc *
nqpplug.hoc *
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc *
params.hoc
plot.py
plotavg.py
plotbatch.sh
plotbatchcluster.sh
plotdeletions.py
plotntes.py
powchgtest.hoc *
pyhoc.py
python.hoc *
pywrap.hoc *
ratlfp.dat *
redE2.hoc *
run.hoc
runsim.sh
setup.hoc *
shufmua.hoc *
sim.hoc
simctrl.hoc *
spkts.hoc *
stats.hoc *
syncode.hoc *
vsampenplug.hoc *
writedata.hoc
xgetargs.hoc *
                            
// $Id: basestdp.hoc,v 1.229 2012/02/01 12:30:33 samn Exp $

//* params/variables

declare("TargRate",8) // target peak location
declare("syl","o[2]") // list of sywvs
declare("taul",new List(),"incl",new List(),"maxwl",new List(),"wgl",new List())
declare("nqrat","o[1]","lnqp",new List(),"mynqp","o[1]")
declare("ltab","o[1]","contab","o[1]","wtab","o[2]","mtab","o[1]","dvtab","o[1]")
declare("witer",0,"updinc",2e3,"uprob",0.1)
declare("EEinc",0.01,"EIinc",0.01,"IEinc",0.01,"IIinc",0.01,"skipI",0)
declare("vrecw",new Vector(),"nqrec","o[1]")
declare("MODW",1,"MINW",1,"MODINC",1,"MININC",0.01)
declare("SPECTY",3,"drawsm",0,"smhz",1,"PRESM",1) // whether to draw smoothed spectra in mkdrawspec
declare("dosend",1,"fith","o[3]","dosendex",0,"nqwrex","o[1]")
declare("maxERate",2,"maxIRate",200)
declare("lasterr",0)

declare("restorewg",0)

CTYP.append(new String("E"))
CTYP.append(new String("I"))

//* updateplast(col)
proc updateplast () { local i,j,a localobj col,xo,vwg,vtau,vinc,vmaxw,vidx
  a=allocvecs(vwg,vtau,vinc,vmaxw,vidx) col=$o1
  for ltr(xo,col.ce) {
    xo.getdvi(vidx)
    vrsz(vidx.size,vwg,vtau,vinc,vmaxw)
    xo.getplast(vwg,vtau,vinc,vmaxw)
    for vtr(&i,vidx,&j) if(vmaxw.x(j)) {
      if(ice(col.ce.o(i))) {
        vmaxw.x(j) = plastEImaxw
        vinc.x(j) = plastEIinc
      } else {
        vmaxw.x(j) = plastEEmaxw
        vinc.x(j) = plastEEinc
      }
    }
    xo.setplast(vwg,vtau,vinc,vmaxw)
  }
  dealloc(a)
}

if(restorewg) {
  nqplast=new NQS("/u/samn/intfnips/data/11may30_basestdp145_TrainW10_TargTrainRate8_LD2e3__plastnq.nqs")
  setplastnq(nqplast,col)
  resetplast_INTF6 = 0 // make sure starting params kept
}

// jrtm_INTF6 = updinc

declare("vcheck",new Vector())
declare("myncl",new List(),"myspkl",new List(),"myspktyl",new List(),"vice",new Vector(col.allcells))

//* declares

declare("wgnq","o[9]","myv","o[10]")
declare("nqp","o[9][10]","nqps","o[9][10]","nqchg","o[9]","nqchgs","o[9]")
declare("Ts","d[10]","Te","d[10]")

//* mkdrawspec(startt,endt,celltype[,norm,sm])  - makes pmtm spectra
obfunc mkonespec () { local i,ct,k,nrm,sampr,sm,ts,te localobj ls,nqp,nqps
  ls=new List()
  {ts=$1*1e3 te=$2*1e3 ct=$3}
  if(numarg()>3) nrm=$4 else nrm=0
  if(numarg()>4) sm=$5 else sm=201
  sampr=1e3/binsz
  if(ct<0) sampr=1e3/vdt_INTF6
  for CDX=0,numcols-1 {
    if(myv==nil) myv=new Vector() else myv.resize(0)
    if(ct==-1) {
      {myv.copy(nqLFP[CDX].v,ts/vdt_INTF6,te/vdt_INTF6-1) myv.sub(myv.mean)}
    } else {
      {myv.copy(nqCTY[CDX].v[ct],ts/binsz,te/binsz-1) myv.sub(myv.mean)}
    }
    nqp=getspecnq(myv,sampr,SPECTY,PRESM)
    {nqps=new NQS() nqps.cp(nqp) boxfilt(nqps.v[1],sm)}
    {ls.append(nqp) ls.append(nqps)}
  }
  return ls
}

//* mk4specs(start1,end1,start2,end2,start3,end3,start4,end4,celltype[,norm,sm])  - makes pmtm spectra for 
// the 4 periods of the sim: baseline, zip on, zip on + learning, zip on + recall
proc mk4specs () { local i,ct,k,nrm,sampr,sm,smHZ,a localobj vec
  a=allocvecs(vec)
  {Ts[0]=$1*1e3 Te[0]=$2*1e3}
  {Ts[1]=$3*1e3 Te[1]=$4*1e3}
  {Ts[2]=$5*1e3 Te[2]=$6*1e3}
  {Ts[3]=$7*1e3 Te[3]=$8*1e3}
  ct=$9
  if(numarg()>9) nrm=$10 else nrm=0
  if(numarg()>10) smHZ=$11 else smHZ=smhz
  sampr=1e3/binsz
  if(ct<0) sampr=200
  for i=0,3 if(Te[i] > Ts[i]) {
    if(myv[i]==nil) myv[i]=new Vector() else myv[i].resize(0)
    if(ct==-1) {
      {myv[i].copy(nqLFP.v,Ts[i]/vdt_INTF6,Te[i]/vdt_INTF6-1)}
      vec.resize(myv[i].size/(1e3/vdt_INTF6/sampr))
      myv[i].downsampavg(vec,(1e3/vdt_INTF6/sampr)) // downsample to 200 Hz
      myv[i].resize(0)
      myv[i].copy(vec)
      myv[i].sub(myv[i].mean)
    } else {
      {myv[i].copy(nqCTY[CDX].v[ct],Ts[i]/binsz,Te[i]/binsz-1) myv[i].sub(myv[i].mean)}
    }
    {nqsdel(nqp[CDX][i]) nqp[CDX][i]=getspecnq(myv[i],sampr,SPECTY,PRESM)}
    sm=0
    while(nqp[CDX][i].v.x(sm) < smHZ) sm += 1
    //print nqp[CDX][i].v.x(sm),smHZ,sm
    {nqsdel(nqps[CDX][i]) nqps[CDX][i]=new NQS()}
    {nqps[CDX][i].cp(nqp[CDX][i]) boxfilt(nqps[CDX][i].v[1],sm)}
    if(g!=nil) {
      if(nrm) nqps[CDX][i].v[1].div(nqps[CDX][i].v[1].sum)
      if(ct>=0) {
        print CTYP.o(ct).s,i,nqps[CDX][i].v.x(nqps[CDX][i].v[1].max_ind),nqps[CDX][i].v[1].max
      } else print "LFP",i,nqps[CDX][i].v.x(nqps[CDX][i].v[1].max_ind),nqps[CDX][i].v[1].max
    }
  }
  dealloc(a)
}

//* dr4spec(celltype) 
proc dr4spec () {
  initg()
  mk4specs(1,PreDur,\
           PreDur+1,PreDur+ZipDur-1,\
           PreDur+ZipDur+1,PreDur+ZipDur+LearnDur-1,\
           PreDur+ZipDur+LearnDur+1,PreDur+ZipDur+LearnDur+PostDur-1,$1)
  for i=0,3 if(g!=nil) {
    if(drawsm) {
      if(nqps[CDX][i]==nil) continue
      nqps[CDX][i].gr("pow","f",0,i+1,1)
    } else {
      if(nqp[CDX][i]==nil) continue
      nqp[CDX][i].gr("pow","f",0,i+1,1)
    }
  }
  fing()
}

//* dr4chg(celltype)
proc dr4chg () { local i,p1,p2
  initg()
  if(numarg()<=1) {
    mk4specs(1,PreDur,\
             PreDur+1,PreDur+ZipDur-1,\
             PreDur+ZipDur+1,PreDur+ZipDur+LearnDur-1,\
             PreDur+ZipDur+LearnDur+1,PreDur+ZipDur+LearnDur+PostDur-1,$1)
  }
  {nqsdel(nqchg[CDX]) nqchg[CDX]=new NQS("f","pow")  }
  nqp[CDX][0].verbose=nqp[CDX][2].verbose=0
  for i=0,int(nqp[CDX][0].v.max)-1 {
    if(!nqp[CDX][0].select("f","[]",i,i+1)) continue
    if(!nqp[CDX][2].select("f","[]",i,i+1)) continue
    p2 = nqp[CDX][2].getcol("pow").sum
    p1 = nqp[CDX][0].getcol("pow").sum
    if(p1>0) nqchg[CDX].append((i+1)/2.0,(p2-p1)/p1) else  nqchg[CDX].append((i+1)/2.0,0)
  }
  if(g!=nil) nqchg[CDX].gr("pow","f",0,1,1)
  fing()
  nqp[CDX][0].verbose=nqp[CDX][2].verbose=1
}
//* getratety(start,end,ty)
func getratety () { local startt,endt,ty,ns,r
  startt=$1 endt=$2 ty=$3
  if(endt-startt<=0) return 0
  if(snq[CDX]==nil) snq[CDX]=SpikeNQS(vit[CDX].tvec,vit[CDX].vec,col[CDX])
  ns=snq[CDX].select("type",ty,"t","[]",startt,endt)
  r = 1e3 * ns / (col[CDX].numc[ty] * (endt-startt) )
  return r
}
//* pravgrates2 (start1,end1,start2,end2,start3,end3,start4,end4) 
proc pravgrates2 () { local i,ct,r localobj s
  {Ts[0]=$1*1e3 Te[0]=$2*1e3}
  {Ts[1]=$3*1e3 Te[1]=$4*1e3}
  {Ts[2]=$5*1e3 Te[2]=$6*1e3}
  {Ts[3]=$7*1e3 Te[3]=$8*1e3}
  s=new String2()
  for col.ctt(&ct) {
    sprint(s.s,"%s:\t",CTYP.o(ct).s)
    for i=0,3 {
      r = getratety(Ts[i],Te[i],ct)
      sprint(s.t,"%0.2f\t",r)
      strcat(s.s,s.t)
    }
    print s.s
  }
}

//* savenqspec - saves nqp,nqps spectra, uses strv
proc savenqspec () { local i,st,et,ct localobj str
  str=new String() 
  st=plaststartT_INTF6/1e3
  et=plastendT_INTF6/1e3
  for case(&ct,E2,I2,E4,I4,CTYPi,CTYPi+1) {
    if(g!=nil) {
      g.erase
      gvmarkflag=0 
    }
    mkdrawspec(0,TrainStart/1e3,TrainStart/1e3,LearnDur+TrainStart/1e3,LearnDur+TrainStart/1e3,tstop/1e3,ct)
    for i=0,2 {    
      {sprint(str.s,"./data/%s_%s_nqp_%d.nqs",strv,CTYP.o(ct).s,i) nqp[0][i].sv(str.s)}
      {sprint(str.s,"./data/%s_%s_nqps_%d.nqs",strv,CTYP.o(ct).s,i) nqps[0][i].sv(str.s)}
    }
  }
}
//* ldnqspec - saves nqp,nqps spectra, uses strv
obfunc ldnqspec () { local i,st,et,ct localobj str,nqo,nq,nqs
  {str=new String() nqo=new NQS("str","nq")  nqo.odec("nq") nqo.strdec("str")}
  for case(&ct,E2,I2,E4,I4,CTYPi,CTYPi+1) {
    for i=0,3 {    
      {sprint(str.s,"./data/%s_%s_nqp_%d.nqs",strv,CTYP.o(ct).s,i) nq=new NQS(str.s)}
      nqo.append(str.s,nq)
      nqsdel(nq)
      {sprint(str.s,"./data/%s_%s_nqps_%d.nqs",strv,CTYP.o(ct).s,i) nqs=new NQS(str.s)}
      nqo.append(str.s,nqs)
      nqsdel(nqs)      
    }
  }
  return nqo
}
//* Vec2Txt - save Vector to text file
func Vec2Txt () { local i localobj fp
  fp=new File()
  fp.wopen($s2)
  if(!fp.isopen) return 0
  for i=0,$o1.size-1 fp.printf("%g\n",$o1.x(i))
  fp.close()
  return 1
}

//* run and init nqs objects
proc myrun () { local i localobj xo
  run()
  {nqsdel(snq[CDX]) snq[CDX]=SpikeNQS(vit[CDX].tvec,vit[CDX].vec,col[CDX])}
  for CDX=0,numcols-1 {snq[CDX].marksym="O"}
  CDX=0
  {nqsdel(nqplast) nqplast=getplastnq(col)}
  {nqsdel(mynqp) mynqp=lnq2nqs(lnqp)}
  pravgrates2(1,PreDur,\
              PreDur+1,PreDur+ZipDur-1,\
              PreDur+ZipDur+1,PreDur+ZipDur+LearnDur-1,\
              PreDur+ZipDur+LearnDur+1,PreDur+ZipDur+LearnDur+PostDur-1)
  initAllMyNQs()
}

//* mysv - save output after myrun
proc mysv () { local cdx localobj s
  s=new String()
  {sprint(s.s,"/u/samn/intfzip/data/%s_snq.nqs",$s1) snq.tog("DB") snq.sv(s.s)}
  if(nqrat!=nil) {
    {sprint(s.s,"/u/samn/intfzip/data/%s_nqrat.nqs",$s1) nqrat.tog("DB") nqrat.sv(s.s)}
  }
  if(mynqp.size>0) {
    {sprint(s.s,"/u/samn/intfzip/data/%s_mynqp.nqs",$s1) mynqp.tog("DB") mynqp.sv(s.s)}
  }
  // save NQS with LFPs
  {cdx=0 sprint(s.s,"/u/samn/intfzip/data/%s_%s_LFP.nqs",$s1,col[cdx].name)}
  {nqLFP[cdx].tog("DB") nqLFP[cdx].sv(s.s)}
  {nqsdel(wgnq) wgnq=mkwgnq(col) sprint(s.s,"/u/samn/intfzip/data/%s_wgnq.nqs",strv) wgnq.sv(s.s)}
}
//* myrunsv(simstr) - run & save output
proc myrunsv () { 
  myrun()
  mysv($s1)
}
//* snq2vit(snq,vit) -  copy snq into a vitem
proc snq2vit () { local i localobj snq,vit  
  snq=$o1 vit=$o2 snq.tog("DB")
  vrsz(0,vit.vec,vit.tvec)
  vit.tvec.copy(snq.getcol("t"))
  vit.vec.copy(snq.getcol("id"))
}
//* myrd(simstr[,rdr]) - read output from data saved with mysv
proc myrd () { local rdr,cdx localobj s
  s=new String() if(numarg()>1) rdr=$2 else rdr=0
  {nqsdel(snq) sprint(s.s,"/u/samn/intfzip/data/%s_snq.nqs",$s1) snq=new NQS(s.s) snq2vit(snq,vit)}
  if(rdr){
    {nqsdel(nqrat) sprint(s.s,"/u/samn/intfzip/data/%s_nqrat.nqs",$s1) nqrat=new NQS(s.s)}
    {nqsdel(mynqp) sprint(s.s,"/u/samn/intfzip/data/%s_mynqp.nqs",$s1) mynqp=new NQS(s.s)}
  }
  cdx=0
  sprint(s.s,"/u/samn/intfzip/data/%s_%s_LFP.nqs",$s1,col[cdx].name)
  {nqsdel(nqLFP[cdx]) nqLFP[cdx]=new NQS(s.s)}
  {nqsdel(wgnq) sprint(s.s,"/u/samn/intfzip/data/%s_wgnq.nqs",$s1) wgnq=new NQS(s.s)}
}

//* settunerc - setup recording of spikes used in tuning
proc settunerec () { local i localobj xo,nc
  for i=0,CTYPi-1 myspktyl.append(new Vector())
  for ltr(xo,col.ce,&i) {
    xo.flag("out",1) // make sure NetCon events enabled from this cell
    myncl.append(nc=new NetCon(xo,nil))
    myspkl.append(new Vector())
    nc.record(myspkl.o(i)) // record each cell separately
    vice.x(i)=ice(xo.type)
  }
}

//* mksyl - setup lists of weight vectors
proc mksyl () { local i,dvt localobj vw1,vw2
  for i=0,1 syl[i]=new List()
  for i=0,col.allcells-1 {
    dvt=col.ce.o(i).getdvi()
    vw1=new Vector(dvt)
    vw2=new Vector(dvt)
    col.ce.o(i).getsywv(vw1,vw2)
    syl[0].append(vw1)
    syl[1].append(vw2)
  }
}
//* mkstdpl - 
proc mkstdpl () { local i,dv localobj vtau,vinc,vmaxw,vwg,xo
  for ltr(xo,col.ce,&i) {
    dv=xo.getdvi
    vtau=new Vector(dv)
    vinc=new Vector(dv)
    vmaxw=new Vector(dv)
    vwg=new Vector(dv)
    if(!ice(xo.type)) { // only STDP from E->X cells
      xo.getplast(vwg,vtau,vinc,vmaxw)
      vcheck.append(i)
    }
    taul.append(vtau)
    incl.append(vinc)
    maxwl.append(vmaxw)
    wgl.append(vwg)
  }
}
//* mkdvtab - make table with dvi
proc mkdvtab () { local i localobj col,xo
  col=$o1
  dvtab=new List()
  for ltr(xo,col.ce,&i) {
    dvtab.append(new Vector(xo.getdvi))
    xo.getdvi(dvtab.o(dvtab.count-1))
  }  
}

//* conn2tab - make lookup tables with connectivity info
obfunc conn2tab () { local i,j,k,id1,id2 localobj ltab,col,nqc,contab,wtab1,wtab2,mtab,vc
  col=$o1 ltab=new List() vc=new Vector(col.allcells)
  for i=0,3 ltab.append(new Matrix(col.allcells,col.allcells))
  {contab=ltab.o(0) wtab1=ltab.o(1) wtab2=ltab.o(2) mtab=ltab.o(3)}
  if(col.connsnq==nil) {
    print "conn2tab ERR: col.connsnq is nil"
    return nil
  }
  nqc=col.connsnq
  nqc.sort("del") // make sure order in NQS corresponds to getdvi order
  for i=0,nqc.v.size-1 {
    id1=nqc.v[0].x(i) // from id1
    id2=nqc.v[1].x(i) // to id2
    contab.x(id1,id2)=1 // is there a connection?
    wtab1.x(id1,id2)=nqc.v[4].x(i) // weight 1
    wtab2.x(id1,id2)=nqc.v[5].x(i) // weight 2
    mtab.x(id1,id2) = vc.x(id1) // index into div vector -- assumes order in connsnq according to div
    vc.x(id1) += 1
  } 
  return ltab
}

//* wrex2nq(col) - return an NQS with external input info
obfunc wrex2nq () { local i,sy,r,w,row,a localobj enq,col,nsl,ncl,nc,vsy,vm
  a=allocvecs(vsy,vm) vsy.append(GA,GA2,AM2,NM2) vm.append(0,0,1,0,3,4,2)
  col=$o1 ncl=col.cstim.ncl nsl=col.cstim.nsl
  enq = new NQS("id","wGA","wGA2","wAM2","wNM2","rGA","rGA2","rAM2","rNM2","bal")
  enq.v.indgen(0,col.allcells-1,1)
  enq.pad()
  for i=1,enq.m-1 enq.v[i].fill(0)
  for ltr(nc,ncl,&i) {
    row = nc.syn.id
    for vtr(&sy,vsy) if(nc.weight(sy)) {
      w = enq.v[vm.x(sy)].x(row) = nc.weight(sy)
      r = enq.v[vm.x(sy)+4].x(row) = 1e3/nsl.o(i).interval
      if(sy==AM2 || sy==NM2) {
        enq.v[9].x(row) += r*w
      } else {
        enq.v[9].x(row) -= r*w
      }
      break
    }
  }
  dealloc(a)
  return enq
}

//* nq2wrex(enq,col) - sets external weight/rate params based on NQS
proc nq2wrex () { local i,sy,row,a localobj enq,col,nsl,ncl,nc,vsy,vm
  a=allocvecs(vsy,vm) vsy.append(GA,GA2,AM2,NM2) vm.append(0,0,1,0,3,4,2)
  enq = $o1 col=$o2 ncl=col.cstim.ncl nsl=col.cstim.nsl
  for ltr(nc,ncl,&i) {
    row = nc.syn.id
    for vtr(&sy,vsy) if(nc.weight(sy)) {
      nc.weight(sy) = enq.v[vm.x(sy)].x(row) 
      nsl.o(i).interval = 1e3 / enq.v[vm.x(sy)+4].x(row) 
      break
    }
  }
  dealloc(a)
}

//* reinforcement learning update
proc RLUpdate () { local i,j,k,md,df,fctr,inc,inc0,idx,trg,ety,cidx,pkx,pky,pkd,a,poid,pkdr,err,epiE,epiI\
                    localobj xo,vs,ce,nqp,vec,tvec,vh,vwg,vtau,vinc,vmaxw,vrate,vE,vI
  print "t:", t, ", witer:",witer
  a=allocvecs(vec,tvec,vh,vwg,vtau,vinc,vmaxw,vrate,vE,vI)
  ce=col.ce
  for i=0,CTYPi-1 if(col.numc[i]) myspktyl.o(i).resize(0) //setup per-type counts
  vrate.resize(ce.count)
  err = epiE = epiI = avgE = avgI = 0
  for ltr(xo,ce,&i) {
    myspktyl.o(xo.type).append(myspkl.o(i))
    vrate.x(i)=1e3*myspkl.o(i).size/updinc
    if(ice(xo.type)) vI.append(vrate.x(i)) else vE.append(vrate.x(i))
  }
  if(vE.mean() > maxERate) epiE=1
  if(vI.mean() > maxIRate) epiI=1
  vrsz(0,vec,tvec)
  for case(&i,E2,E4,E5R,E5B,E6,&j) tvec.append(myspktyl.o(i))
  vh = tvec.histogram(witer*updinc,(witer+1)*updinc,binsz)
  vh.sub(vh.mean)  
  lnqp.append(nqp=getspecnq(vh,1e3/binsz,SPECTY,PRESM)) // get spectrum from E MUA
  pkx = nqp.v.x(nqp.v[1].max_ind)
  pky = nqp.v[1].max
  pkd = TargRate - pkx // difference in peak frequency, + = too slow, - = too fast
  pkdr = abs(pkd) / TargRate

  for i=0,nqp.v.size-1 {
    if(abs(nqp.v.x(i)-TargRate) < 3) {
      err -= nqp.v[1].x(i)
    } else err += nqp.v[1].x(i)
  }

  // err += pkd^2

  if(err <= lasterr && !epiE && !epiI) {
    col.ce.o(0).dopelearn(1)  // LTP
    if(batch_flag==0) print pkx, pky, lasterr, err, " LTP"
  } else { 
    col.ce.o(0).dopelearn(-1) // LTD
    if(batch_flag==0) print pkx, pky, lasterr, err, " LTD"
  } 

  lasterr = err

  witer += 1
  for i=0,myspkl.count-1 myspkl.o(i).resize(0) // reset spike counts for cells to 0
  for(i=CTYPi-1;i>=0;i-=1) if(col.numc[i]) {
    j=1e3*myspktyl.o(i).size/(col.numc[i]*updinc)
    if(batch_flag==0) print CTYP.o(i).s, " " , j , " avg Hz "
    nqrat.append(t,witer,i,j,pkx,pky,err)
    myspktyl.o(i).resize(0) // reset spike counts for types to 0
  }
  dealloc(a)
  cvode.event(t+updinc,"RLUpdate()") // set next update weights event
}

//* updateSTDP - update the STDP params
proc updateSTDP () { local i,j,k,md,df,fctr,inc,inc0,idx,trg,ety,cidx,pkx,pky,pkd,a,poid,pkdr,err\
                    localobj xo,vs,ce,nqp,vec,tvec,vh,vwg,vtau,vinc,vmaxw,vrate
  print "t:", t, ", witer:",witer
  a=allocvecs(vec,tvec,vh,vwg,vtau,vinc,vmaxw,vrate)
  ce=col.ce
  for i=0,CTYPi-1 if(col.numc[i]) myspktyl.o(i).resize(0) //setup per-type counts
  vrate.resize(ce.count)
  for ltr(xo,ce,&i) {
    myspktyl.o(xo.type).append(myspkl.o(i))
    vrate.x(i)=1e3*myspkl.o(i).size/updinc
  }
  vrsz(0,vec,tvec)
  for case(&i,E2,E4,E5R,E5B,E6,&j) tvec.append(myspktyl.o(i))
  vh = tvec.histogram(witer*updinc,(witer+1)*updinc,binsz)
  vh.sub(vh.mean)  
  lnqp.append(nqp=getspecnq(vh,1e3/binsz,SPECTY,PRESM)) // get spectrum from E MUA
  pkx = nqp.v.x(nqp.v[1].max_ind)
  pky = nqp.v[1].max
  pkd = TargRate - pkx // difference in peak frequency, + = too slow, - = too fast
  pkdr = abs(pkd) / TargRate
  err = pkd^2
  if(batch_flag==0) print pkx, pky, pkd, pkdr
  if(MODINC || MODW) {
    if(t >= plaststartT_INTF6 && t <= plastendT_INTF6) for vtr(&i,vcheck) if(myspkl.o(i).size) { xo = ce.o(i)    
      vrsz(xo.getdvi,vwg,vtau,vinc,vmaxw)
      xo.getplast(vwg,vtau,vinc,vmaxw) // get current weight gains    
      for vtr(&j,dvtab.o(i),&k) {
        inc = 0
        if(vice.x(i)) { // presynaptic I cell
          if(ISTDP_INTF6) {
            if(vice.x(j)) { // postsynaptic I cell
              if(skipI) continue
              if(pkd > 0) inc = IIinc else if(pkd < 0) inc = -IIinc
            } else { // postsynaptic E cell
              if(pkd > 0) inc = IEinc else if(pkd < 0) inc = -IEinc
            }
          }
        } else { // presynaptic E cell
          if(ESTDP_INTF6) {
            if(vice.x(j)) { // postsynaptic I cell
              if(skipI) continue
              if(pkd > 0 && vrate.x(j)<maxIRate) inc = EIinc else if(pkd < 0) inc = -EIinc
            } else { // postsynaptic E cell
              if(pkd > 0) inc = -EEinc else if(pkd < 0 && vrate.x(j)<maxERate) inc = EEinc
            }
          }
        }
        if(MODINC) vinc.x(k) = MAXxy(MININC,vinc.x(k) + inc*pkdr)
        if(MODW) vmaxw.x(k) = MAXxy(MINW,vmaxw.x(k) + inc*pkdr)
      }
      xo.setplast(vwg,vtau,vinc,vmaxw) // reset plasticity params    
    }
  }
  for vtr(&i,vrecw) {
    for j=0,col.allcells-1 if(contab.x(i,j)) {
//      idx = mtab.x(i,j)
//      nqrec.append(i,j,ce.o(i).type,ce.o(j).type,syl[0].o(i).x(idx),syl[1].o(i).x(idx),witer)
    }
  }
  witer += 1
  for i=0,myspkl.count-1 myspkl.o(i).resize(0) // reset spike counts for cells to 0
  for(i=CTYPi-1;i>=0;i-=1) if(col.numc[i]) {
    j=1e3*myspktyl.o(i).size/(col.numc[i]*updinc)
    if(batch_flag==0) print CTYP.o(i).s, " " , j , " avg Hz "
    nqrat.append(t,witer,i,j,pkx,pky,err)
    myspktyl.o(i).resize(0) // reset spike counts for types to 0
  }
  dealloc(a)
  cvode.event(t+updinc,"updateSTDP()") // set next update weights event
}

//* mysend - starts off the update q
proc mysend () { local sz localobj xo
  for ltr(xo,lnqp) nqsdel(xo)
  lnqp.remove_all()
  sz = (LearnDur/(updinc/1e3))
  {nqsdel(nqrat) nqrat=new NQS("t","witer","ty","rate","pkx","pky","err") nqrat.clear(sz)}
  {nqsdel(nqrec) nqrec=new NQS("id1","id2","ty1","ty2","witer","wgain","inc","maxw")}
  lasterr=0
  if(DopeRL) {
    cvode.event(updinc,"RLUpdate()")
  } else cvode.event(updinc,"updateSTDP()") 
}

//* updateEX
proc updateEX () { local i
  print "updateEX at t = ", t
  for i=1,4 nqwrex.v[i].mul(0.99)
  nq2wrex(nqwrex,col)
  cvode.event(t+updinc,"updateEX()") // set next update weights event
}

//* mysendex - starts off the update EX q
proc mysendex () { local sz localobj xo
  nqsdel(nqwrex)
  nqwrex = wrex2nq(col)
  cvode.event((PreDur+ZipDur)*1e3,"updateEX()") 
}

//* prnqrat(nqrat) - print average peak rates in 20 iters
proc prnqrat () { local sdx localobj nqrat
  nqrat=$o1
  sdx=0
  nqrat.verbose=0
  while(sdx<witer) {
    nqrat.select("witer","[]",sdx,sdx+20,"ty",E5R)
    print nqrat.getcol("pkx").mean,"+/-",nqrat.getcol("pkx").stderr
    sdx += 20
  }
  nqrat.verbose=1
}
//* drtargrate(nqrat) -
proc drtargrate () { localobj nqrat
  nqrat=$o1
  nqrat.tog("DB")
  nqrat.select("ty",E5R)
  {gvmarkflag=0 nqrat.gr("pkx","t",0,2,1)}
  {gvmarkflag=1 nqrat.gr("pkx","t",0,2,3) gvmarkflag=0}
  drline(0,TargRate,tstop,TargRate,g,3,4)
  drline(TrainStart,0,TrainStart,100,g,9,9)
  drline(TrainStop,0,TrainStop,100,g,9,9)
}
//* drpoprates(nqrat) - draw population rates vs iteration
proc drpoprates () { local ct,i localobj nqrat
  nqrat=$o1
  nqrat.verbose=0
  drline(0,TargRate,tstop,TargRate,g,3,4)
  for case(&ct,E2,E4,E5R,E5B,E6,I2,I2L,I4,I4L,&i) {
    nqrat.select("ty",ct)
    nqrat.gr("rate","t",0,i+1,1)
  }
  nqrat.verbose=1
  drline(TrainStart,0,TrainStart,100,g,9,9)
  drline(TrainStop,0,TrainStop,100,g,9,9)
}
//* drlnqp(startidx,endidx)
proc drlnqp () { local sidx,eidx,i
  sidx=$1 eidx=$2
  for i=sidx,eidx {
    lnqp.o(i).gr("pow","f",0,(i+1)%10,1)
  }
}
//* mynqp2lnqp(mynqp) - copy the nqp NQS objects in $o1 to a list and return it
obfunc mynqp2lnqp () { local i localobj lnqp,mynqp,nqp
  mynqp=$o1
  mynqp.tog("DB")
  lnqp=new List()
  for i=0,mynqp.v.size-1 lnqp.append(nqp=mynqp.get("nqp",i).o)
  return lnqp
}
//* lnq2nqs(ls) - copy the ls NQS objects to a NQS and return it
obfunc lnq2nqs () { local i localobj lnq,nq,xo
  lnq=$o1
  {nq=new NQS("i","nqp") nq.odec("nqp") nq.clear(lnq.count)}
  for ltr(xo,lnq,&i) nq.append(i,xo)
  return nq
}
//* mkavgnqp(lnqp,startidx,endidx) - return an nqs with average +/ stderr of power at frequencies
// lnqp is a list of nqp objects
obfunc mkavgnqp () { local i,sidx,eidx localobj nqa,nqp,lnqp,va,ve
  lnqp=$o1 sidx=$2 eidx=$3
  nqa=new NQS("f","avg","err")
  for i=sidx,eidx {
    nqp=lnqp.o(i)
    
  }
  return nqa
}
//* svimg1(basename) - saves idraw, for when screen capture not working
proc svimg1 () {
  {sprint(tstr,"gif/%s.id",$s1) save_idraw(tstr)}
}
//* svimg2(basename) - saves id as gif & gets rid of intermediate pdf
proc svimg2 () {
  {sprint(tstr,"epstopdf gif/%s.id",$s1) system(tstr)}
  {sprint(tstr,"convert  gif/%s.pdf gif/%s.gif",$s1,$s1) system(tstr)}
  {sprint("rm gif/%s.pdf",$s1) system(tstr)}
  {sprint("rm gif/%s.id",$s1) system(tstr)}
}

//* calls

if(dosend) fith=new FInitializeHandler("mysend()")
if(dosendex) fith[1]=new FInitializeHandler("mysendex()")

settunerec()
mkstdpl()
mkdvtab(col)

ltab=conn2tab(col)
contab=ltab.o(0)
wtab[0]=ltab.o(1)
wtab[1]=ltab.o(2)
mtab=ltab.o(3)