Computational Surgery (Lytton et al. 2011)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:140881
Figure 2 in Neocortical simulation for epilepsy surgery guidance: Localization and intervention, by William W. Lytton, Samuel A. Neymotin, Jason C. Wester, and Diego Contreras in Computational Surgery and Dual Training, Springer, 2011
Reference:
1 . Lytton WW, Neymotin SA, Wester JC, Contreras D (2011) Neocortical simulation for epilepsy surgery guidance: Localization and intervention Computational Surgery and Dual Training
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 fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Epilepsy;
Implementer(s): Lytton, William [billl at neurosim.downstate.edu];
/
b11aug17
data
readme.txt
intf6_.mod *
misc.mod *
nstim.mod *
stats.mod *
vecst.mod *
col.hoc
declist.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
gload.hoc
grvec.hoc
init.hoc
labels.hoc
local.hoc *
misc.h *
mosinit.hoc
network.hoc
nqs.hoc
nqs_utils.hoc *
nqsnet.hoc *
nrnoc.hoc *
params.hoc
run.hoc
setup.hoc *
simctrl.hoc *
spkts.hoc *
stats.hoc
stdgui.hoc *
syncode.hoc
                            
// $Id: spkts.hoc,v 1.86 2010/07/10 02:32:11 samn Exp $

// load_file("spkts.hoc")
// ancillary programs for handling vectors

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


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

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

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

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

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

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

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

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

burst_time=0
burst_maxfreq = 30
calc_ave = 0

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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



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

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

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


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

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

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

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

Lytton WW, Neymotin SA, Wester JC, Contreras D (2011) Neocortical simulation for epilepsy surgery guidance: Localization and intervention Computational Surgery and Dual Training

References and models cited by this paper

References and models that cite this paper

Babloyantz A, Destexhe A (1986) Low-dimensional chaos in an instance of epilepsy. Proc Natl Acad Sci U S A 83:3513-7 [PubMed]

Benifla M, Otsubo H, Ochi A, Snead OC, Rutka JT (2006) Multiple subpial transections in pediatric epilepsy: indications and outcomes. Childs Nerv Syst 22:992-8 [PubMed]

Brown SP, Hestrin S (2009) Intracortical circuits of pyramidal neurons reflect their long-range axonal targets. Nature 457:1133-6 [PubMed]

Buonomano DV (2009) Harnessing chaos in recurrent neural networks. Neuron 63:423-5 [PubMed]

Carnevale NT, Hines ML (2006) The NEURON Book

Clusmann H, Kral T, Gleissner U, Sassen R, Urbach H, Blumcke I, Bogucki J, Schramm J (2004) Analysis of different types of resection for pediatric patients with temporal lobe epilepsy. Neurosurgery 54:847-59; discussion 859-60 [PubMed]

Cohen-Gadol AA, Stoffman MR, Spencer DD (2003) Emerging surgical and radiotherapeutic techniques for treating epilepsy. Curr Opin Neurol 16:213-9 [PubMed]

Crampin EJ, Halstead M, Hunter P, Nielsen P, Noble D, Smith N, Tawhai M (2004) Computational physiology and the Physiome Project. Exp Physiol 89:1-26 [PubMed]

Duarte NC, Becker SA, Jamshidi N, Thiele I, Mo ML, Vo TD, Srivas R, Palsson BØ (2007) Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci U S A 104:1777-82 [PubMed]

Dubitzky W (2006) Understanding the computational methodologies of systems biology. Brief Bioinform 7:315-7 [PubMed]

Groh A, Meyer HS, Schmidt EF, Heintz N, Sakmann B, Krieger P (2010) Cell-type specific properties of pyramidal neurons in neocortex underlying a layout that is modifiable depending on the cortical area. Cereb Cortex 20:826-36 [PubMed]

   [48 reconstructed morphologies on NeuroMorpho.Org]

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]

Jahr CE, Stevens CF (1990) A quantitative description of NMDA receptor-channel kinetic behavior. J Neurosci 10:1830-7 [PubMed]

Jahr CE, Stevens CF (1990) Voltage dependence of NMDA-activated macroscopic conductances predicted by single-channel kinetics. J Neurosci 10:3178-82 [PubMed]

Kossoff EH, Ritzl EK, Politsky JM, Murro AM, Smith JR, Duckrow RB, Spencer DD, Bergey GK (2004) Effect of an external responsive neurostimulator on seizures and electrographic discharges during subdural electrode monitoring. Epilepsia 45:1560-7 [PubMed]

