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

 Download zip file   Auto-launch 
Help downloading and running models
Accession:141505
This model is an extension of a model (<a href="http://senselab.med.yale.edu/ModelDB/ShowModel.asp?model=138379">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.
Reference:
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;
Gene(s):
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 neurosim.downstate.edu]; Neymotin, Sam [samn at neurosim.downstate.edu]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
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;
/
neuroprosthesis
README
infot.mod *
intf6_.mod *
intfsw.mod *
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
batch.hoc
boxes.hoc
bsmart.py
col.hoc
comparecausality.py
comparerasters.py
declist.hoc
decmat.hoc *
decnqs.hoc *
decvec.hoc
default.hoc *
drline.hoc *
filtutils.hoc
flexinput.hoc
grvec.hoc
infot.hoc *
init.hoc
intfsw.hoc
labels.hoc
local.hoc *
misc.h *
mosinit.hoc
network.hoc
nload.hoc
nqs.hoc
nqsnet.hoc
nrnoc.hoc
params.hoc
pyhoc.py
ratlfp.dat *
run.hoc
runsim
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
xgetargs.hoc *
                            
// $Id: syncode.hoc,v 1.412 2011/07/29 18:56:31 billl Exp $

print "Loading syncode.hoc..."

//* setup
{load_file("grvec.hoc")}
{load_file("labels.hoc")}
{load_file("nqs.hoc")}
strdef syn1,syn2
thresh = -20
objref cp,svs,ncq,intf,vq
objref ivspks,vspks,wvspks,ncl[1][1][1]
objref sp[3], c[1], ce, nc[1], cells[10] // enough room for 10 cell types
objref vite
double numc[CTYPi],ix[CTYPi],ixe[CTYPi],wmat[CTYPi][CTYPi][STYPi],wd0[CTYPi][CTYPi][STYPi]
Incol=2
ncl = new List()
sp = new NQS()

//* iterator ixt(),ctt(),stt()
// ctt() iterates over all cells; eg for ctt() print ii
iterator ctt () { local jj // global ii
  if (argtype(2)==3) $&2=0 else i1=0
  for jj=0,CTYPi-1 if (numc[jj]!=0) { 
    if (ce!=nil) if (ix[jj]<ce.count()) { 
      if (argtype(1)==1) $o1=ce.o(ix[jj]) else XO=ce.o(ix[jj]) 
    }
    if (argtype(1)==3) $&1=jj else ii=jj
    iterator_statement
    if (argtype(2)==3) $&2+=1 else i1+=1
  }
}

//* iterator cttr() -- reverse iterator of ctt
iterator cttr () { local jj // global ii
  if (argtype(2)==3) $&2=0 else i1=0
  for(jj=CTYPi-1;jj>=0;jj-=1) if (numc[jj]!=0) { 
    if (ce!=nil) if (ix[jj]<ce.count()) { 
      if (argtype(1)==1) $o1=ce.o(ix[jj]) else XO=ce.o(ix[jj]) 
    }
    if (argtype(1)==3) $&1=jj else ii=jj
    iterator_statement
    if (argtype(2)==3) $&2+=1 else i1+=1
  }
}

// stt() iterates over all synapse types
iterator stt () { local ia,ja,ka,ib,jb,im,jm // global ii,jj,kk
  i1=ib=jb=0 im=jm=CTYPi-1 
  if (argtype(1)==0) if ($1>=0) ib=im=$1
  if (argtype(2)==0) if ($2>=0) jb=jm=$2
  if (argtype(4)==3) $&4=0
  for ia=ib,im for ja=jb,jm if (numc[ia]!=0 && numc[ja]!=0) {
    if (ce!=nil) if (ix[ia]<ce.count()) {XO=ce.o(ix[ia]) YO=ce.o(ix[ja])}
    for ka=0,STYPi-1 if (wmat[ia][ja][ka]!=0) {
      if (argtype(3)==3) { $&1=ia $&2=ja $&3=ka } else { ii=ia jj=ja kk=ka }
      iterator_statement
      if (argtype(4)==3) $&4+=1 else i1+=1
    }
  }
}

ixi=ixj=0 // ixi steps through cell#s and ixj is consecutive
for ii=0,CTYPi-1 ixe[ii]=-1
// ctype(#,TYPE)
func ctype () { local ii
  if (numarg()==2) {
    return ($1>=ix[$2] && $1<=ixe[$2]) 
  } else {
    for ii=0,CTYPi-1 if ($1>=ix[ii] && $1<=ixe[ii]) return ii 
    return -1
  }
}
// ixt(type[,col,vec]); ixt(type,&x); ixt(type,col,&x)
iterator ixt () { local ty,col,a,lfl,ixa,i localobj o
  ty=$1 lfl=0
  if (argtype(2)==0) col=$2 else { col=-1 if (argtype(2)==3) lfl=2 }
  if (argtype(3)==1) { o=$o3 vrsz(numc[ty],o,"O") } else if (argtype(3)==3) lfl=3
  i=lfl
  if (col<0) {
    for ({ixa=ix[ty] ixj=0};ixa<=ixe[ty];{ixa+=1 ixj+=1}) { XO=ce.o(ixa)
      if (lfl) $&i=ixa else ixi=ixa
      iterator_statement
    }
  } else {
    if (col>=ncols) printf("ixt ERR: Only %d columns (%d)\n",ncols,col) else {
      a=ix[ty]+col*numc[ty]/ncols
      for ({ixa=a ixj=0};ixa<a+numc[ty]/ncols;{ixa+=1 ixj+=1}) { XO=ce.o(ixa)
        if (XO.col!=col) printf("ixt ERR: XO.col!=col, %d %d\n",XO.col,col)
        if (lfl) $&i=ixa else ixi=ixa
        iterator_statement
      }
    }
  }
  if (o!=nil) stat(o)
}

// for ltr(XO,cvode.netconlist("", "", "")) print XO.precell,XO.postcell,XO.syn
//* synapse linking -- old routines not using NQS
//** geolink(s1,s2,s3) s1 for presyns and s2 for postsyns, connect s3's
// only checks for INTF as a possible pre/post point process
// connects a layer to another assuming 1 dim netgeom
proc geolink() { local i,ii,a,sav,nowrap,noself
  if (numarg()==0) { print "\tgeolink(s1,s2,s3,DEFCON)"
    print "  make connections from s1 to s2 connecting onto s3 type PPs."
    print "  DEFCON struct defines connection"
    return
  }
  // default to yes wrap; yes within-column connect unless $s1==$s2
  nowrap = !$o4.wrap  noself=!$o4.self
  if (strcmp($s1,$s2)==0) $o4.self=0
  tmpobj=new List($s1)
  // object containing a receive/event PP need fflag=1 and intf as name of PP
  if (tmpobj.object(0).fflag) { // presynaptic object flag
    sprint(temp_string_,"ncl.append(new NetCon(XO.intf, YO.%s))",$s3)
  } else {
    sprint(temp_string_,"XO.soma ncl.append(new NetCon(&v(.5), YO.%s))",$s3)
  }
  for ltr (XO,tmpobj,&y) { // presyn list
    $o4.grot(y)
    for lvtr (YO,&x,new List($s2),$o4.scr) { // postsyn list and vector
      if (x==1) { // connect
        print "connecting ",XO," to ",YO
        execute(temp_string_)
      }
    }
  }
  XO=nil YO=nil
}

//** glink() same as geolink but takes lists statt strings
proc glink() { local i,ii,a,sav,nowrap,noself
  if (numarg()==0) { print "\tgeolink(l1,l2,DEFCON)"
    print "  make connections from items in l1 to items in l2."
    print "  DEFCON struct defines connection"
    return
  }
  // default to yes wrap; yes within-column connect unless $s1==$s2
  nowrap = !$o4.wrap  noself=!$o4.self
  // object containing a receive/event PP need fflag=1 and intf as name of PP
  if ($o1.object(0).fflag) { // presynaptic object flag
    sprint(temp_string_,"ncl.append(new NetCon(XO.intf, YO.%s))",$s3)
  } else {
    sprint(temp_string_,"XO.soma ncl.append(new NetCon(&v(.5), YO.%s))",$s3)
  }
  for ltr (XO,$o1,&y) { // presyn list
    $o4.grot(y)
    for lvtr (YO,&x,$o2,$o4.scr) { // postsyn list and vector
      if (x==1) { // connect
        // print "connecting ",XO," to ",YO
        execute(temp_string_)
      }
    }
  }
  XO=nil YO=nil
}

//** synlink(s1,s2,s3,[gmax,del]) s1 for presyns and s2 for postsyns, connect s3's
//  eg synlink("PYR","INH","AMPA")
// provides full connectivity
proc synlink() { local gm,dl
  if (numarg()==0) { print "\tsynlink(s1,s2,s3)"
    print "  make connections from all of type s1 to one of type s2 "
    print "     only connecting onto s3 type PPs."
    print "  Matching done with regexps."
    return
  }
  if (numarg()==5) { gm=$4 dl=$5 } else { gm=0 dl=1 }
  tmplist = new List($s3) // list of postsyn point processes 
  for ltr (XO,new List($s1)) { // presyn list
    for ltr (YO,tmplist) {
      if (pp_loc(YO,$s2,syn2)) { 
        sfunc.head(syn2,"\\.",syn2)
        sprint(syn1,"%s",XO)
        if (strcmp(syn1,syn2)!=0) { // don't allow self connects
          print "connecting ",XO," to ",syn2
          XO.soma ncl.append(new NetCon(&v(.5), YO, thresh, gm, dl))
        }
      }
    }
  }
}
      
// run NQS.mo(1) first to set up PRIDv ... vectors

//** simple netconnect routines: netconn(), netgen()
proc netconn () { local pre,post,syn,pri,poi
  pre=$1 post=$2 
  if (numarg()==5) { syn=$3 pri=$4 poi=$5 } else { syn=0 pri=$3 poi=$4 }
  cells[pre].object(pri).conn(cells[post].object(poi).syns[syn])
}

