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: network.hoc,v 1.160 2011/08/17 15:17:08 billl Exp $

//* Numbers and connectivity params

declare("colW",7,"colH",2,"torus",0)
declare("numcols",colW*colH)
declare("dbgcols",0) // whether to debug columns by making them have the same wiring and inputs
declare("colr",2) // maximal trans-column projection distance; 0 within col; 1 next col etc
declare("colnq","o[5]","lcol",new List())
{sprint(tstr,"o[%d]",numcols) declare("col",tstr)}
{sprint(tstr,"o[%d][%d]",colH,colW) declare("gcol",tstr)} // 2D column grid
double div[CTYPi][CTYPi][colr+1]//div[i][j]==# of outputs from type i->j
double wmat[CTYPi][CTYPi][STYPi][colr+1] // wmat[i][j][k]==weight from type i->j for synapse k
double delm[CTYPi][CTYPi]//avg. delay from type i->j
double deld[CTYPi][CTYPi]//delay variance from type i->j
double conv[CTYPi][CTYPi][colr+1]
dosetpmat=name_declared("pmat")==0
{sprint(tstr,"d[%d][%d][%d]",CTYPi,CTYPi,colr+1) declare("pmat",tstr)}
double prumat[CTYPi][CTYPi] //pruning matrix:prumat[i][j] specifies ratio (0-1) of synapses to prune
double sprmat[CTYPi][CTYPi] //sprouting matrix:sprmat[i][j] specifies ratio (0-1) to sprout i->j pathway with
double synloc[CTYPi][CTYPi]//location of synapses

declare("EEGain",4*15/11,"EIGain",16,"IEGain",8,"IIGain",8)
// declare("NMAMR",0.1,"EENMGain",1,"EIGainInterC",0.125*4,"EEGainInterC",0.1625*4)
declare("NMAMR",0.1,"EENMGain",1,"EIGainInterC",0.5,"EEGainInterC",0.6)

//* prdiv() - print div
proc prdiv () { local ii,jj
  for ii=0,CTYPi-1 for jj=0,CTYPi-1 if(div[ii][jj][0]) {
    printf("div[%s][%s][0]=%g\n",CTYP.o(ii).s,CTYP.o(jj).s,div[ii][jj][0])
  }
}

// %con (con/pre) = %div (div/post)
DEAD_DIV_INTF6=0
declare("jcn",1)
declare("disinhib",0) //iff==1 , turn off inhibition, by setting wmat[I%d][...]==0 in inhiboff()
declare("scale",4)//16//8//4
declare("pmatscale",1/scale) // scale for pmat - allows keeping it fixed while changing # of cells in network

batch_flag=declare("dstr",datestr,"setdviPT",NORM)
declare("params","not batch","ofile",output_file)
declare("dvseed",534023) // seed for wiring

dosetcpercol=name_declared("cpercol")==0 // whether to set values in cpercol or use user-supplied values
{sprint(tstr,"d[%d]",CTYPi) declare("cpercol",tstr)} // cells of a specific type per column
declare("vcpercol",new Vector(CTYPi))
declare("E5BNumF",1,"E5RNumF",1) // factors for # of E5 cells
declare("newkmjnums",0) // use #s based on KMJ #s in  /u/samn/vcsim/data/Cell_Numbers.xlsx columns R, T, V
declare("E2E2Inter",1) // whether to connect E2->E2 intercol
declare("E5RFctr",1.5) // used to control level of recurrent of E5
declare("E2RFctr",1) // used to control level of recurrence of E2

//* setcpercol - set # of cells per column
proc setcpercol () { local i // (/u/samn/vcsim/notebook.dol_1:24562)(notebook.dol_1:24492)
  if(dosetcpercol) { // if user didn't supply values (default), set # of cells of a type per column
    if(newkmjnums) {

      // based on KMJ #s in  /u/samn/vcsim/data/Cell_Numbers.xlsx columns R, T, V

      cpercol[E2] = 169 * scale
      cpercol[I2] = 48 * scale
      cpercol[I2L] = 8 * scale

      cpercol[E4] = 83 * scale
      cpercol[I4] = 24 * scale
      cpercol[I4L] = 4 * scale

      cpercol[E5R] = 93 * scale
      cpercol[E5B] = 32 * scale
      cpercol[I5] = 36 * scale
      cpercol[I5L] = 6 * scale

      cpercol[E6] = 218 * scale
      cpercol[I6] = 62 * scale
      cpercol[I6L] = 11 * scale

    } else {
      cpercol[E2]  = 150 * scale
      cpercol[E4] =   30 * scale
      cpercol[E5B] =  int(17 * scale * E5BNumF)
      cpercol[E5R] =  int(65 * scale * E5RNumF)
      cpercol[E6] =   60 * scale
      cpercol[I2L] =  13 * scale
      cpercol[I2]  =  25 * scale
      cpercol[I4L] =  14 * scale 
      cpercol[I4]  =  20 * scale
      cpercol[I5L] =  13 * scale
      cpercol[I5]  =  25 * scale
      cpercol[I6L] =  13 * scale
      cpercol[I6] =   25 * scale
    }
  }
  {vcpercol.resize(CTYPi) vcpercol.fill(0)} // store the values in a vector
  for i=0,CTYPi-1 vcpercol.x(i)=cpercol[i]
}

