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: network.hoc,v 1.126 2010/12/29 03:14:37 cliffk Exp $

print "Loading network.hoc..."

//* Numbers and connectivity params

{declare("colW",1,"colH",1,"torus",1)}
{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",15,"IEGain",4*15/11,"IIGain",4*15/11)}
{declare("NMAMR",0.1,"EENMGain",1,"EIGainInterC",0.125,"EEGainInterC",0.325*0.5)}

//* 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",1)}//16//8//4
{declare("pmatscale",1/scale)} // scale for pmat - originally 1/scale - 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

//* 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 
//  pmat[E2][E2][1]=0.14 // turned off 9/23/10
  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
  pmat[E5B][E5B][1]=0.25
  pmat[E5B][E5B][2]=0.1
  pmat[E5B][E5R][0]=0.017
  pmat[E5B][E5R][1]=0.25
  pmat[E5B][E5R][2]=0.1
  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
  pmat[E5R][E5B][1]=0.25
  pmat[E5R][E5R][0]=0.191
  pmat[E5R][E5R][1]=0.14
  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
  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
  wmat[E5B][E5B][AM2][2]=0.47 * EEGainInterC
  wmat[E5B][E5R][AM2][0]=0.24
  wmat[E5B][E5R][AM2][1]=0.47 * EEGainInterC
  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
  wmat[E5R][E5R][AM2][0]=0.66
  wmat[E5R][E5R][AM2][1]=0.47 * EEGainInterC
  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")
  }
}

//* clear sprmat entries to 0
proc clrsprmat () { local i,j
  for i=0,CTYPi for j=0,CTYPi sprmat[i][j]=0
}

//* unkill/prune all cells
proc unkp () {
  for i=0,allcells-1 {
    ce.o(i).flag("dead",0)
    ce.o(i).prune(0)
  }
}

//* kill cells who's ids are in $o1
proc dokill () { local id
  for vtr(&id,$o1) ce.o(id).flag("dead",1)
}

//* 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
  }
}

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




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)