proc netgen () { ncl.append(new NetCon($o1, $o2)) }

//** syncopy() copies from one set of syns to another for colocalization
proc syncopy () { 
  if (numarg()==0) { printf("syncopy(to_type,[]): copy from type, post/type, pre/post/type\n") return }
  sprint(temp_string_,"XO.precell.soma ncl.append(new NetCon(&v(.5),XO.postcell.%s",$s1)
  if (numarg()==1) tmplist = cvode.netconlist("", "" , $s2)
  if (numarg()==2) tmplist = cvode.netconlist("", $s2, $s3)
  if (numarg()==3) tmplist = cvode.netconlist($s2, $s3, $s4)
  for ltr(XO,tmplist) execute(temp_string_)
}

//** syn1to1(PRE,POST,SYN)
proc syn1to1 () { local pre,post,syn
  pre=$1 post=$2
  if (numarg()==3) syn=$3 else syn=0
  if (cells[pre].count !=cells[post].count) { 
    printf("\tERROR: 1-to-1 connectivity requires same # of %s (=%d) and %s (=%d).\n",\
           CTYP.object(pre).s,cells[pre].count,CTYP.object(post).s,cells[post].count) return }
  for ltr2(XO,YO,cells[pre],cells[post]) {
    printf("SRC: %s -> TRG: %s (%s)\n",XO,YO,YO.syns(syn))
    XO.conn(YO.syns[syn])
  }
}

//* synconn() synapse linking -- NQS routines
//** synconn(preid,postid,pre#,post#,%div)
// eg synconn(PYRi,PYRi,AMPAi,pmat[PYR][PYR])
// provides % connectivity based on C/pre==%div==%conv==D/post
// S==C*post==D*pre %conv=S/(post*pre)=C/pre=D/post
objref convec,convec1,convec2  // vectors to count how much div is from each nrn
convec = new Vector(1e3)
convec1 = convec.c // a scratch vector
convec2 = convec.c // copy of convec to blot out self-connect and redundent connects

proc synconn () { local preid,posid,pdiv,con,div,ii,jj,prn,pon,targ,sz,styp1,styp2
  if (numarg()==0) { print "\tsynconn(preid,postid,prn,pon,pdiv)" return }
  preid=$1 posid=$2 prn=$3 pon=$4 pdiv=$5
  CODEv=sp.v[sp.fi("CODE")] PRv=sp.v[sp.fi("PR")] POv=sp.v[sp.fi("PO")]
  sz=PRv.size
  if (pdiv==1) {
    if (preid==posid) div=pon-1 else div=pon
    con=div
  } else {
    con=int(pdiv*prn+1) div=int(pdiv*pon)
  }
  if (isobj(CTYP,"List")) {
    printf("%s->%s:\tdiv=%d,conv=%d (%d syns)\n",CTYP.object(preid).s,CTYP.object(posid).s,div,con,prn*div)
  } else {
    printf("%d->%d:\tdiv=%d,conv=%d (%d syns)\n",preid,posid,div,con,prn*div)
  }
  if (prn*div==0) return
  sp.pad(sz+prn*div)
  if (pdiv==1) {
    convec.indgen(0,pon-1,1)
    for ii=0,prn-1 {
      if (preid==posid) {convec.indgen(0,pon-1,1) convec.remove(ii)}
      POv.copy(convec,sz+ii*div)
      PRv.fill(ii,sz+ii*div,sz+(ii+1)*div-1)
    }
  } else {
    convec.resize(pon) convec.fill(0) // counter for convergence
    for ii=0,prn-1 { // sources
      convec2.copy(convec) 
      if (preid==posid) convec2.x[ii]=1e10 // block self-connect
      for jj=1,div POv.set(sz+ii*div+jj-1,pickpost(pon,con)) // pick desired target
      PRv.fill(ii,sz+ii*div,sz+(ii+1)*div-1)
    }
  }
  CODEv.fill(mkcodf(preid,posid,0,0,0),sz,PRv.size-1)
}

proc syncnn () { local preb,pree,posb,pose,pdiv,con,div,ii,jj,prn,pon,targ,sz,styp1,styp2
  if (numarg()==0) { print "\tsyncnn(preb,pree,posb,pose,pdiv)" return }
  preb=$1 pree=$2 posb=$3 pose=$4  pdiv=$5
  CODEv=sp.v[sp.fi("CODE")] PRv=sp.v[sp.fi("PR")] POv=sp.v[sp.fi("PO")]
  sz=PRv.size
  prn=pree-preb+1  pon=pose-posb+1
  if (pdiv==1) {
    if (preid==posid) div=pon-1 else div=pon
    con=div
  } else {
    con=int(pdiv*prn+1) div=int(pdiv*pon)
  }
  if (isobj(CTYP,"List")) {
    printf("%s->%s:\tdiv=%d,conv=%d (%d syns)\n",CTYP.object(preid).s,CTYP.object(posid).s,div,con,prn*div)
  } else {
    printf("%d->%d:\tdiv=%d,conv=%d (%d syns)\n",preid,posid,div,con,prn*div)
  }
  if (prn*div==0) return
  sp.pad(sz+prn*div)
  if (pdiv==1) {
    convec.indgen(0,pon-1,1)
    for ii=0,prn-1 {
      if (preid==posid) {convec.indgen(0,pon-1,1) convec.remove(ii)}
      POv.copy(convec,sz+ii*div)
      PRv.fill(ii,sz+ii*div,sz+(ii+1)*div-1)
    }
  } else {
    convec.resize(pon) convec.fill(0) // counter for convergence
    for ii=0,prn-1 { // sources
      convec2.copy(convec) 
      if (preid==posid) convec2.x[ii]=1e10 // block self-connect
      for jj=1,div POv.set(sz+ii*div+jj-1,pickpost(pon,con)) // pick desired target
      PRv.fill(ii,sz+ii*div,sz+(ii+1)*div-1)
    }
  }
  PRv.add(preb) POv.add(posb)
  CODEv.fill(mkcodf(preid,posid,0,0,0),sz,PRv.size-1)
}

//** synconn2() uses elimrepeats() and shuffle methods
proc synconn2 () { local preid,posid,pdiv,con,div,ii,jj,prn,pon
  if (numarg()==0) { print "\tsynconn2(preid,postid,prn,pon,pdiv)" return }
  preid=$2 posid=$3 prn=$4 pon=$5 pdiv=$6
  $o1.clear()
  PRv=$o1.v[$o1.fi("PR")] POv=$o1.v[$o1.fi("PO")]
  con=int(pdiv*prn+1) div=int(pdiv*pon)
  if (prn*div==0) return
  printf("%s->%s:\tdiv=%d,conv=%d (%d syns)\n",CTYP.object(preid).s,CTYP.object(posid).s,div,con,prn*div)
  if (pdiv==1) {
    $o1.pad(prn*div)
    convec.indgen(0,pon-1,1)
    for ii=0,prn-1 {
      POv.copy(convec,ii*div)
      PRv.fill(ii,ii*div,(ii+1)*div-1)
    }
  } else {
    $o1.pad(1.5*prn*div)
    rdm.discunif(0,prn-1)  PRv.setrand(rdm)
    rdm.discunif(0,pon-1)  POv.setrand(rdm)
    $o1.elimrepeats("PR","PO")
    $o1.shuffle()
    $o1.pad(prn*div)
  }
  $o1.fill("CODE",1,preid) $o1.fill("CODE",2,posid)
}

//** synconn3() doesn't worry about eliminating repeats
proc synconn3 () { local preid,posid,pdiv,con,div,ii,jj,prn,pon
  if (numarg()==0) { print "\tsynconn2(preid,postid,prn,pon,pdiv)" return }
  preid=$2 posid=$3 prn=$4 pon=$5 pdiv=$6
  $o1.clear()
  PRv=$o1.v[$o1.fi("PR")] POv=$o1.v[$o1.fi("PO")]
  con=int(pdiv*prn+1) div=int(pdiv*pon)
  if (prn*div==0) return
  printf("%s->%s:\tdiv=%d,conv=%d (%d syns)\n",CTYP.object(preid).s,CTYP.object(posid).s,div,con,prn*div)
  $o1.pad(prn*div)
  rdm.discunif(0,prn-1)  PRv.setrand(rdm)
  rdm.discunif(0,pon-1)  POv.setrand(rdm)
  $o1.fill("CODE",1,preid) $o1.fill("CODE",2,posid)
}

//*** pickpost() tries to reduce divergence variance
// pickpost(postlist,maxcon,YO)
// maxcon == -1 means to ignore convergence
// MUST do convec.resize() and convec.fill(0) before using
func pickpost () { local ran, maxcon, maxpo, min, indx
  maxcon = $2  // max convergence to be allowed
  maxpo = $1   // number of postsyn choices
  min = convec2.min  // convec should start all 0's
  if (min >= maxcon) { // all full up
    printf("Pickpost full WARNING: %d %d\n",min,maxcon) vlk(convec2) }
  convec1.indvwhere(convec2,"==",min)  // look for all the smallest to fill up first
  inx = convec1.x[rdm.discunif(0,convec1.size-1)]
  convec.x[inx]+=1
  convec2.x[inx]=1e10 // block from reconnecting here
  return inx
}