//* setpmat()
proc setpmat () { local pre,po
  if(!dosetpmat) return // if pmat setup by user (in notebook), then don't reset its values
  for ii=0,CTYPi-1 for jj=0,CTYPi-1 for kk=0,1 pmat[ii][jj][kk]=0
  pmat[E2][E2][0]=0.187 * E2RFctr
  if(E2E2Inter) pmat[E2][E2][1]=0.14 * E2RFctr
  pmat[E2][E4][0]=0.024
  pmat[E2][E5B][0]=0.024
  pmat[E2][E5R][0]=0.057
  pmat[E2][E6][0]=0
  pmat[E2][I2L][0]=0.51
  pmat[E2][I2][0]=0.43
  pmat[E2][I2][1]=0.14
  pmat[E4][E2][0]=0.145
  pmat[E4][E4][0]=0.243
  pmat[E4][E5B][0]=0.122
  pmat[E4][E5R][0]=0.116
  pmat[E4][E6][0]=0.032
  pmat[E4][I4L][0]=0.51
  pmat[E4][I4][0]=0.43
  pmat[E4][I4][1]=0.14
  pmat[E5B][E2][0]=0.018
  pmat[E5B][E2][1]=0.25
  pmat[E5B][E2][2]=0.1
  pmat[E5B][E4][0]=0.007
  pmat[E5B][E5B][0]=0.07 * E5RFctr
  pmat[E5B][E5B][1]=0.25 * E5RFctr
  pmat[E5B][E5B][2]=0.1 * E5RFctr
  pmat[E5B][E5R][0]=0.017 * E5RFctr
  pmat[E5B][E5R][1]=0.25 * E5RFctr
  pmat[E5B][E5R][2]=0.1 * E5RFctr
  pmat[E5B][E6][0]=0.07
  pmat[E5B][I2L][1]=0.14
  pmat[E5B][I2L][2]=0.07
  pmat[E5B][I5L][0]=0.51
  pmat[E5B][I5L][1]=0.14
  pmat[E5B][I5L][2]=0.07
  pmat[E5B][I5][0]=0.43
  pmat[E5B][I5][1]=0.14
  pmat[E5B][I5][2]=0.07
  pmat[E5R][E2][0]=0.022
  pmat[E5R][E4][0]=0.007
  pmat[E5R][E5B][0]=0.08 * E5RFctr
  pmat[E5R][E5B][1]=0.25 * E5RFctr
  pmat[E5R][E5R][0]=0.191 * E5RFctr
  pmat[E5R][E5R][1]=0.14 * E5RFctr
  pmat[E5R][E6][0]=0.032
  pmat[E5R][I5L][0]=0.51
  pmat[E5R][I5][0]=0.43
  pmat[E5R][I5][1]=0.14
  pmat[E6][E2][0]=0
  pmat[E6][E4][0]=0
  pmat[E6][E5B][0]=0.028
  pmat[E6][E5R][0]=0.006
  pmat[E6][E6][0]=0.028
  pmat[E6][I6L][0]=0.51
  pmat[E6][I6][0]=0.43
  pmat[E6][I6][1]=0.14
  pmat[I2L][E2][0]=0.35
  pmat[I2L][E5B][0]=0.5
  pmat[I2L][E5R][0]=0.35
  pmat[I2L][E6][0]=0.25
  pmat[I2L][I2L][0]=0.09
  pmat[I2L][I2][0]=0.53
  pmat[I2L][I5][0]=0.53
  pmat[I2L][I6][0]=0.53
  pmat[I2][E2][0]=0.44
  pmat[I2][I2L][0]=0.34
  pmat[I2][I2][0]=0.62
  pmat[I4L][E4][0]=0.35
  pmat[I4L][I4L][0]=0.09
  pmat[I4L][I4][0]=0.53
  pmat[I4][E4][0]=0.44
  pmat[I4][I4L][0]=0.34
  pmat[I4][I4][0]=0.62
  pmat[I5L][E2][0]=0.35
  pmat[I5L][E5B][0]=0.35
  pmat[I5L][E5R][0]=0.35
  pmat[I5L][E6][0]=0.25
  pmat[I5L][I2][0]=0.53
  pmat[I5L][I5L][0]=0.09
  pmat[I5L][I5][0]=0.53
  pmat[I5L][I6][0]=0.53
  pmat[I5][E5B][0]=0.44
  pmat[I5][E5R][0]=0.44
  pmat[I5][I5L][0]=0.34
  pmat[I5][I5][0]=0.62
  pmat[I6L][E2][0]=0.35
  pmat[I6L][E5B][0]=0.25
  pmat[I6L][E5R][0]=0.25
  pmat[I6L][E6][0]=0.35
  pmat[I6L][I2][0]=0.53
  pmat[I6L][I5][0]=0.53
  pmat[I6L][I6L][0]=0.09
  pmat[I6L][I6][0]=0.53
  pmat[I6][E6][0]=0.44
  pmat[I6][I6L][0]=0.34
  pmat[I6][I6][0]=0.62
}