Lesser RP, Crone NE, Webber WR (2010) Subdural electrodes. Clin Neurophysiol 121:1376-92

Lesser RP, Crone NE, Webber WR (2011) Using subdural electrodes to assess the safety of resections. Epilepsy Behav 20:223-9 [PubMed]

Lesser RP, Kim SH, Beyderman L, Miglioretti DL, Webber WR, Bare M, Cysyk B, Krauss G, Gordon (1999) Brief bursts of pulse stimulation terminate afterdischarges caused by cortical stimulation. Neurology 53:2073-81 [PubMed]

Lytton W, Hellman K, Sutula T (1996) Computer network model of mossy fiber sprouting in dentate gyrus Epilepsia-aes Proceedings 37 S 5:117

Lytton WW, Hellman KM, Sutula TP (1998) Computer models of hippocampal circuit changes of the kindling model of epilepsy. Artif Intell Med 13:81-97 [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, 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

Milton JG (2010) Epilepsy as a dynamic disease: a tutorial of the past with an eye to the future. Epilepsy Behav 18:33-44 [PubMed]

Morrell F, Whisler WW, Bleck TP (1989) Multiple subpial transection: a new approach to the surgical treatment of focal epilepsy. J Neurosurg 70:231-9 [PubMed]

Motamedi GK, Salazar P, Smith EL, Lesser RP, Webber WR, Ortinski PI, Vicini S, Rogawski MA (2006) Termination of epileptiform activity by cooling in rat hippocampal slice epilepsy models. Epilepsy Res 70:200-10 [PubMed]

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]

Osorio I, Frei MG, Manly BF, Sunderam S, Bhavaraju NC, Wilkinson SB (2001) An introduction to contingent (closed-loop) brain electrical stimulation for seizure blockage, to ultra-short-term clinical trials, and to multidimensional statistical analysis of therapeutic efficacy. J Clin Neurophysiol 18:533-44 [PubMed]

Polsky A, Mel BW, Schiller J (2004) Computational subunits in thin dendrites of pyramidal cells. Nat Neurosci 7:621-7 [Journal] [PubMed]

   CA1 pyramidal neuron: as a 2-layer NN and subthreshold synaptic summation (Poirazi et al 2003) [Model]

Richardson KA, Schiff SJ, Gluckman BJ (2005) Control of traveling waves in the Mammalian cortex. Phys Rev Lett 94:028103-92

Rutecki P (1990) Anatomical, physiological, and theoretical basis for the antiepileptic effect of vagus nerve stimulation. Epilepsia 31 Suppl 2:S1-6

Schramm J, Aliashkevich AF, Grunwald T (2002) Multiple subpial transections: outcome and complications in 20 patients who did not undergo resection. J Neurosurg 97:39-47 [PubMed]

Schulz R, Hoppe M, Boesebeck F, Gyimesi C, Pannek HW, Woermann FG, May T, Ebner A (2011) Analysis of reoperation in mesial temporal lobe epilepsy with hippocampal sclerosis. Neurosurgery 68:89-97; discussion 97 [PubMed]

Spencer SS, Schramm J, Wyler A, O'Connor M, Orbach D, Krauss G, Sperling M, Devinsky O, Elger (2002) Multiple subpial transection for intractable partial epilepsy: an international meta-analysis. Epilepsia 43:141-5 [PubMed]

Steriade M (2004) Neocortical cell classes are flexible entities. Nat Rev Neurosci 5:121-34 [PubMed]

Tonnesen J, Sorensen AT, Deisseroth K, Lundberg C, Kokaia M (2009) Optogenetic control of epileptiform activity. Proc Natl Acad Sci U S A 106:12162-7

Ty L, Yorke J (1975) Period three implies chaos Amer Math Monthly 82:985-992

Wyler AR (1997) Recent advances in epilepsy surgery: temporal lobectomy and multiple subpial transections. Neurosurgery 41:1294-301; discussion 1301-2 [PubMed]

Wyler AR, Hermann BP, Somes G (1995) Extent of medial temporal resection on outcome from anterior temporal lobectomy: a randomized prospective study. Neurosurgery 37:982-90; discussion 990-1 [PubMed]

Yogarajah M, Focke NK, Bonelli SB, Thompson P, Vollmar C, McEvoy AW, Alexander DC, Symms MR, (2010) The structural plasticity of white matter networks following anterior temporal lobe resection. Brain 133:2348-64 [PubMed]

(41 refs)