//** smap()
// excitatory cells project to AMPA and inhibitory cells to GABA
// assumes PRIDv ... defined by sp.mo(1)
objref pro,poo // pre and post pointers
func smap () { local ct,sy,conv,prdx,podx,prx,pox,delx,w0x,w1x localobj nc,ty
  ty=new Union()
  ct=cp.fi("PRID") sy=cp.fi("STYP")
  sp.resize("NC")
  sp.odec("NC")
  sp.pad()
  sp.mo(1)
  for ii=0,PRv.size-1 {
    uncodf(CODEv.x[ii],&prdx,&podx) prx=PRv.x[ii] pox=POv.x[ii]
    delx=DELv.x[ii] w0x=WT0v.x[ii] w1x=WT1v.x[ii]
    pro=c[prdx].object(prx) poo=c[podx].object(pox)
    NCv.x[ii]=sp.fcdo.append(nc=smap1(prdx))-1
    for kk=0,nc.wcnt-1 nc.weight[kk]=0
    x=cp.fetch(ct,prdx,sy)
    syntyp(x,ty)
    nc.weight[ty.x]=w0x
    if (ty.x[1]>-1) nc.weight[ty.x[1]]=w1x
    nc.delay=delx
  }
  return ii
}

//*** smap1(SYN#) poo has postsyn and pro has pre
obfunc smap1 () { localobj si
  if (poo.fflag) { YO=poo 
  } else {
    YO=poo.po[$1]
    if (isobj(YO,"List")) {
      snsr($1,poo)
      YO=YO.object(YO.count-1)
    }
  }
  if (pro.fflag) {            si=new NetCon(pro,    YO)
  } else {           pro.soma si=new NetCon(&v(0.5),YO) }
  return si
}

//*** snsr() handles situation where multiple PPs must be hung onto postsyn objref
proc snsr () {
  printf("PROBLEM: replace snsr() which uses an execute\n")
  if (isobj($o2.po[$1],"List")) {
    sprint(tstr,"%s.soma %s.po[%d].append(new %s(0.5))",$o2,$o2,$1,STYP.object($1).s)
  } else {
    sprint(tstr,"%s.soma %s.po[%d] = new %s(0.5)",$o2,$o2,$1,STYP.object($1).s)
  }
  execute(tstr)
}

//** umbrella(preid,postid,max,width[,shift])
// set lateral projections from pre to post out to width
// sets direct connection (jj==ii) iff preid!=posid
proc umbrella () { local ii,jj,preid,posid,width,sh,max,sz,styp1,styp2
  preid=$1 posid=$2 max=$3 width=$4 
  if (numarg()>=5) sh=$5 else sh=0
  sz=PRv.size styp1=styp2=-1
  if (width) { 
    for ii=0,max-1 for jj=ii-width+sh,ii+width+sh {
      if (jj>=0 && jj<max && (preid!=posid || jj!=ii)) {
        PRv.append(ii) POv.append(jj) 
        PRIDv.append(preid) POIDv.append(posid)
      }
    }
  } else { // just within column connections
    for ii=0,max-1 { PRv.append(ii) POv.append(ii) PRIDv.append(preid) POIDv.append(posid) }
  }
  sp.pad()
  // WID0v.fill(styp1,sz,PRv.size-1)
  // WID1v.fill(styp2,sz,PRv.size-1)
}

//** umbflow(preid,postid,max,width[,shift])
// like umbrella but does 'flow' boundary conditions -- reflects back on sides
proc umbflow () { local ii,jj,ja,preid,posid,width,sh,max,sz,styp1,styp2
  preid=$1 posid=$2 max=$3 width=$4 
  if (numarg()>=5) sh=$5 else sh=0
  sz=PRv.size styp1=styp2=-1
  if (width) { 
    for ii=0,max-1 for jj=ii-width+sh,ii+width+sh {
      if (preid!=posid || jj!=ii) { ja=jj
        if (jj<0) ja=-jj
        if (jj>max-1) ja=2*max-jj-2
        PRv.append(ja) POv.append(ii) 
        PRIDv.append(preid) POIDv.append(posid)
      }
    }
  } else { // just within column connections
    for ii=0,max-1 { PRv.append(ii) POv.append(ii) PRIDv.append(preid) POIDv.append(posid) }
  }
  sp.pad()
  // WID0v.fill(styp1,sz,PRv.size-1)
  // WID1v.fill(styp2,sz,PRv.size-1)
}

//** nqdiv(PRID,POID,PR,PO0[,PO1,...]) 
// nqdiv(PRID,POID,PR,POBEG,..,POEND)
proc nqdiv () { local i,beg,end
  if (numarg()==0) { 
    print "nqdiv(PRID,POID,PR,PO0[,PO1,...])\nnqdiv(PRID,POID,PR,POBEG,\"..\",POEND)" return }
  if (numarg()==6 && argtype(5)==2) {
    beg=$4 end=$6
    for i=beg,end sp.append("CODE",EQU1,$1,"CODE",EQU2,$2,"PR",$3,"PO",i)
  } else for i=4,numarg() sp.append("CODE",EQU1,$1,"CODE",EQU2,$2,"PR",$3,"PO",$i)
  sp.pad()
}

//** nqconv(POID,PRID,PO,PR0[,PR1,...]) 
// nqconv(POID,PRID,PO,PRBEG,..,PREND)
proc nqconv () { local i,beg,end
  if (numarg()==0) { 
    print "nqconv(POID,PRID,PO,PR0[,PR1,...])\nnqconv(POID,PRID,PO,PRBEG,\"..\",PREND)" return }
  if (numarg()==6 && argtype(5)==2) {
    beg=$4 end=$6
    for i=beg,end sp.append("PRID",$2,"POID",$1,"PR",i,"PO",$3)
  } else for i=4,numarg() sp.append("PRID",$2,"POID",$1,"PR",$i,"PO",$3)
  sp.pad()
}

//** proc nq1to1(PRID,POID,num)
proc nq1to1 () { umbrella($1,$2,$3,0) }

//* direct weight setting routines
//** scale gmax
proc synscgmax () {
  for ltr(XO,cvode.netconlist("" , "", $s1)) { XO.weight *= $2 }
}
//** set/get gmax
proc synsetgmax () {
  if (numarg()==0) { print "synsetgmax(pre,post,type,wt)"
  } else if (numarg()==2) {for ltr(XO,cvode.netconlist("" , "", $s1)) XO.weight=$2
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist("" , $s1, $s2)) XO.weight=$3
  } else if (numarg()==4) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) XO.weight=$4
  } else { print "ERROR: too many args" }
  if (i1==0) printf("WARNING: nothing set in synsetgmax(%s ...)\n",$s1)
}

// synnormgmax() -- like synsetgmax except normalizes the wt by convergence
proc synnormgmax () {
  for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) {
    XO.weight=$4/cvode.netconlist($s1,XO.postcell,$s3).count
  }
}

proc syngetgmax () {
  if (numarg()==1) {
    if (strcmp($s1,"help")==0) { printf("type,post/type,pre/post/type\n") } else {
      for ltr(XO,cvode.netconlist("" , "", $s1)) printf("%g ",XO.weight) }
  } else if (numarg()==2) {for ltr(XO,cvode.netconlist("" , $s1, $s2)) printf("%g ",XO.weight)
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) printf("%g ",XO.weight)
  } else for ltr(XO,cvode.netconlist("","","")) printf("%g ",XO.weight)
  print ""
}
 
//** set/get delay
proc synsetdel () {
  if (numarg()==2) {for ltr(XO,cvode.netconlist("" , "", $s1)) XO.delay=$2
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist("" , $s1, $s2)) XO.delay=$3
  } else if (numarg()==4) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) XO.delay=$4
  } else { print "ERROR: too many args" }
  if (i1==0) printf("WARNING: nothing set in synsetdel(%s ...)\n",$s1)
}

proc syngetdel () {
  if (numarg()==1) {for ltr(XO,cvode.netconlist("" , "", $s1)) printf("%g ",XO.delay)
  } else if (numarg()==2) {for ltr(XO,cvode.netconlist("" , $s1, $s2)) printf("%g ",XO.delay)
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) printf("%g ",XO.delay)
  } else for ltr(XO,cvode.netconlist("","","")) printf("%g ",XO.delay)
  print ""
}
 
//** set/get thresh
proc synsetthr () {
  if (numarg()==1) {for ltr(XO,cvode.netconlist("" , "", "")) XO.threshold=$1
  } else if (numarg()==2) {for ltr(XO,cvode.netconlist($s1 , "", "")) XO.threshold=$2
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist($s1, "", $s2)) XO.threshold=$3
  } else if (numarg()==4) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) XO.threshold=$4
  } else { print "ERROR: too many args" }
  if (i1==0) printf("WARNING: nothing set in synsetthr(%s ...)\n",$s1)
}

proc syngetthr () {
  if (numarg()==1) {for ltr(XO,cvode.netconlist($s1 , "", "")) printf("%g ",XO.threshold)
  } else if (numarg()==2) {for ltr(XO,cvode.netconlist($s1, "", $s2)) printf("%g ",XO.threshold)
  } else if (numarg()==3) {for ltr(XO,cvode.netconlist($s1 , $s2, $s3)) printf("%g ",XO.threshold)
  } else for ltr(XO,cvode.netconlist("","","")) printf("%g ",XO.threshold)
  print ""
}
 
//** synremove: remove post, pre/post, pre/post/type; synshow
proc synremove () {
  if (numarg()==0) { printf("synremove: remove post, pre/post, pre/post/type\n") return }
  if (numarg()==1) tmplist = cvode.netconlist("", $s1 , "")
  if (numarg()==2) tmplist = cvode.netconlist("", $s1, $s2)
  if (numarg()==3) tmplist = cvode.netconlist($s1 , $s2, $s3)
  if (tmplist.count>0) for ltr(XO,tmplist) ncl.remove(ncl.index(XO))
  tmplist = nil // need to remove these references too
  if (numarg()==1) tmplist = cvode.netconlist("", $s1 , "")
  if (numarg()==2) tmplist = cvode.netconlist("", $s1, $s2)
  if (numarg()==3) tmplist = cvode.netconlist($s1 , $s2, $s3)
  if (tmplist.count>0) for ltr(XO,tmplist) printf("ERROR: %s removed from ncl but still exists\n",XO)
  tmplist = nil
}
  