//* scalepmat(fctr) - multiply values in pmat by fctr
proc scalepmat () { local fctr,from,to,cl
  fctr=$1
  for from=0,CTYPi-1 for to=0,CTYPi-1 for cl=0,1 pmat[from][to][cl] *= fctr
}

//* pmat2nq - return an NQS with info in pmat
obfunc pmat2nq () { local i,j,k localobj nqpmat
  nqpmat=new NQS("froms","tos","from","to","cold","pij")
  {nqpmat.strdec("froms") nqpmat.strdec("tos")}
  for i=0,CTYPi-1 for j=0,CTYPi-1 for k=0,colr if(pmat[i][j][k]) {
    nqpmat.append(CTYP.o(i).s,CTYP.o(j).s,i,j,k,pmat[i][j][k])
  }  
  return nqpmat
}

//* nq2pmat - load NQS ($o1) into pmat
proc nq2pmat () { local i,j,k localobj nq,vf,vto,vc,vpij
  {nq=$o1 nq.tog("DB") vf=nq.getcol("from") vto=nq.getcol("to") vc=nq.getcol("cold") vpij=nq.getcol("pij")}
  for i=0,CTYPi-1 for j=0,CTYPi-1 for k=0,colr pmat[i][j][k]=0
  for i=0,vf.size-1 pmat[vf.x(i)][vto.x(i)][vc.x(i)]=vpij.x(i)
  print "loaded " , nq , " into pmat"
}

//* synapse locations DEND SOMA AXON
proc setsynloc () { local from,to
  for from=0,CTYPi-1 for to=0,CTYPi-1 {
    if(ice(from)) {
      if(IsLTS(from)) {
        synloc[from][to]=DEND // distal [GA2] - from LTS
      } else {
        synloc[from][to]=SOMA // proximal [GA] - from FS
      }
    } else {
      synloc[from][to]=DEND // E always distal. use AM2,NM2
    }
  }
}

//* setdelmats -- setup delm,deld
proc setdelmats () { local from,to,ii,jj
  for from=0,CTYPi-1 for to=0,CTYPi-1 {
    if(synloc[from][to]==DEND) {
      delm[from][to]=4
      deld[from][to]=1
    } else {
      delm[from][to]=2.0
      deld[from][to]=0.2
    }
  }
  // snum=0
  // for ii=0,CTYPi-1 for jj=0,CTYPi-1 snum+=int(pmat[ii][jj][0]*numc[ii]*numc[jj]+1)
}