proc synshow () { 
  sprint(temp_string_,"x=XO.%s",$s2)
  for ltr(XO,cvode.netconlist("" , "", $s1)) { execute(temp_string_) printf("%g ",x) }
  print ""
}

//* weight setting routines using NQS
//** stwt(PREID,POSTID,WT0[,WT1,norm])
// stwtnorm() causes dividing by individual convergence values
proc stwt () { local w0,w1
  if (sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2) ==0) {
    printf("WARNING NO CONNECTS FROM %d TO %d\n",$1,$2) return }
  w0=$3
  if (numarg()>=5) {  w0=$3/$5 w1=$4/$5 } else if (numarg()>=4) { w1=$4 }
  if (numarg()>=4) sp.fillin("WT0",w0,"WT1",w1) else sp.fillin("WT0",w0)
}

//** strwt(PREID,POSTID,WT0,WT1,psdev[,norm])
proc strwt () { local w0,w1,psdev
  cnt=sp.select("CODE",EQU1,$1,"CODE",EQU2,$2)
  if (cnt==0) {printf("WARNING NO CONNECTS FROM %d TO %d\n",$1,$2) return }
  if (numarg()>=6) {  w0=$3/$6 w1=$4/$6 } else {  w0=$3 w1=$4 }
  psdev=$5
  rdm.lognormal(w0,psdev*psdev*w0*w0)
  sp.out.v[sp.fi("WT0")].setrand(rdm)
  if (w1!=0) {
    rdm.lognormal(w1,psdev*psdev*w1*w1)
    sp.out.v[sp.fi("WT1")].setrand(rdm)
  }
  sp.delect()
}

//** strwt2(NQS,WT0,WT1,psdev[,norm])
proc strwt2 () { local w0,w1,psdev
  if (numarg()>=5) {  w0=$2/$5 w1=$3/$5 } else {  w0=$2 w1=$3 }
  psdev=$4
  rdm.lognormal(w0,psdev*psdev*w0*w0)
  $o1.v[$o1.fi("WT0")].setrand(rdm)
  if (w1!=0) {
    rdm.lognormal(w1,psdev*psdev*w1*w1)
    $o1.v[$o1.fi("WT1")].setrand(rdm)
  }
}

//** strdel(PREID,POSTID,del,psdev)
proc strdel () { local del0,psdev
  if (numarg()==4) {
    cnt=sp.select("CODE",EQU1,$1,"CODE",EQU2,$2)
    if (cnt==0) {printf("WARNING NO CONNECTS FROM %d TO %d\n",$1,$2) return }
    del0=$3 psdev=$4
    rdm.lognormal(del0,psdev*psdev*del0*del0)
    sp.out.v[sp.fi("DEL")].setrand(rdm)
    sp.delect()
  } else { del0=$2 psdev=$3
    rdm.lognormal(del0,psdev*psdev*del0*del0)
    $o1.v[sp.fi("DEL")].setrand(rdm)
  }
}

//** clrwt(PRID,POID,%clr)
proc clrwt () { local n
  n=round($3*sp.select("CODE",EQU1,$1,"CODE",EQU2,$2))
  sp.out.v[sp.fi("WT0")].fill(0,0,n)
  sp.delect()
}

//** chksp() confirm correspondance between sp and ncl -- defunct if using multiple ncl
proc chksp () { local prid,poid,pr,po
  for ltr(XO,ncl,&x) {
    prid=sp.v[0].x[x] poid=sp.v[1].x[x] pr=sp.v[2].x[x] po=sp.v[3].x[x] 
    if (XO.pre.id!=pr || XO.pre.type!=prid ||  XO.syn.id!=po || XO.syn.type!=poid) {
      printf("%d %d %d %d %d %d %d %d\n",XO.pre.id,pr,XO.pre.type,prid,XO.syn.id,po,XO.syn.type,poid)  }}
}

//** swtmap() -- redund code to permit it to run relatively fast
// no longer overloaded to set delay, delmap() or to clr clrwts()
proc swtmap () { local ii,ct,sy,conv,prx,pox,delx,w0x,w1x localobj nc,ty
  ty=new Union()
  ct=cp.fi("PRID") sy=cp.fi("STYP")
  if (PRv.size!=$o1.count) { 
    printf("swtmap ERR: size discrepency: %d vs %d\n",PRv.size,$o1.count)
    return }
  for ii=0,PRv.size-1 {
    prdx=PRIDv.x[ii] podx=POIDv.x[ii] prx=PRv.x[ii] pox=POv.x[ii]
    delx=DELv.x[ii] w0x=WT0v.x[ii] w1x=WT1v.x[ii]
    nc=$o1.object(ii)
    poo=nc.syn
    x=cp.fetch(ct,prdx,sy)
    syntyp(x,ty)
    nc.weight[ty.x]=w0x
    if (ty.x[1]>-1) nc.weight[ty.x[1]]=w1x
    nc.delay=delx
  }
}

//** wmul(PRID,POID,SYNID,MUL) multiplies set of syns by a factor
// resets from weights stored in sp
for ii=0,CTYPi-1 for jj=0,CTYPi-1 for kk=0,STYPi-1 wd0[ii][jj][kk]=1
proc wmul () { local ii,pr1,po1,sy1,w1,wd,sy2
  pr1=$1 po1=$2 sy1=$3 w1=$4
  if (wd0[pr1][po1][sy1]==0) {
    if (w1!=0) printf("\nWARNING: %s->%s(%s) set to 0 can't reset\n",\
                      CTYP.o(pr1).s,CTYP.o(po1).s,STYP.o(sy1).s)
    return
  }
  wd=w1/wd0[pr1][po1][sy1]
  if (wd!=1) {
    XO=ncl[pr1][po1][0]
    for ii=0,XO.count-1 XO.o(ii).weight[sy1]*=wd
    wd0[pr1][po1][sy1]=w1
  } else printf(".")
} 

//** wset(PRID,POID,SYNID,WT) resets weights to gaussian dist similar to strwt()
// uses sp to determine how many weights are needed
proc wset () { local num,prx,pox,sox,wet,a,err,psdev localobj v1
  prx=$1 pox=$2 sox=$3 wet=$4 
  if (numarg()==5) psdev=$5*$5 else psdev=0
  err=0 a=allocvecs(v1)
  if (cp.fetch("PRID",prx,"STYP")==EX && !(sox==AM || sox==NM)) err=1
  if (cp.fetch("PRID",prx,"STYP")==IX && !(sox==GA || sox==GB)) err=1
  if (err) { printf("cell/syn discrepancy: %d -> %d\n",prx,sox) return }
  if (Incol==2) num=sp.select(-1,"CODE",EQU1,prx,"CODE",EQU2,pox) else {
    num=sp.select(-1,"CODE",EQU1,prx,"CODE",EQU2,pox,"CODE",EQU3,Incol)
  }
  v1.resize(num)
  if (psdev) { 
    rdm.normal(wet,psdev*wet*wet) 
    v1.setrand(rdm)
  } else v1.fill(wet)
  vwtmap(v1,sp,sox)
  dealloc(a)
}

//** wbim(PRID,POID,SYNID,PSDEV,%A,WTA,%B,WTB,...) bimodal weight setting routines
// uses sp to determine how many weights are needed
proc wbim () { local fsz,i,prx,pox,sox,wet,a,err,psdev,pcw,ptot localobj v1,v2
  prx=$1 pox=$2 sox=$3 psdev=$4*$4 pcw=$5 wet=$6 ptot=0
  err=0 a=allocvecs(v1,v2)
  if (cp.fetch("PRID",prx,"STYP")==EX && !(sox==AM || sox==NM)) err=1
  if (cp.fetch("PRID",prx,"STYP")==IX && !(sox==GA || sox==GB)) err=1
  if (err) { printf("cell/syn discrepancy: %d -> %d\n",prx,sox) return }
  fsz=sp.select(-1,"CODE",EQU1,prx,"CODE",EQU2,pox)
  for i=5,numarg() {
    pcw=$i i+=1 wet=$i // 
    if (pcw*fsz != int(pcw*fsz)){
      printf("wbim WARNING: %d->%d non-int %gx%g=%g\n",prx,sox,pcw,fsz,pcw*fsz) }
    v2.resize(int(pcw*fsz))
    rdm.normal(wet,psdev*wet*wet) 
    v2.setrand(rdm)
    v1.append(v2)
    ptot+=pcw // keep track of percent filled
  }
  if (ptot!=1) printf("wbim WARNING: %d->%d only %g filled\n",prx,pox,ptot)
  if (v1.size!=fsz) {
    printf("wbim WARNING: %d->%d size error %d->%d\n",prx,sox,v1.size,fsz)
    v1.resize(fsz) }
  vwtmap(v1,sp,sox)
  dealloc(a)
}

//** vwtmap(VEC,NQS) -- copy weights from a vector onto ncl
//   vwtmap(NQS,VEC) -- copy weights from ncl into vector
func vwtmap () { local wix
  if (numarg()==3) wix=$3 else wix=0
  if (isojt($o2,sp)) { // copy from vector to sp.fcdo weights 
    if ($o1.size!=$o2.ind.size){
      printf("vwtmap ERR: size discrepency: %d vs %d\n",$o1.size,$o2.ind.size)
      return 0}
    for ii=0,$o1.size-1 $o2.fcdo.o[$o2.ind.x[ii]].weight[wix]=$o1.x[ii]
    return $o1.size
  } else { // copy from ncl weights to vector
    if ($o1.ind.size!=$o2.size){printf("vwtmap WARNING: resizing vec\n") $o2.resize($o1.ind.size)}
    for ii=0,$o2.size-1 $o2.x[ii]=$o1.o[$o1.ind.x[ii]].weight[wix]
    return $o2.size
  }
}

//** wtstat(WIX) -- check id of weights directly from the ncs
// assume that have already been selected in sp.out
proc wtstat () { local wix,a,sz,fnc localobj v1
  a=allocvecs(v1)
  if (!eqojt(sp.cob,sp.out)) print "WARNING: no sp.select() done"
  sz=sp.size(1) fnc=sp.fi("NC")
  v1.resize(sz) v1.resize(0)
  if (numarg()==1) wix=$1 else wix=0
  // for sp.qt(XO,"NC") v1.append(XO.weight[wix]) // 5x slower
  for ii=0,sz-1 v1.append(sp.fcdo.o(sp.out.v[fnc].x[ii]).weight[wix])
  stat(v1)
  dealloc(a)
}

// need to set delays 
proc delmap () { local ii,deli,nc0,nc1
  nc0=$o1.fi("NC0") nc1=$o1.fi("NC1","NOERR") deli=$o1.fi("DEL")
  for ii=0,$o1.v.size-1 {
    ncl.object($o1.v[nc0].x[ii]).delay=$o1.v[deli].x[ii]
    if (nc1!=-1) if ($o1.v[nc1].x[ii]!=-1) {
      ncl.object($o1.v[nc1].x[ii]).delay=$o1.v[deli].x[ii] }
  }
}

// clrwtmap
proc clrwts () { local ii,jj,nc0,wt0,deli,nc1,wt1,typ,sy2,wid0,wid1,clr
  nc0=$o1.fi("NC0") wt0=$o1.fi("WT0") deli=$o1.fi("DEL") 
  wid0=$o1.fi("WID0") typ=$o1.fi("TYPE")
  nc1=$o1.fi("NC1","NOERR") wt1=$o1.fi("WT1","NOERR") wid1=$o1.fi("WID1","NOERR")
  tobj=ncl.object($o1.v[nc0].x[ii])
  for ii=0,$o1.v.size-1 { 
    tobj=ncl.object($o1.v[nc0].x[ii])
    for jj=0,tobj.wcnt-1 tobj.weight[jj]=0 // clear
    tobj.delay=1 tobj.threshold=0 
    if (nc1!=-1) if ($o1.v[nc1].x[ii]!=-1) for jj=0,tobj.wcnt-1 tobj.weight[jj]=0 // clear
  }
}

// swtchk(sp,"WT0","NC0") compare weights to sp weights
func swtchk () { local ii,jj,nc0,wt0,n
  nc0=$o1.fi($s3) wt0=$o1.fi($s2)  n=0
  for ii=0,$o1.v.size-1 { 
    if ($o1.v[nc0].x[ii]==-1) continue
    tobj=ncl.object($o1.v[nc0].x[ii])
    if ($o1.v[wt0].x[ii]==tobj.weight) n+=1 else {
      printf("Mismatch %d: %g %g\n",ii,$o1.v[wt0].x[ii],tobj.weight) }
  }
  tobj=nil
  return n/ii
}

// synchk() look for internal consistency
proc synchk () { 
  for ltr(XO,cvode.netconlist("","","")) if (isassigned(XO.postcell) && XO.weight<0) {
    printf("Error for %s with wt %g\n",XO.syn,XO.weight)
  }
}

//** getwt() useful in eg sp.spr("<NC0>.c.apply('getwt')")
func getwt () { return NetCon[$1].weight }

spnormflag = 0 // set this flag to normalize weights by number of projections
//** actumb(PRID,POID,WIDTH) -- clears umbrella and then set sp.ind to chosen umbrella
proc actumb () { local width,flag,divisor
  width=$3
  if (width==-1) flag=1 else flag=0
  sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2) 
  if (! flag) { // width=-1 is flag for full connection
    if (sp.fi("WT1")!=-1) sp.fillin("WT0",0,"WT1",0) else sp.fillin("WT0",0)
    sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2,"DIST","()",-width,width)
  }
  if (! spnormflag) { // just set the values
    if (sp.fi("WT1")!=-1 && numarg()==5) sp.fillin("WT0",$4,"WT1",$5) else sp.fillin("WT0",$4)
  } else { // need to calculate convergence for individual cells here
  }
}

//** inline(PRID,POID,WT) sets the in-column weights
proc inline () { local a
  sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2,"PR",EQV,"PO")
  if (numarg()==4) sp.fillin("WT0",$3,"WT1",$4) else sp.fillin("WT0",$3)
}

// stwtnorm(PREID,POSTID,WT0[,WT1,CONVVEC])
func stwtnorm () { local cv,a,b,ret
  a=b=allocvecs(2) b+=1
  if (sp.select("CODE",EQU1,$1,"CODE",EQU2,$2)==0) {
    printf("WARNING NO CONNECTS FROM %d TO %d\n",$1,$2) return 0 }
  sp.uniq("PO") // only retains unique values of PO => cv>0
  mso[a].copy(sp.out.v[sp.fi("PO")])
  for vtr(&x,mso[a]) {
    cv=sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2,"PO",x) // do select without copy to out
    mso[b].append(cv) 
    if (cv>0) {
      if (numarg()==4) sp.fillin("WT0",$3/cv,"WT1",$4/cv) else sp.fillin("WT0",$3/cv)
    }
  }
  ret=mso[b].mean // mean convergence
  if (numarg()==4) if (argtype(4)==1) $o4.copy(mso[b])
  if (numarg()==5) if (argtype(5)==1) $o5.copy(mso[b])
  dealloc(a)
  return ret
}

// wtnorm(PREID,POSTID) -- normalize the values according to convergence
func wtnorm () { local cv,a,b,ret,wt0,wt1
  wt0=sp.fi("WT0") wt1=sp.fi("WT1") 
  if (sp.select(-1,"CODE",EQU1,$1,"CODE",EQU2,$2)==0) {
    printf("WARNING NO CONNECTS FROM %d TO %d\n",$1,$2) return }
  a=b=allocvecs(2) b+=1
  sp.uniq("PO") // only retains unique values of PO => cv>0
  mso[a].copy(sp.out.v[sp.fi("PO")])
  for vtr(&x,mso[a]) {
    cv=sp.select("CODE",EQU1,$1,"CODE",EQU2,$2,"PO",x)
    mso[b].append(cv) 
    if (cv>0) {
      sp.out.v[wt0].div(cv)
      if (wt1>-1) sp.out.v[wt1].div(cv)
      sp.delect()
    }
  }
  ret=mso[b].mean // mean convergence
  dealloc(a)
  return ret
}



proc setdist () { sp.spr("DIST","<PR>.c.sub(<PO>).apply('abs')") }

proc stwtvar() { local a,idr
  if (numarg()==5) idr=1 else idr=0
  a=allocvecs(1)
  mso[a].resize(sp.ind.size)
  rdm.uniform($3*(1-$4),$3*(1+$4))
  mso[a].setrand(rdm)
  if (spnormflag) mso[a].div(sp.ind.size)
  sp.wt[idr].indset(sp.ind,mso[a])
  dealloc(a)
}

//** snc() -- set the actual netcons to the params in sp
DEL=DELD=1
proc snc () { local ii
  for ii=0,sp.size-1 {
    nc[sp.nc.x[ii]].threshold=0
    nc[sp.nc.x[ii]].weight=sp.wt.x[ii]
    nc[sp.nc.x[ii]].delay=DEL+DELD*sp.dist.x[ii]
    if (sp.nc[1].x[ii]>-1) {
      nc[sp.nc[1].x[ii]].weight=sp.wt[1].x[ii]
      nc[sp.nc[1].x[ii]].delay=DEL+DELD*sp.dist.x[ii]
    }
  }
}
//** snci() -- take a vec of indices (eg sp.ind) as arg
proc snci () { local ii,jj
  for jj=0,$o1.size-1 {
    ii=$o1.x[jj]
    nc[sp.nc.x[ii]].weight=sp.wt.x[ii]
    nc[sp.nc.x[ii]].delay=DEL+DELD*sp.dist.x[ii]
    if (sp.nc[1].x[ii]>-1) {
      nc[sp.nc[1].x[ii]].weight=sp.wt[1].x[ii]
      nc[sp.nc[1].x[ii]].delay=DEL+DELD*sp.dist.x[ii]
    }
  }
}

//* informational procs
//** proc print_pp_location(PP), from doc/html/refman/nocmodl.html
proc print_pp_location() { local x //arg1 must be a point process
   x = $o1.get_loc()
   sectionname(section)
   printf("%s located at %s(%g)\n", $o1, section, x)
   pop_section()
}
//** pp_loc(PP,LOC,SCRATCH) returns true if point process PP is located in LOC (regexp match)
func pp_loc () { local x //arg1 must be a point process
   x = 0
   $o1.get_loc()
   if (numarg()==3) { sectionname($s3) }
   ifsec $s2 { x = 1 }
   pop_section()
   return x
}
 

//* ndivo, ncono, sdivo, scono: objects; ndivs, ncons, sdivs, scons: strings
//** for use with INTF
iterator divr () { local ii
  for ii=0,ncl.count-1 {
    XO=ncl.object(ii)
    if (object_id(sfunc.obj(XO.pre.p))==object_id(cells[$1].object($2))) {
      iterator_statement
    }
  }
}

iterator conr () { local ii
  for ii=0,ncl.count-1 {
    XO=ncl.object(ii)
    if (object_id(sfunc.obj(XO.syn.p))==object_id(cells[$1].object($2))) {
      iterator_statement
    }
  }
}

//** iterators
// eg for syt("ns[0]","ns[1]") print XO,YO,nco
iterator syt () { local num,err
  err=0
  canobj($o1,"XO") canobj($o2,"YO") // canobj -- canonical object
  if ((num=cvode.netconlist(XO,"",YO).count)!=1) { 
    printf("syt() ERROR num==%d (%s,%s)\n",num,XO,YO) err=1 }
  if (!err) { 
    nco=cvode.netconlist(XO,"",YO).object(0)
    iterator_statement
  }
  XO=nil YO=nil nco=nil
}