//* weight params
//** delay all 2+/-0.02 within column for now
proc setwmat () { local from,to,sy,gn,c
  for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=0,STYPi-1 for c=0,colr wmat[from][to][sy][c]=0

  wmat[E2][E2][AM2][0]=0.78
  wmat[E2][E2][AM2][1]=0.47 * EEGainInterC / 5.0
  wmat[E2][E4][AM2][0]=0.36
  wmat[E2][E5B][AM2][0]=0.36
  wmat[E2][E5R][AM2][0]=0.93
  wmat[E2][E6][AM2][0]=0
  wmat[E2][I2L][AM2][0]=0.23

  wmat[E2][I2][AM2][0] = 0.23
  wmat[E2][I2][AM2][1] = 1.5 * EIGainInterC

  wmat[E4][E2][AM2][0]=0.58
  wmat[E4][E4][AM2][0]=0.95
  wmat[E4][E5B][AM2][0]=1.01
  wmat[E4][E5R][AM2][0]=0.54
  wmat[E4][E6][AM2][0]=2.27
  wmat[E4][I4L][AM2][0]=0.23

  wmat[E4][I4][AM2][0] = 0.23
  wmat[E4][I4][AM2][1] = 1.5 * EIGainInterC

  wmat[E5B][E2][AM2][0]=0.26
  wmat[E5B][E2][AM2][1]=0.47 * EEGainInterC
  wmat[E5B][E2][AM2][2]=0.47 * EEGainInterC
  wmat[E5B][E4][AM2][0]=0.17
  wmat[E5B][E5B][AM2][0]=0.71
  wmat[E5B][E5B][AM2][1]=0.47 * EEGainInterC * 1.5
  wmat[E5B][E5B][AM2][2]=0.47 * EEGainInterC
  wmat[E5B][E5R][AM2][0]=0.24
  wmat[E5B][E5R][AM2][1]=0.47 * EEGainInterC * 1.5
  wmat[E5B][E5R][AM2][2]=0.47 * EEGainInterC
  wmat[E5B][E6][AM2][0]=0.49

  wmat[E5B][I2L][AM2][1]=1.5 * EIGainInterC
  wmat[E5B][I2L][AM2][2]=1.5 * EIGainInterC

  wmat[E5B][I5L][AM2][0]=0.23
  wmat[E5B][I5L][AM2][1]=1.5 * EIGainInterC
  wmat[E5B][I5L][AM2][2]=1.5 * EIGainInterC

  wmat[E5B][I5][AM2][0]=0.23
  wmat[E5B][I5][AM2][1]=1.5 * EIGainInterC
  wmat[E5B][I5][AM2][2]=1.5 * EIGainInterC

  wmat[E5R][E2][AM2][0]=0.67
  wmat[E5R][E4][AM2][0]=0.48
  wmat[E5R][E5B][AM2][0]=0.88
  wmat[E5R][E5B][AM2][1]=0.47 * EEGainInterC * 1.5
  wmat[E5R][E5R][AM2][0]=0.66 * 1.5
  wmat[E5R][E5R][AM2][1]=0.47 * EEGainInterC * 1.5
  wmat[E5R][E6][AM2][0]=0.28
  wmat[E5R][I5L][AM2][0]=0.23
  wmat[E5R][I5][AM2][0]=0.23
  wmat[E5R][I5][AM2][1]=1.5 * EIGainInterC

  wmat[E6][E2][AM2][0]=0
  wmat[E6][E4][AM2][0]=0
  wmat[E6][E5B][AM2][0]=0.53
  wmat[E6][E5R][AM2][0]=0.08
  wmat[E6][E6][AM2][0]=0.53
  wmat[E6][I6L][AM2][0]=0.23
  wmat[E6][I6][AM2][0]=0.23
  wmat[E6][I6][AM2][1]=1.5 * EIGainInterC

  wmat[I2L][E2][GA2][0]=0.83
  wmat[I2L][E5B][GA2][0]=0.83
  wmat[I2L][E5R][GA2][0]=0.83
  wmat[I2L][E6][GA2][0]=0.83
  wmat[I2L][I2L][GA2][0]=1.5
  wmat[I2L][I2][GA2][0]=1.5
  wmat[I2L][I5][GA2][0]=0.83
  wmat[I2L][I6][GA2][0]=0.83

  wmat[I2][E2][GA][0]=1.5
  wmat[I2][I2L][GA][0]=1.5
  wmat[I2][I2][GA][0]=1.5

  wmat[I4L][E4][GA2][0]=0.83
  wmat[I4L][I4L][GA2][0]=1.5
  wmat[I4L][I4][GA2][0]=1.5

  wmat[I4][E4][GA][0]=1.5
  wmat[I4][I4L][GA][0]=1.5
  wmat[I4][I4][GA][0]=1.5

  wmat[I5L][E2][GA2][0]=0.83
  wmat[I5L][E5B][GA2][0]=0.83
  wmat[I5L][E5R][GA2][0]=0.83
  wmat[I5L][E6][GA2][0]=0.83
  wmat[I5L][I2][GA2][0]=0.83
  wmat[I5L][I5L][GA2][0]=1.5
  wmat[I5L][I5][GA2][0]=1.5
  wmat[I5L][I6][GA2][0]=0.83

  wmat[I5][E5B][GA][0]=1.5
  wmat[I5][E5R][GA][0]=1.5
  wmat[I5][I5L][GA][0]=1.5
  wmat[I5][I5][GA][0]=1.5

  wmat[I6L][E2][GA2][0]=0.83
  wmat[I6L][E5B][GA2][0]=0.83
  wmat[I6L][E5R][GA2][0]=0.83
  wmat[I6L][E6][GA2][0]=0.83
  wmat[I6L][I2][GA2][0]=0.83
  wmat[I6L][I5][GA2][0]=0.83
  wmat[I6L][I6L][GA2][0]=1.5
  wmat[I6L][I6][GA2][0]=1.5

  wmat[I6][E6][GA][0]=1.5
  wmat[I6][I6L][GA][0]=1.5
  wmat[I6][I6][GA][0]=1.5

  //set NMDA weights
  for from=0,CTYPi-1 for to=0,CTYPi-1 for c=0,colr wmat[from][to][NM2][c]=NMAMR*wmat[from][to][AM2][c]
  //gain control
  for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=AM,GA2 for c=0,colr if(wmat[from][to][sy][c] > 0) {
    if(ice(from)) {
      if(ice(to)) {
        gn = IIGain
      } else {
        gn = IEGain
      }
      if(IsLTS(from) && !IsLTS(to)) gn *= 0.5
    } else {
      if(ice(to)) {
        gn = EIGain 
        if(IsLTS(to)) gn *= 0.5
      } else {
        gn = EEGain
        if(sy==NM || sy==NM2) gn *= EENMGain // E->E NMDA gain
      }
    }
    wmat[from][to][sy][c] *= gn 
  }
}