// nca makes list backwards -- syn, then postcell, then precell
iterator nca () { local ii
  tmplist.remove_all
  if (numarg()==0) { cvode.netconlist("", "", "",tmplist)}
  if (numarg()==1) { cvode.netconlist("", "",   $s1,tmplist)}
  if (numarg()==2) { cvode.netconlist("",  $s2, $s1,tmplist)}
  if (numarg()==3) { cvode.netconlist($s3, $s2, $s1,tmplist)}
  for ii=0,tmplist.count-1 {
    XO=tmplist.object(ii)
    iterator_statement
  }
}

func wtvec () {
  revec(vec)
  for nca($s1) vec.append(XO.weight)
  vlk(vec)
  return vec.size
}

func dlyvec () {
  revec(vec)
  for nca($s1) vec.append(XO.delay)
  vlk(vec)
  return vec.size
}

iterator prtr () {  local ii
  tmplist.remove_all
  for ii=0,cvode.netconlist("", $s1, "").count-1 { 
    nco=cvode.netconlist("", $s1, "").object(ii)
    iterator_statement 
  }
}
iterator potr () {  local ii
  for ii=0,cvode.netconlist($s1,"","").count-1 {
    nco=cvode.netconlist($s1,"","").object(ii)
    iterator_statement 
  }
}
iterator sytr () {  local ii
  for ii=0,cvode.netconlist("","",$s1).count-1 { 
    nco=cvode.netconlist("","",$s1).object(ii)
    iterator_statement 
  }
}

proc sdivo () {for ltr(XO,cvode.netconlist($o1, "", "")) prsxo() }
func ndivo () { return cvode.netconlist($o1, "", "").count }
func ncono () { return cvode.netconlist("", $o1, "").count }
proc scono () {for ltr(XO,cvode.netconlist("", $o1, "")) prsxo() }
func ndivs () { 
  if (numarg()==1) return cvode.netconlist($s1, "", "").count else {
                   return cvode.netconlist("" , "", "").count }
}
func ncons () { return cvode.netconlist("", $s1, "").count }
proc sdivs () { for ltr(XO,cvode.netconlist($s1, "", "")) prsxo() }
proc scons () { 
  if (numarg()==0) for ltr(XO,cvode.netconlist("", "", ""))  prsxo()
  if (numarg()==1) for ltr(XO,cvode.netconlist("", $s1, ""))  prsxo()
  if (numarg()==2) for ltr(XO,cvode.netconlist("", $s1, $s2)) prsxo()
  if (numarg()==3) for ltr(XO,cvode.netconlist($s1, $s2,$s3)) prsxo()
}
 
// print pre,post,syntype,threshold,weight,delay
proc prsxo () {
  if (isobj(XO.precell,"NULLobject")) {
    printf("%s:%s->%s (%s) (t%g,w%g,d%g)\n",XO,XO.pre,XO.postcell,XO.syn,XO.threshold,XO.weight,XO.delay)
  } else {
    printf("%s:%s->%s (%s) (t%g,w%g,d%g)\n",XO,XO.precell,XO.postcell,XO.syn,XO.threshold,XO.weight,XO.delay)
  }
}

//*  avwtsp(), actwsc(), avwt() to get proper weight values
obfunc sysel () { local pr,po,sy,co,vb,num
  if (numarg()==1) sy=$1 else { // already selected
    pr=$1 po=$2 sy=$3 co=$4
    num=sp.select("CODE",EQU1,pr,"CODE",EQU2,po,"CODE",EQU3,co)
    if (num==0) { vec0.resize(0) 
      return vec0} // kludge to get an error returned
  }
  if ((pr==IN && sy==GA) || sy==AM) return sp.getcol("WT0",0) else return sp.getcol("WT1",0)
}

//** avwtsp(PR,PO,COL) looks at NQS db to determine weights -- 7x faster than avwt() below
//   avwtsp(sp) if do select before checking this
obfunc avwtsp () {  local a,s1,s2,pr,po,co,wd localobj sp0,v1,o,vr
  pr=$1 po=$2 co=$3 sp0=sp o=new Union()
  a=allocvecs(v1)
  if (name_declared("netstem")==4) {
    sprint(tstr,"%s%sS%02d.spnqs",netstem,CTYP.o(po).s,scale)
    if (! strm(tstr,sp0.file)) { printf("Reading %s (was %s)\n",tstr,sp0.file)
      sp0.rd(tstr)
    }
  }
  if (pr==IN) {s1=GA s2=GB} else {s1=AM s2=NM}
  if (numarg()==4) sp0=$o4.cob else {
    sp0.verbose=0
    vr=sysel(pr,po,s1,co)
    if (vr.size==0 && verbose) {printf("%s->%s EMPTY\n",CTYP.o(pr).s,CTYP.o(po).s) return sp0}}
  sp0.verbose=0
  wd=ncq.ay(pr,po,s1,co).x
  vr.mul(wd)
  stat(vr,v1) // 0size,1max,2min,3mean,4sdev
  o.x=v1.x[3]
  printf("%s->%s %s %s: max: %g; min: %g; mean: %g; sdev: %g\n",\
         CTYP.o(pr).s,CTYP.o(po).s,STYP.o(s1).s,INCOL.o(co).s,\
         v1.x[1],v1.x[2],v1.x[3],v1.x[4])
  wd=ncq.ay(pr,po,s2,co).x
  vr=sysel(s2)  vr.mul(wd)  stat(vr,v1)
  o.x[1]=v1.x[3]
  printf("%s->%s %s %s: max: %g; min: %g; mean: %g; sdev: %g\n",\
         CTYP.o(pr).s,CTYP.o(po).s,STYP.o(s2).s,INCOL.o(co).s,\
         v1.x[1],v1.x[2],v1.x[3],v1.x[4])
  sp0.verbose=1
  dealloc(a)
  return o
}

//** actwsc(PR,PO,SY,COCODE,COL,[,tMIN,tMAX])
func actwsc () { local ret,f1,a,wd,we,pr,po,min,max,u,x,y,sy,co,col,prspks,th,conv,diff localobj v1,v2,vr,q,p
  if (numarg()==0) {
    printf(" actwsc(PR,PO,SY,COCODE,COL,[,tMIN,tMAX])\n") return 0 }
  pr=$1 po=$2 sy=$3 co=$4 col=$5
  // min,max: just look at activity during this period
  if (numarg()>5) {min=$6 max=$7} else {min=0 max=tstop}
  if (min==max) f1=0 else f1=1
  a=allocvecs(v1,v2) q=new NQS() q.verbose=0 p=printlist
  if (f1) {
    // find name for raster plot should be eg 'SU'
    for ltr(XO,printlist) if (strc(XO.name,CTYP.o(pr).s)) {x=i1 break} 
    for ixt(pr,col) v2.append(ixi) // get the values for this col
    q.resize("time",p.o(x).tvec,"ind",p.o(x).vec)
    u=q.select("time","[]",min,max,"ind","[]",v2.min,v2.max) // u=number of spikes in period
    min=q.applf(".min","time") max=q.applf(".max","time") 
    q.out.v.sort()  y=q.out.v.uniq() // y= uniq spikers
    if (max-min<10) diff=10 else diff=max-min
    // prspks = y/diff/numc[pr]*10 // number of spikes in 1 cell spike per 10ms
    prspks = y/diff/numc[pr]*ncols*10 // normalize by no. of cols
    x=prspks*(conv=cp.fetch("PRID",pr,"POID",po,"CONV")) // incoming spks over 10ms
  }
  sp.verbose=0  vr=sysel(pr,po,sy,co) sp.verbose=1
  if (vr.size==0) {dealloc(a) return -1}
  vrsz(ce,v1,"O")
  for ixt(po) v1.append(XO.VTH-XO.RMP)
  stat(v1,v2) th=v2.x[3] // threhold stats
  wd=ncq.ay(pr,po,sy,co).x // current scaling factor
  we=4 // 4 cells must fire at about same time to get thresh-level drive
  if (f1) {
    printf("%s:%d,%d spks/%.2fms ~ %.2f spks/10ms (%d/%d)\n",\
           CTYP.o(pr).s,u,y,max-min,prspks,x,conv)
    printf("Thresh: %.2f+/-%.2f PSP: %.2f+/-%.2f\n",th,v2.x[4],vr.mean,vr.stdev)
    printf("\tFor 1/%d continuous activation: ",we)
    // wd *= threshold/spks_in/<PSP>/weighting // wd is scaling_factor
    printf("\nwmul(%s,%s,%s,%s,%g)\n",\
           CTYP.o(pr).t,CTYP.o(po).t,STYP.o(sy).t,INCOL.o(co).t,\
           th/x/vr.mean/we*wd)
    ret=th/x/vr.mean/we*wd
  } else {
    printf("%s->%s %s,%s: %g\n",\
           CTYP.o(pr).s,CTYP.o(po).s,STYP.o(sy).s,INCOL.o(co).s,\
           vr.mean*wd/th*100) // percent of threshold
    ret=vr.mean*wd/th*100
  }
  dealloc(a)
  return ret
}

//** avwt() goes through ncl list to sum weights
proc avwt () {  local pr,po,sy,co localobj nc,o
  pr=$1 po=$2 sy=$3 co=$4
  o=ncq.ay(pr,po,sy,co)
  vrsz(o.o.count,ind,"O")
  for ltr(nc,o.o) ind.append(nc.weight[sy])
  ind.mul(o.x) // mutiply by scale factor
  if (ind.mean!=0) {
    printf("%s%s %s: ",CTYP.o(pr).s,CTYP.o(po).s,STYP.o(sy).s) 
    stat(ind)
  }
}