// %con (con/pre) = %div (div/post)

//* prune using values in prumat
proc pruc () { local i,j
  for i=0,CTYPi-1 for j=0,CTYPi-1{
      if(div[i][j][0] && numc[i] && numc[j] && prumat[i][j]){
        printf("Warning: pruning random %.2f%% of %s->%s syns\n",prumat[i][j]*100,CTYP.o(i).s,CTYP.o(j).s)
        for ixt(i) XO.prune(prumat[i][j],j)
      }
  }
}

//* get sprouting value assuming 0% sprouting == 50% pruning
func getspr () { local pr
  pr = $1
  return ((0.5-pr)/.5)*100
}

//* turn off pruning
proc pruoff () { local i,j
 for i=0,CTYPi-1 for j=0,CTYPi-1 prumat[i][j]=0
 for i=0,allcells-1 INTF6[i].prune(0)
}

//* set all entries in pruning matrix to $1
proc setpru () { local from,to,val
  val=$1
  pruoff() // first turn off pruning
  for from=0,CTYPi-1 for to=0,CTYPi-1 prumat[from][to]=val
}

//* print prumat
proc prumatpr () { local i,j
  for i=0,CTYPi-1 { for j=0,CTYPi-1{
      printf("%.2f  ",prumat[i][j])
   }
   printf("\n")
  }
}

//* the Killers
//** ColKillType(column id, type id, flag==1 to kill, 0 to unkill)
proc ColKillType () { local i,cidx,f,ty
  cidx=$1 ty=$2 f=$3
  for i=col[cidx].ix[ty],col.ix[ty]+col.numc[ty]-1 col[cidx].ce.o(i).flag("dead",f)
}

//** unkill() will set all the dead to undead
proc unkill () { for ii=0,numcols-1 col[ii].ce.o(0).flag("dead",0,1) }

//** kill2(column,flag==1 to kill, 0 to unkill) -- only effects cells in L2 of specified column
proc kill2 () { local f,c,i
  c=$1 f=$2
  for case(&i,E2,I2,I2L) ColKillType(c,i,f)
}

//** kill5(column,flag==1 to kill, 0 to unkill) -- only effects cells in L5 of specified column
proc kill5 () { local f,c,i
  c=$1 f=$2
  for case(&i,E5R,E5B,I5,I5L) ColKillType(c,i,f)
}

//** getkillids - gets ids of cells to kill in $o1 but excludes cells that are stim'ed
//$1=cell type to kill,$2=prct of cells to kill,$o3=vq stim nqs,$4=out vector of kill ids,$5=rnd seed
func getkillids () { local killcnt,i,j,ct,prct localobj vq,vkid,rd
  ct=$1 prct=$2 vq=$o3 vkid=$o4 killcnt=int(prct*numc[ct]) vkid.resize(0) j=0 i=ix[ct]
  rd=new Random() rd.ACG($5)
  while(j<killcnt){
    i=rd.discunif(ix[ct],ixe[ct])
    if(!vq.v[0].contains(i)){
      j+=1
      vkid.append(i)
    }
    i+=1
  }
  return killcnt
}

//* read .net file
strdef netfile
{sp = new NQS() cp = new NQS()}

//* CREATE CELLS
// %con (con/pre) = %div (div/post)
n=ty=id=0