//* grvec addendum -- new_printlist_
//** new_printlist_nc(obj,id[,NAME]) records spikes from obj to vitem in printlist
// with NAME creates a new vitem
warn_flag=1
obfunc new_printlist_nc () { local id,thresh,yoset localobj xo,yo,ox,cvn
  ox=$o1 id=$2 thresh=-20 yoset=0
  if (!isojt(ncl,tmplist)) ncl=new List()
  if (numarg()>=3) {
    if (argtype(3)==2) { 
      printlist.append(yo=new vitem($s3,100))
      yoset=1
    } else if (argtype(3)==0) { thresh=$3
    } else printf("new_printlist_nc arg 3 should be name or thresh\n")
  }
  if (numarg()>=4) thresh=$4
  if (!yoset) if (printlist.count==0) { 
    printlist.append(yo=new vitem("SPKS",100))
  } else { 
    yo=printlist.object(printlist.count-1) // continue with the last one
  }
  if (isobj(ox,"NetCon")) {  xo=ox // use this NetCon
  } else {
    cvn=cvode.netconlist(ox, "", "")
    if (cvn.count==0) { // no netcons available
      if (! warn_flag) printf(".") else { 
        warn_flag=0 // only warn first time
        printf("WARNING (new_printlist_nc) creating NetCon for %s ",ox) }
        if (ox.fflag) {  xo=new NetCon(ox, nil)
        } else { ox.soma xo=new NetCon(&v(x), nil) }
        xo.threshold=thresh
        ncl.append(xo)
    } else {
      xo=cvn.object(0) // the first of the postsyn netcons
      if (0) if (xo.threshold!=thresh && isassigned(xo.precell)) {
        printf("Warning: new_printlist_nc resetting thresh for %s\n",xo)
        xo.threshold=thresh
      }
    }
  }
  xo.record(yo.tvec,yo.vec,id)
  return yo
}

//** proc new_printlist_nc2(VITEM,CELL or NETCON,ID#) -- don't bother looking for an existing synapse
obfunc new_printlist_nc2 () { local id,thresh,yoset localobj xo,ox
  ox=$o1 xo=$o2 id=$3 thresh=-20 yoset=0
  if (isobj(xo,"NetCon")) { // use this NetCon
  } else {
    if ($o2.fflag) {  xo=new NetCon($o2, nil)
    } else { $o2.soma xo=new NetCon(&v(x), nil) }
    xo.threshold=thresh
    ncl.append(xo)
  }
  xo.record(ox.tvec,ox.vec,id)
  return xo
}

//** fipre(NETCON) find a presynaptic index for a netcon
func fipre () { local nu
  if (isassigned($o1.pre)) { // a point-process
    return objnum($o1.pre)
  } else {
    if ($o1.preloc==-1) {
      printf("fipre() ERR: %s\n",$o1) err() }
    sscanf(secname(),"%*[^[][%d]",&x) // find number of the cell
    // eg grvecstr="PYR2 sINT sTC  sRE" // use locations in this string as id
    if (sfunc.len(grvecstr)==0) {
      return x
    } else {
      sscanf(secname(),"%[^[]",tstr)
      nu=sfunc.substr(grvecstr,tstr)/5
      return nu*50+x
    }
    pop_section()
  }
}

//** fspks (index,vec) put spike times for index in vec
// fspks (index,vec,num)  use printlist.object(num).vec
// fspks (index,vec,spkvec,tvec)
proc fspks () { local a,b,ix
  a=b=allocvecs(2) b+=1
  if (numarg()==2) { ix=$1 XO=printlist.object(0).vec YO=printlist.object(0).tvec }
  if (numarg()==3) { ix=$1 XO=printlist.object($3).vec YO=printlist.object($3).tvec }
  if (numarg()==4) { ix=$1 XO=$o3 YO=$o4 }
  revec($o2)
  mso[a].indvwhere(XO,"==",ix)
  $o2.index(YO,mso[a])
  dealloc(a)
}

// fspks1(index,vec)
proc fspks1 () { local a,b,ix
  if (numarg()==2) { ix=$1 XO=printlist.object($1).vec YO=printlist.object($1).tvec }
  if (numarg()==4) { ix=$1 XO=$o3 YO=$o4 }
  $o2.resize(XO.size)
  $o2.xing(XO,YO,thresh) // times
}

//** spktimes(plist#,cell#)  prints out spike times from a tvec/ivec printlist item
// spktimes(tvec,vec,cell#)
func spktimes () { local sz,wh,a,i localobj va,vb,tv,vv,o
  a=allocvecs(va,vb)
  i=0 // i serves as flag for saving to a vector
  if (numarg()==4) {
    tv=$o1 vv=$o2 wh=$3 i=4
  } else if (numarg()==3) {
    if (argtype(1)==1) {
      tv=$o1 vv=$o2 wh=$3
    } else {
      o=printlist.object($1) wh=$2 vv=o.vec tv=o.tvec
      i=3
    }
  } else {
    if (numarg()==2) {
      o=printlist.object($1) wh=$2
    } else if (numarg()==1) {
      o=printlist.object(0)  wh=$1
    } else wh=-1
    vv=o.vec tv=o.tvec
  }
  sz=vv.size
  if (wh==-1) { vlk(tv,vv) 
  } else {
    va.indvwhere(vv,"==",wh)
    vb.index(tv,va)
    sz=vb.size
    if (i) $oi.copy(vb) else vlk(vb)
  }
  dealloc(a)
  return sz
}

//** ceqi() cvode.event_queue_info(3,...)
proc ceqi () { local flag,ii
  flag=$1
  tmplist.remove_all revec(ind,vec)
  printf("\n\n\t**** t=%g ****\n\n",t)
  if (flag==3) {
    cvode.event_queue_info(3, ind, vec, tmplist)
    for ltr(XO,tmplist,100) printf("t=%g flag=%d nc=%s\n",ind.x[i1],vec.x[i1],XO)
  } else if (flag==2) {
    cvode.event_queue_info(2, ind, tmplist)
    for ltr(XO,tmplist,100) { printf("t=%g %s->%s:",ind.x[i1],XO.pre,XO.syn)
      for ii=0,XO.wcnt-1 printf("%g ",XO.weight[ii])
      print ""
    }
  } else if (flag==0) {
    for ltr(XO,ncl,10) { printf("%s->%s:",XO.pre,XO.syn)
      for ii=0,XO.wcnt-1 printf("%g ",XO.weight[ii])
      print ""
    }
  }
}

//** new_printlist_ac(obj,name[,NAME,num]) 
// adds params from an artificial cell with built-in vec dumping
obfunc new_printlist_ac () { local i,max,num,svloc,sz,tsz localobj lo,st
  npacsz=1.2 st=new String()
  if(name_declared("vdt_INTF")==5) {
    vdt = vdt_INTF
    st.s="INTF"
  } else if(name_declared("vdt_INTF6")==5) {
    vdt = vdt_INTF6
    st.s="INTF6"
  } else print "new_printlist_ac warning: no INTF/INTF6 found!"
  lo = new Union()
  if (numarg()>=4) num=$4
  if (numarg()==4) { sprint(lo.s,"%s%d.%s",$s3,$4,$s2) 
  } else {           sprint(lo.s,"%s.%s",$o1,$s2)  }
  if (tstop<=1e4) tsz=tstop else tsz=1e4
  sz=npacsz*tsz/vdt
  vite = new vitem(lo.s,sz)
  vite.o=$o1
  if ($o1.recflag==0) {  // record tvec
    vite.tvec=new Vector(sz)
    $o1.initrec("time",vite.tvec)
  } else {
    for ltr(YO,printlist) if (eqobj($o1,YO.o)) vite.tvec=YO.tvec
    if (!isassigned(vite.tvec)) printf("ERR %s %s not found in printlist\n",st.s,$o1)
  }
  vite.tvflag=1
  $o1.initrec($s2,vite.vec)
  printlist.append(vite)
  return vite
}

objref pfih
psgchk = 1
proc psend () { local tt
  for ltr(XO,printlist) {
    if (strm(XO.var,"[()][()]$")) {
      revec(XO.vec,XO.tvec)
      for (tt=psgchk;tt<=tstop;tt+=psgchk) {
        sprint(tstr,"%s.append(%s)",XO.vec,XO.var)
        cvode.event(tt,tstr)
        sprint(tstr,"%s.append(t)",XO.tvec)
        cvode.event(tt,tstr)
      }
    }
  }
}

// eg new_printlist_fc("intf.M1()")
// obsolete since intf.mod no longer defines M1(),M2(),...
obfunc new_printlist_fc () {
  if (! isassigned(pfih)) { pfih = new FInitializeHandler("psend()") }
  if (! strm($s1,"\\(\\)$")){print "Should be func eg intf.M1()" return }
  XO = new vitem($s1,0)
  printlist.append(XO)
  return XO
}
  

//* misc and commentary
// con=absolute convergence number, div=absolute div number
// con = %con * pre
// div * pre = con * post = S (total synapses)
// %div = div/post = S/(pre*post) = con/pre = %con
// div = %con * post
// maxdiv = int((1 + var) * div)
 
//** vari returns randomly chosen $1+/-$2, $2 is a percent
func frani () { return int(rdm.uniform($1,$2+1)) }
func uvari () { return $1 - $1*$2 + (rdm.uniform(0,1) * 2 * $1*$2) }
func gvari () { return $1 - (rdm.normal(0,1) * $1*$2) }

//** clearsyns() 
proc clearsyns () { 
  for ltr(XO,ncl) XO.active(0)
  ncl.remove_all 
  for ltr(XO,new List("NetCon")) printf("**** CLEARSYN WARNING %s STILL EXISTS ****\n",XO)
}

// findstr() look thru list for a string and return the index
func findstr () { local flag
  flag = -1 // failure
  for ltr(XO,$o2) {
    if (strcmp($s1,XO.s) == 0) flag=i1
  }
  return flag
}

//** syntyp() reads syntyps out of cp directory
proc syntyp () { local x 
  x=$1
  if (x==EX) {$o2.x=AM $o2.x[1]=NM  // excit
  } else if (x==IX) {$o2.x=GA $o2.x[1]=GB // inhib
  } else $o2.x[1]=-1
}

//** excitp(PRE) -> T if cell is excitatory
func excitp () { 
  if ($1==SU || $1==DP || $1==TC) return 1
  if ($1==RE || $1==IN) return 0
  printf("ID not classified in excitp\n")
  return -1
}

//* Saving and reading
//** svnet() rdnet()
proc svnet () { local ii
  ii=0
  if (numarg()==1) filename=$s1 else {
    sprint(filename,"data/%s%c.net",datestr,97+ii)
    while (tmpfile.ropen(filename)) sprint(filename,"data/%s%c.net",datestr,97+ii+=1)
    printf("nqsnet output to %s\n",filename)
  }
  XO=usefiles("nqs.hoc","syncode.hoc","nqsnet.hoc","labels.hoc")
  cp.comment=XO.s
  sp.comment=XO.s
  cp.sv(filename)
  sp.sv(filename,1)
  tmpfile.close()
}

proc rdnet () {
  filename=$s1
  if (!isobj(cp,"NQS")) {cp=new NQS() sp=new NQS()}
  printf("reading network file %s\n",filename)
  cp.rd(filename)
  sp.rd() // continue reading
  tmpfile.close()
}

//** svstate() rdstate() rerun()
proc svstate () {
  if (! isobj(svs,"SaveState")) svs=new SaveState()
  svs.save
  tmpfile.wopen($s1)
  svs.fwrite(tmpfile)
}

proc rdstate () {
  if (! isobj(svs,"SaveState")) svs=new SaveState()
  tmpfile.ropen($s1)
  svs.fread(tmpfile)
  svs.restore
}

proc initrr () { printf("WARNING: initrr() -- init for rerun() -- not defined\n") }

proc rerun () { local tstp,told
  tstp=tstop
  if (numarg()>0) rdstate($s1) else svs.restore
  if (numarg()==1) if (argtype(1)==0) tstp=$1
  if (numarg()>1) tstp=$2
  told=t
  initrr()
  if (t!=told) {
    printf("ERROR: initrr() appears to have reset t (finitialize()?) (%d->%d)\n",told,t) 
    return
  }
  sprint(tstr,"cvode.solve(%g)",t+tstp)
  time(tstr)
  finish()
}

// mkvspks(NUMCELLS,NUMGROUPS,INVL1,SPREAD,[APPEND])
// produces ivspks index vector and vspks time vector for inserting activations
// groups are not stable -- they are reconstituted at each spiking time
// NUMCELLS   indices will be drawn from 0..NUMCELLS-1 
// NUMGROUPS  will give the separate groupings for simultaneous activation across a group
//            there will be NUMCELLS/NUMGROUPS cells in each group
// INVL1      is the average interval in ms between the inputs
// SPREAD     is the +/- SPREAD/2 sloppiness in this interval 
// APPEND     is a flag saying to add on to existing ivspks/vspks vectors
proc mkvspks () { local gcnt,numc,numg,inv1,inv2,a,ii,j,apflg localobj iv,tv
  numc=$1 numg=$2 inv1=$3 inv2=$4 apflg=0
  if (numarg()>=5) if ($5) apflg=1 // 5th arg is append flag
  if (!apflg) revec(ivspks,vspks)
  gcnt = int(numc/numg) // number of cells in each group
  a=allocvecs(iv,tv)
  for ({ii=20 j=0};ii<tstop-20;{ii+=inv1 j+=1e-8}) {
    iv.indgen(0,numc-1,1) // assume that start from intf#0
    shuffle(iv)
    iv.resize(gcnt) // now have a set of cells that need to spike at
    iv.add(j) // add a little bit to ensure that sortindex below preserves forward time
    tv.resize(gcnt)
    rdm.uniform(-inv2/2,inv2/2)
    tv.setrand(rdm)
    tv.add(ii)
    ivspks.append(iv) vspks.append(tv)
    if (j>=1-1e-8) printf("mkvspks ERR: j>=1\n") // should never happen
  }
  ivspks.sortindex(iv)
  tv.index(vspks,iv)  vspks.copy(tv)
  tv.index(ivspks,iv) ivspks.copy(tv)
  ivspks.rnd()
  wvspks.resize(ivspks.size)
  dealloc(a)
}

obfunc netconlist () { 
  if (numarg()==1) {
    if (argtype(1)==1) return cvode.netconlist("","",$o1)
    if (argtype(1)==2) return cvode.netconlist("",$s2,"")
  }
  return cvode.netconlist("","","") 
}

//* mkncq() -- this may not be useful since will slow down smap() still more
proc smap () { } // stub
proc mkncq () { local pr,po,sy,co localobj o
  if (!isassigned(ncq)) {
    ncq=new NQS("PRID","POID","SYN","INCOL","SCALE","NCL")
    ncq.odec("NCL")
    ncq.useslist("PRID",CTYP) ncq.useslist("POID",CTYP) 
    ncq.useslist("SYN",STYP)  ncq.useslist("INCOL",INCOL) 
  }
  ncq.clear()
  sp.out.mo(1) // for smap
  // no syn specificity for ARTCELLs -- assume no NM w/out AM, no GB w/out GA
  for case(&pr,DP,SU,IN) for case(&po,DP,SU,IN) for case(&co,INC,RIT,LFT) {
    if (sp.select("CODE",EQU1,pr,"CODE",EQU2,po,"CODE",EQU3,co) == 0) {
      // printf("%s->%s (%s) EMPTY\n",CTYP.o(pr).s,CTYP.o(po).s,INCOL.o(co).s)
    } else {
      if (pr==IN) sy=GA else sy=AM
      o=ncl[pr][po][co]
      ncq.append(pr,po,sy  ,co,1,o) // AM,GA
      ncq.append(pr,po,sy+1,co,1,o) // NM,GB
      smap(pr,po,sy,co,o)
    }
  }
}

proc smap () { local prdx,w0x,w1x,pr,po,sy,co localobj o
  pr=$1 po=$2 sy=$3 co=$4 o=$o5
  for ii=0,PRv.size-1 { // need prior sp.out.mo(1)
    w0x=WT0v.x[ii] w1x=WT1v.x[ii] 
    nc=new NetCon(ce.o(PRv.x[ii]),ce.o(POv.x[ii]),-10,DELv.x[ii],0)
    o.append(nc)
    for kk=0,nc.wcnt-1 nc.weight[kk]=0
    if (pr==IN) { 
      nc.weight[GA]=w0x
      if (w1x>0) nc.weight[GB]=w1x
    } else {
      nc.weight[AM]=w0x
      if (w1x>0) nc.weight[NM]=w1x
    }
  }
}

proc smap () { printf("NOW using dedicated smap() routines in network.hoc") }

//* nrstim - stim cells with noise/random input
// nrstim(celltype,[,maxt,clear,seed,nqs])
// maxt specifies time after which input turned off
// clear specifies whether to clear nqs passed in (or vq if declared)
// global params are pvar,pwght,prate,ppos described below
{declare("pvar",2)}     //variance of noise input
{declare("pwght",0)}    //mean value of noise input
{declare("prate",1000)} //sampling rate for random input
{declare("ppos",1)}     //use abs value only -- should be 1 otherwise intf cant handle it!!!
func nrstim () { local a,tt,i,j,ti,maxt,ct localobj rdp,vr,vt,vi,nq
  rdp=new Random() nq=nil
  a=allocvecs(vr,vt,vi)
  if(numarg()>1) maxt=$2 else maxt=tstop
  if(numarg()>4) nq=$o5 else if(name_declared("vq")) nq=vq
  if(nq==nil) {
    printf("nrstim ERRA: input nq==nil!\n") dealloc(a) return 0
  }
  if(numarg()>2) {
    if($3) nq.clear() 
  } else nq.clear()  
  if(numarg()>3) rdp.ACG($4) else rdp.ACG(1234)
  rdp.normal(pwght,pvar)
  vt.indgen(0,maxt,1000/prate) vi.resize(vt.size) vr.resize(vt.size)
  for i=ix[$1],ixe[$1] {
    vr.setrand(rdp) if(ppos) vr.abs()
    vi.fill(i)   nq.v[0].append(vi)    nq.v[1].append(vt)    nq.v[2].append(vr)
  }
  nq.sort("time") nq.sort("ind")
  intf.initvspks(nq.v,nq.v[1],nq.v[2])
  dealloc(a)
  return 1
}

//* subthstim(celltype,% of type,start,end,multiplier,bias,[nqs -- default is vq]) 
//this function takes the stim nqs , selects random % of cells (0-100), and between
//start and end times multiplies the stim values by multiplier and adds bias
//if the 7th argument is absent, the vq NQS should be declared, non-nil
//& should have been initialized with nrstim or setnrstim
func subthstim () { local startt,endt,fctr,ct,prc,bias,idx,a localobj vid,nq
  ct=$1 prc=$2 startt=$3 endt=$4 fctr=$5 bias=$6 a=allocvecs(vid) nq=nil
  if(numarg()>6) nq=$o7 else if(name_declared("vq")) nq=vq
  if(nq==nil) {
    printf("subthstim ERRA: input nq is nil!!\n")    dealloc(a)
    return 0
  }
  vid.resize(numc[ct]) vid.fill(0)
  vid.fill(1,0,int(numc[ct]*prc/100)) vid.shuffle()
  idx=nq.select(-1,"ind","[]",ix[ct],ixe[ct],"time","[]",startt,endt) printf("ni=%d\n",idx)
  for vtr(&idx,nq.ind) if(nq.v[2].x(idx)>0 && vid.x(nq.v[0].x(idx)-ix[ct])>0) {
    nq.v[2].x(idx) = nq.v[2].x(idx)*fctr+bias
  }
  dealloc(a)
  return 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]

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)