//* sprcells() sprout cells in specific pathways using sprmat, $1=seed for rand generator
//max div is still 0.75*poty
func sprcells () { local id,a,prty,poty,sz,ls,mx,i localobj vid,vnewid,vnewdel,rd,vd,vtmp
  ls=$1 a=allocvecs(vid,vnewid,vnewdel,vd,vtmp) rd=new Random() rd.ACG(ls)
  for prty=0,CTYPi-1 for poty=0,CTYPi-1 if(sprmat[prty][poty]) for id=ix[prty],ixe[prty] {
    ce.o(id).getdvi(vid) ce.o(id).getdvi(0.2,vd)
    sz=div[prty][poty][0]*sprmat[prty][poty] mx=0.75*numc[poty]
    if(vd.x(poty)>=mx)continue//already @ max size
    while(sz+vd.x(poty)>mx) sz-=1
    vrsz(sz*4,vtmp,vnewid) rd.discunif(ix[poty],ixe[poty]) vtmp.setrand(rd)
    vtmp.uniq(vnewid) vtmp.resize(0)
    for i=0,vnewid.size-1 if(!vid.contains(vnewid.x(i))) vtmp.append(vnewid.x(i))
    vtmp.resize(sz)
    if(vtmp.size) {
      vnewdel.resize(vtmp.size)
      rd.uniform(delm[prty][poty]-deld[prty][poty],delm[prty][poty]+deld[prty][poty])
      vnewdel.setrand(rd)
      ce.o(id).setdvi(vtmp,vnewdel,2)
    }
  }
  dealloc(a)
  return 1
}

//** gethublims(col,hubtype,hubfactor,numhubs,mode) 
// get a matrix of size CTYPi X CTYPi, specifying div with mat.x(hubtype,othertype)
// and conv with mat.x(othertype,hubtype)
// hubtype = type of hub. hubfactor = desired ratio of hub div/conv vs non-hub div/conv
// numhubs = # of hubs. col = COLUMN for which to set hubs.
// mode == 0 <-- hub div(conv) is set to hubfactor * original div(conv)
// mode == 1 <-- hub div(conv) is set so that final hub div = hubfactor * final non_hub div (same for conv)
//  formula is based on:  m / ((N-H*m) / (C-H)) = F , and then solving for m
//   m = div for the hubs,  F = desired ratio of final hub div to final non-hub div
//   N = # of synapses (links),  C = total # of postsynaptic cells (including hubs) , H = # of hubs
//   similarly done for conv , but replace N with appropriate values
//   (/u/samn/intfcol/notebook.dol_1:21933)
obfunc gethublims () { local ct,mode,from,to,lim,nc,nhubs,fctr localobj col,mat
  {col=$o1 ct=$2 fctr=$3 nhubs=$4 mode=$5 mat=new Matrix(CTYPi,CTYPi)}
  for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to]) {
    {nc=col.numc[to] if(ct==to)nc-=1} // deduct for self-link
    if(mode==0) {
      lim = int( 0.5 + col.div[ct][to]*fctr )
    } else {
      lim = int( 0.5 + col.div[ct][to]*col.numc[ct]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) )
    }
    mat.x(ct,to) = MINxy(lim, nc) // at most div to all postsynaptic cells
  }
  for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct]) {
    {nc=col.numc[from] if(ct==from)nc-=1} // deduct for self-link
    if(mode==0) {
      lim = int( 0.5 + col.conv[from][ct]*fctr )
    } else {
      lim = int( 0.5 + col.div[from][ct]*col.numc[from]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) )
    }
    mat.x(from,ct) = MAXxy(MINxy(lim, nc),1) // at most conv from all presynaptic cells, but at least 1
  }
  return mat
}

//** addhubs(column,cell-type,numhubs,scaling factor,skipI[,seed,allowz,hubmode,verbose])
// add hubs to the network by stealing wires from other neurons
// $o1 == column
// $2 == cell type of hub
// $3 == number of hubs to add
// $4 == scaling factor (should be > 1.0) for conv,div of hub
// $5 == skip div/conv of I cells
// $6 == seed - optional
// $7 == allowz - whether to allow pulling all links from/to another cell
// $8 == hubmode - which mode to use for gethublims (see above)
// $9 == verbose - optional
// function returns a Vector containing the ids of the cells selected as hubs (within column ids)
obfunc addhubs () { local a,ct,fctr,nhubs,idx,jdx,lseed,hubid,szorig,cursz,preid,poid,lim,skipI,to,from,vrb,changed,allowz,hmode\
                 localobj col,ce,vin,vout,nq,vd,vc,vdd,vdt,vddt,vpicked,vhubid,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz,mhlim
  col=$o1 ct=$2 nhubs=$3 fctr=$4 skipI=$5
  if(numarg()>5) lseed=$6 else lseed=1234
  if(numarg()>6) allowz=$7 else allowz=1
  if(numarg()>7) hmode=$8 else hmode=0
  if(numarg()>8) vrb=$9 else vrb=0
  {ce=col.ce hashseed_stats(lseed,lseed,lseed)}
  a=allocvecs(vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz)
  vrsz(col.allcells,vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vdsz,vcsz,vtmp)
  mhlim=gethublims(col,ct,fctr,nhubs,hmode)
  //vin,vout = input/output markers. vd,vc = div/conv.
  //vdd div/conv delays, vdt div/conv temp. vddt=div/conv delay temp
  //vpicked=which cells already picked as hubs
  vhubid=new Vector()
  {vhubid.indgen(col.ix[ct],col.ixe[ct],1) vhubid.shuffle() vhubid.resize(nhubs)}
  if(vrb) vlk(vhubid)
  for idx=0,vhubid.size-1 vpicked.x(vhubid.x(idx))=1 
  for idx=0,vhubid.size-1 { hubid=vhubid.x(idx) 
    if(vrb) printf("hub%d id = %d\n",idx+1,hubid)
    {ce.o(hubid).getdvi(vd,vdd,vw1,vw2,vprob,vsyn) ce.o(hubid).getconv(vc)}//IDs of post/presynaptic cells
    {ce.o(hubid).getconv(1.2,vcsz) vdsz.resize(CTYPi) vdsz.fill(0)}//counts of post/pre types
    for jdx=0,vd.size-1 vdsz.x(ce.o(vd.x(jdx)).type)+=1
    {vout.fill(0) vin.fill(0)}     //init as 0
    for jdx=0,vd.size-1 vout.x(vd.x(jdx))=1 //mark current postsynaptic cells
    for jdx=0,vc.size-1 vin.x(vc.x(jdx))=1  //mark current presynaptic cells
    for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to] && (!skipI || !ice(to))) {
      cursz=szorig=vdsz.x(to) // update divergence
      if(vrb) print "\torig div -> " , CTYP.o(to).s, " = " , szorig
      {lim=mhlim.x(ct,to) changed=1}
      while(cursz<lim && changed==1) { changed=0
        for(preid=col.ix[ct];preid<=col.ixe[ct] && cursz<lim;preid+=1) {// pick same presynaptic type
          if(vpicked.x(preid)) continue //dont take from other hubs
          ce.o(preid).getdvi(vdt,vddt,vw1,vw2,vprob,vsynt) 
          vtmp.fill(0)
          for jdx=0,vdt.size-1 vtmp.x(ce.o(vdt.x(jdx)).type)+=1
          if(!allowz && vtmp.x(to)<=1)continue//dont want to turn div of another cell to 0
          for jdx=0,vdt.size-1 { poid=vdt.x(jdx) // go thru postsynaptic cells looking for target type            
            if(ce.o(poid).type==to && poid!=hubid && vout.x(poid)==0) { cursz+=1
              {vd.append(poid) vdd.append(vddt.x(jdx)) vsyn.append(vsynt.x(jdx))}
              {vdt.remove(jdx) vddt.remove(jdx) vsynt.remove(jdx)}
              ce.o(preid).setdvi(vdt,vddt,vsynt) // update presynaptic cell
              vout.x(poid)=changed=1 // this cell synapses on poid
              break
            }
          }
        }
      }
      if(vrb) print "\tnew div -> " , CTYP.o(to).s, " = " , cursz
    }
    ce.o(hubid).setdvi(vd,vdd,vsyn) // update hub dvi
    for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct] && (!skipI || !ice(from))) {
      cursz=szorig=vcsz.x(from) // update convergence
      {lim=mhlim.x(from,ct) changed=1}
      if(vrb) print "\torig conv <- ", CTYP.o(from).s, " = " , szorig
      while(cursz<lim && changed==1) { changed=0
        for(preid=col.ix[from];preid<=col.ixe[from]&&cursz<lim;preid+=1) {
          if(preid==hubid || vin.x(preid)) continue // don't make self or double-connects
          ce.o(preid).getdvi(vdt,vddt,vw1,vw2,vprob,vsynt)
          for jdx=0,vdt.size-1{
            poid = vdt.x(jdx)
            if(vpicked.x(poid)) continue // don't take wires from other hubs            
            if(ce.o(poid).type==ct){ ce.o(poid).getconv(1.2,vtmp)
              if(allowz || vtmp.x(from)>1) { // make sure not to remove all inputs of a type to a cell
                vdt.x( jdx ) = hubid // reassign input to hub
                ce.o(preid).setdvi(vdt,vddt,vsynt) // reset presynaptic cell's div
                vin.x( preid ) = changed = 1 // mark input
                cursz += 1
                break
              }
            }
          }
        }
      }
      if(vrb) print "\tnew conv <- " , CTYP.o(from).s, " = " , cursz
    }    
  }
  {dealloc(a) return vhubid}
}

//* mkcolnqs - make an nqs with current pmat,wmat,delm,deld info for use by a COLUMN for wiring
// "dist" represents distance between columns: dist==0 for intra-COLUMN setup, dist>0 for INTER-COLUMN setup
proc mkcolnqs () { local from,to,sy,idx,d localobj froms,tos,sys
  if(numarg()>0)idx=$1 else idx=0
  {nqsdel(colnq[idx]) colnq[idx]=new NQS("froms","tos","sys","from","to","sy","w","pij","delm","deld","loc","dist")}
  colnq[idx].strdec("froms","tos","sys")
  for from=0,CTYPi-1 { froms=CTYP.o(from)
    for to=0,CTYPi-1 { tos=CTYP.o(to)
      for d=0,colr if(pmat[from][to][d]>0) for sy=0,STYPi-1 if(wmat[from][to][sy][d]>0) { sys=STYP.o(sy)
        colnq[idx].append(froms.s,tos.s,sys.s,from,to,sy,wmat[from][to][sy][d],pmat[from][to][d],delm[from][to],deld[from][to],synloc[from][to],d)
      }
    }
  }
}

//* mkcols - make the COLUMNs
proc mkcols () { local id,x,y,seed
  id=0
  for y=0,colH-1 for x=0,colW-1 {
    if(dbgcols)seed=dvseed else seed=(id+1)*dvseed
    lcol.append(gcol[y][x]=new COLUMN(id,vcpercol,colnq,seed,x,y,setdviPT))
    col[id]=gcol[y][x]
    col[id].verbose=verbose_INTF6
    id+=1
  }
}

//* wirecols - setup inter-COLUMN connectivity with NetCon
proc wirecols () { local x1,y1,x2,y2,dx,dy,maxd,d localobj fromc,toc
  if(numarg()>0) d=$1 else d=colr
  if(torus) { // wraparound
    //alternate coordinates: ( -colW+x   ,  -colH+y )
    //alternate system: -5  -4  -3  -2  -1
    //original system:   0   1   2   3   4
    //layed out as a line: -5  -4  -3  -2  -1  0   1   2   3   4
    //only need to compare in normal system, and 1 alternate coordinate vs original (and vice versa)
    for y1=0,colH-1 for x1=0,colW-1 for y2=0,colH-1 for x2=0,colW-1 {
      if(y1==y2 && x1==x2) continue // skip self-self    
      dx=MINxy(abs(x1-x2), MINxy(abs((-colW+x1)-x2), abs(x1-(-colW+x2))) )
      dy=MINxy(abs(y1-y2), MINxy(abs((-colH+y1)-y2), abs(y1-(-colH+y2))) )
      if((maxd=MAXxy(dx,dy)) > d) continue // skip too far
      gcol[y1][x1].wire2col(gcol[y2][x2],colnq,maxd,ncl) // unidirectional wiring
    }
  } else { // no wrap-around
    for y1=0,colH-1 for x1=0,colW-1 for y2=0,colH-1 for x2=0,colW-1 {
      if(y1==y2 && x1==x2) continue // skip self-self    
      if((maxd=MAXxy(abs(x1-x2),abs(y1-y2))) > d) continue // skip too far
      gcol[y1][x1].wire2col(gcol[y2][x2],colnq,maxd,ncl) // unidirectional wiring
    }
  }
}

//* intercoloff - turn off all weights between COLUMNs
proc intercoloff () { local i localobj xo
  for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) {
    for i=0,6 xo.weight(i)=0
  }
}

//* intercolmul(from,to,sy,w)
proc intercolsyw () { local from,to,sy,w localobj xo
  setwmat()
  for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) {
    if(xo.pre.flag("inhib")) {
      for case(&sy,GA,GA2) xo.weight(sy) = wmat[xo.pre.type][xo.syn.type][sy][1]
    } else {
      for case(&sy,AM2,NM2) xo.weight(sy) = wmat[xo.pre.type][xo.syn.type][sy][1]
    }
  }
}

//* function calls to setup network

//** # of cells per column
setcpercol() //new numbers (10aug30)

//** setup pmat
if(name_declared("nqpmat")==2) { // read pmat from NQS if available, else set to default
  if(nqpmat!=nil) nq2pmat(nqpmat) else setpmat()
} else setpmat()
if(pmatscale!=1) scalepmat(pmatscale)

//** setup synapse locations,delays,wmat
setsynloc()
setdelmats()
setwmat() // new KMJ version

scrsz=50*1e3
double scr[scrsz]

//** make cells, columns, wire columns
mkcolnqs()
mkcols()
wirecols(1)

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)