Emergence of physiological oscillation frequencies in neocortex simulations (Neymotin et al. 2011)

 Download zip file 
Help downloading and running models
"Coordination of neocortical oscillations has been hypothesized to underlie the "binding" essential to cognitive function. However, the mechanisms that generate neocortical oscillations in physiological frequency bands remain unknown. We hypothesized that interlaminar relations in neocortex would provide multiple intermediate loops that would play particular roles in generating oscillations, adding different dynamics to the network. We simulated networks from sensory neocortex using 9 columns of event-driven rule-based neurons wired according to anatomical data and driven with random white-noise synaptic inputs. ..."
1 . 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 [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 L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell;
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA; Gaba;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Oscillations; Synchronization; Laminar Connectivity;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; GabaA; AMPA; NMDA; Gaba; Gaba; Glutamate;
misc.mod *
nstim.mod *
stats.mod *
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
grvec.hoc *
init.hoc *
labels.hoc *
local.hoc *
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc *
python.hoc *
pywrap.hoc *
simctrl.hoc *
spkts.hoc *
stats.hoc *
syncode.hoc *
xgetargs.hoc *
// $Id: col.hoc,v 1.83 2011/03/23 15:36:33 samn Exp $ 

//* globals
declare("UNIF",0)    //use roughly uniform prob. distrib for div
declare("NEGEXP",1)  //use neg exp. with mean=div prob. distrib in setdvi
declare("PARETO",2)  //use pareto (scalefree) distrib with mean=div prob., shape=2
declare("NORM",3)    //use normal distrib with mean=div, stdev=div/5
declare("DVFIXED",4) //use fixed divergence ( no variability )
declare("PARETOSH",4.5)//shape for pareto (scalefree) distrib -- low value will produce fatter/higher tails

//* template CSTIM
begintemplate CSTIM

double wmatex[1][1],ratex[1][1],EIBalance[1]
public wmatex,ratex,setwm,vq,col,sgrhzdel,sgrdur,setspks,pushspks,savewt,restorewt,turnoff,mulwts,sgrdel,EIBalance,shock
external CTYP,CTYPi,STYPi,AM,AM2,NM,NM2,GA,GA2,case,nil,allocvecs,dealloc,nqsdel
objref this,vq,col,vsgrpp,vsgrsidx,vwt
double iseed[1],sgrhzdel[1],sgrdur[1],sgrdel[1],usens[1]
objref ncl,nsl,rsl,vrse // <-- objects for external inputs via NetCon(ncl)+NetStim(nsl) 
// with Random(rsl), Random sead(vrse)
public ncl,nsl,rsl,vrse,usens,sgrcells,initrands

//* clrwm - reset wmatex to 0
proc clrwm () { local ct,sy
  for ct=0,CTYPi-1 for sy=0,STYPi-1 wmatex[ct][sy]=0
//* new CSTIM(COLUMN,iseed,sgrdur[,sgrhzdel,EIBalance,usens])
// iseed = seed for random # generator, sgrdur = duration of persistent stims
// sgrhzdel = variability of stim rates, EIBalance = make E,I inputs balanced
// usens = whether to use NetStim instead of vq for persistent stims
proc init () { 
  double wmatex[CTYPi][STYPi]
  double ratex[CTYPi][STYPi]
  double EIBalance[1]
  clrwm() col=$o1 col.cstim=this
  if(numarg()>1) iseed=$2 else iseed=1234
  sgrdur=$3 sgrdel=0
  vsgrpp=new Vector(CTYPi) //% X 100 of cells of a type to stim, used in stim,sgrcells
  vsgrsidx=new Vector(CTYPi) //startind index of cell to stim when using topstim
  vq=new NQS("ind","time","wt","sy")
  if(numarg()>3)sgrhzdel=$4 else sgrhzdel=0.2
  if(numarg()>4)EIBalance=$5 else EIBalance=0
  if(numarg()>5)usens=$6 else usens=0
//** setwm - set weights of external inputs to INTF6s
proc setwm () {  local i,sz,ct,sy localobj nq,vct,vsy,vw,vr
  {nq=$o1 nq.tog("DB") vct=nq.getcol("ct")}
  {vsy=nq.getcol("sy") vw=nq.getcol("w") vr=nq.getcol("rate")}
  for ct=0,CTYPi-1 for sy=0,STYPi-1 wmatex[ct][sy]=ratex[ct][sy]=0
  for i=0,nq.v.size-1 {
//** setspks - set time of the external inputs to cells, can still call shock later
proc setspks () { local se,ndx localobj rdm,ce,intf
  if(ce.count<=0) { //need INTFs to stim
    printf("CSTIM: no cells to stim!\n")
  intf=ce.o(0) rdm=new Random() rdm.ACG(iseed) //initialize Random object
  if(vq==nil) vq=new NQS("ind","time","wt","sy") else vq.clear() //NQS storing INTF6 spike times
  intf.clrvspks() //clear INTF6 spike times
  if(sgrdur>0) {
    if(usens) nsstim(rdm) else popstim(rdm) //add random firing of INTF6 cells
  // for ndx=1,MAXNP-1 if(npulse[ndx]!=0) pulsestim(rdm,ndx) //add pulse-like firing of INTF6 cells
  // if(usevpl && VPLGain>0) vplstim(rdm) // VPL inputs
  if(vq.v.size>0) pushspks()
//** pushspks() - finalize spike times to INTF6 cells
// (useful, since after calls to shock/pulsestim, vq is modified)
proc pushspks () { localobj ce,intf
  ce=col.ce intf=ce.o(0)
  vq.sort("time") vq.sort("ind") //finalize INTF6 spike times
  intf.initvspks(vq.v,vq.v[1],vq.v[2],vq.v[3]) // single global call
//** savewt - initialize temporary storage space to allow modifications to vq weight column
proc savewt () {
  if(vwt==nil) vwt=new Vector(vq.v.size)  
//** restorewt - restore weights in vq to original values
proc restorewt () {
  if(vwt==nil) return
//** mulwts(celltype,synapsetype,scaling factor,[skippush]) - scale weights of cell and synapse type by factor
proc mulwts () { local ct,sy,fctr,i,j localobj ce
  ct=$1 sy=$2 fctr=$3 ce=col.ce vq.select(-1,"sy",sy)
  for i=0,vq.ind.size-1 {j=vq.ind.x(i) 
    if(ce.o(vq.v[0].x(j)).type==ct) vq.v[2].x(j)*=fctr
  if(numarg()>3) { if(!$4) pushspks() } else pushspks()
//** turnoff(celltype,synapsetype,[skippush]) - turn off spikes by setting weights to 0
proc turnoff () { 
  if(numarg()>2) mulwts($1,$2,0,$3) else mulwts($1,$2,0)
//** pulsestim(rate[,maxt,clear,seed]) - stim using pulses @ specified rate - doesn't work correctly now
proc pulsestim () { local nspk,tt,i,j,nrate,ti,maxt localobj rdp,vqtmp
  rdp=new Random()   nrate=$1
  if(numarg()>1) maxt=$2 else maxt=tstop
  if(numarg()>2) {
    if($3) vq.clear() 
  } else vq.clear()  
  if(numarg()>3) rdp.ACG($4) else rdp.ACG(1234)
  vqtmp=new NQS()
  tt=1 rdp.poisson(1000/nrate)
  while(tt<=maxt) {
    while(ti <= 0) ti=rdp.repick()    
    if(vqtmp.size(-1)==0) vqtmp.cp(vq) else vqtmp.append(vq)
    print tt
    tt += ti
  vq.cp(vqtmp) nqsdel(vqtmp)
  vq.sort("time") vq.sort("ind")

//** initrands() - initialize Random # objects used by NetCon/NetStim for external inputs
// this function should be called in init -- at the start of each run to ensure same random # stream each time
proc initrands () { local i,se localobj rds
  if(vrse==nil) return  
  for i=0,vrse.size-1 {
    rds = rsl.o(i)       // go thru Random objects 
    se = vrse.x(i)       // seed previously used
    rds.MCellRan4(se,se) // initialize MCellRan4 Random # generator
    rds.negexp(1)        // must be called this way 

//** nsstim(Random) - setup random external inputs via NetStims
proc nsstim () { local a,sy,sdx,idx,sglo,sghi,hzlo,hzhi,ct,se localobj rdm,vsy,vhz,nc,ns,rds
  rdm=$o1  vsy.append(GA) vsy.append(GA2) vsy.append(AM2) vsy.append(NM2)
  sglo=1-sgrhzdel sghi=1+sgrhzdel if(sglo<0) sglo=0
  {nsl=new List() ncl=new List() rsl=new List() vrse=new Vector() se=10}
  for sdx=0,vsy.size-1 { sy=vsy.x(sdx) // go through synapse types
    for ct=0,CTYPi-1 if(col.numc[ct] && wmatex[ct][sy] && ratex[ct][sy]) {//go through cell types, check weights,rates of inputs
      vhz.resize(col.numc[ct]) //set vector to # of cells of type ct
      {hzlo=MAXxy(ratex[ct][sy]*sglo,1) hzhi=MAXxy(ratex[ct][sy]*sghi,1)}//find min,max frequencies
      rdm.discunif(hzlo,hzhi) //set random # generator to frequency btwn hzlo,hzi
      vhz.setrand(rdm) //set rate of inputs to each synapse
      for idx=col.ix[ct],col.ixe[ct] { // for each cell of type ct (idx is id of the cell)

        nsl.append( ns = new NetStim(0.5) ) // NetStim is source
        ns.number = int( 0.5 + vhz.x(idx-col.ix[ct]) * sgrdur / 1e3 ) // # of inputs
        ns.start = sgrdel // approx start of first input
        ns.noise = 1 // make it fully random
        ns.interval = 1e3 / vhz.x(idx-col.ix[ct]) // avg interval btwn the inputs

        rsl.append( rds = new Random() )
        rds.negexp(1) // set random # generator using negexp(1) - avg interval in NetStim.interval
        rds.MCellRan4(se,se)  // seeds are in order, shouldn't matter
        ns.noiseFromRandom(rds) // use random # generator for this NetStim
        vrse.append(se) // save random # seed

        ncl.append( nc = new NetCon(ns,col.ce.o(idx)) ) // set the connection (netstim->postsynaptic INTF6 cell)
        nc.delay = 0
        nc.weight(sy) = wmatex[ct][sy]

        se += 10 // increment seed for NetStim random # generator

//** popstim(Random) - setup random external synaptic inputs 
proc popstim () { local a,se,sy,sz,sdx,idx,jdx,pos,sglo,sghi,hzlo,hzhi,nspk,ct,mxr localobj rdm,vsy,vnspks,vtmp,vfctr
  rdm=$o1  vsy.append(GA) vsy.append(GA2) vsy.append(AM2) vsy.append(NM2)
  pos=mxr=0 sglo=1-sgrhzdel sghi=1+sgrhzdel if(sglo<0) sglo=0
  for ct=0,CTYPi-1 for sy=0,STYPi-1 mxr=MAXxy(mxr,ratex[ct][sy])
  if(EIBalance) { vfctr.resize(col.allcells) // flag for excit/inhib inputs to be balanced
    {rdm.uniform(sglo,sghi) vfctr.setrand(rdm) jdx=pos=0}
    for idx=0,col.allcells-1 { ct=col.ce.o(idx).type
      for sdx=0,vsy.size-1 { sy=vsy.x(sdx)
        if(wmatex[ct][sy] && ratex[ct][sy]) {
          if((nspk=int(ratex[ct][sy]*vfctr.x(idx)*sgrdur/1e3))<1) {nspk=1}
          if(sy==GA||sy==GA2) {
          } else vq.v[2].fill(wmatex[ct][sy],pos,pos+nspk-1)
          pos += nspk
  } else { // default here
    for sdx=0,vsy.size-1 { sy=vsy.x(sdx) // go through synapse types
      for ct=0,CTYPi-1 if(col.numc[ct] && wmatex[ct][sy] && ratex[ct][sy]) {//go through cell types, check weights,rates of inputs
        vnspks.resize(col.numc[ct]) //set vector to # of cells of type ct
        {hzlo=MAXxy(ratex[ct][sy]*sglo,1) hzhi=MAXxy(ratex[ct][sy]*sghi,1)}//find min,max frequencies
        rdm.discunif(hzlo,hzhi) //set random # generator to frequency btwn hzlo,hzi
        {vnspks.setrand(rdm) vnspks.mul(sgrdur/1e3) vnspks.apply("int")} //set # of spikes per intf in vnspks
        jdx = 0 // pointer into vnspks
        for idx=col.ix[ct],col.ixe[ct] { // for each cell of type ct (idx is id of the cell)
          vq.v[0].fill(idx,pos,pos+vnspks.x(jdx)-1) //fill vq.v[0] with ID of cell
          if(sy==GA||sy==GA2) { //is synapse inhibitory? if so, use negative weights
          } else vq.v[2].fill(wmatex[ct][sy],pos,pos+vnspks.x(jdx)-1) //use positive weights for excit. synapses
          vq.v[3].fill(sy,pos,pos+vnspks.x(jdx)-1)//fill corresponding entries in vq.v[3] with synapse type
          pos += vnspks.x(jdx) //increment pointer into vq row
          jdx += 1 //move jdx to next cell
  vq.pad(pos)//set all vq vectors to right size
  rdm.uniform(sgrdel,sgrdel+sgrdur) vq.v[1].setrand(rdm) // set times of inputs in vq.v[1]

//** shock(duration,percent,type[,time,seed,sy,wt])
// shock random % of cells
// single shock so duration means it spreads it out over eg 5 ms in terms of delivery to different cells
// % > 100 means multiple stime eg 200% will stim each cell 2 times in that interval
// sy specifies which synapse to stim, default is AM2
// wt specifies weight of inputs -- default is 1e9 for guaranteed spiking @ AM2 synapses. otherwise
// will provide a wt-dependent stim
proc shock () { local ct,startt,sy,wt
  if(numarg()>3) startt=$4 else startt=0
  if(numarg()>4) mc4seed_stats($5) else mc4seed_stats(5974302451) 
  if(numarg()==7) {
    {sy=$6    wt=$7   sgrcells($1,$2,$3,startt,1,0,sy,wt)}
  } else sgrcells($1,$2,$3,startt)

//* sgrcells(dur,percent,stimcelltype,starttime[,burst#,ISI,sy,wt]) - stim the cells
//  sgrcells(dur,percent,stimcelltype,starttime[,burst#,ISImin,ISImax])
// if numarg == 8, uses specified sy,wt for sub-threshold stim
proc sgrcells () { local ii,sgrdur,sgrpp,sgrpsz,stimcc,startt,a,bst,sy,wt localobj v1,v2,v3,vqtmp,dnq,vid,str
  if (vq==nil) vq=new NQS("ind","time","wt","sy")
  sgrdur=$1 sgrpp=$2 stimcc=$3 startt=$4 str=new String()
  bst=1 v3.resize(bst) // number of spikes/cell -- default 1
  if (argtype(5)==0) { 
    bst=$5 v3.resize(bst) 
    if (argtype(7)==0) { v3.setrnd(4.5, $6, $7) // interval of intervals
    } else if (argtype(6)==0) { v3.fill($6)
    } else v3.fill(4) // default 250 Hz without randomization
  if(numarg()==8) {sy=$7 wt=$8 } else {sy=AM2 wt=1e9} // use specified wt,sy?
  v3.x[0]=0.0 // so first one delivered at start time
  vqtmp=new NQS("ind","time","wt","sy")
      v1.setrnd(6,col.ix[stimcc],col.ixe[stimcc])//unique random ints in range
    } else {
      v1.setrnd(5,col.numc[stimcc])      v1.add(col.ix[stimcc])
  } else {
    if(col.div[stimcc][stimcc][0]>0)sprint(str.s,"to%s",CTYP.o(stimcc).s) else str.s="toE4"
    {dnq.select("type",stimcc) dnq.sort(str.s) vid=dnq.getcol("id")}
    } else v1.copy(vid,0,sgrpsz-1)
  v2.resize(v1.size) // same size as v1
  v2.setrnd(4.5,startt,startt+sgrdur)// at some random time between startt and startt+sgrdur    
  for ii=1,bst {vqtmp.v[0].append(v1) vqtmp.v[1].append(v2.add(v3.x[ii-1]))}
  {vqtmp.pad() vqtmp.v[2].fill(wt) vqtmp.v[3].fill(sy) } 
  vq.append(vqtmp) // need to call pushspks after call to shock/sgrcells

endtemplate CSTIM

//* template COLUMN
begintemplate COLUMN

external CTYPi,STYPi,ice,DEND,SOMA,IsLTS,CTYP,STYP,allocvecs,vrsz,dealloc,nil,NM,NM2,AM,AM2,GA,GA2,vlk,nqsdel

double numc[50]
double ix[50],ixe[50] //start,end IDs of types in a column
double div[50][50]//div[i][j]==# of outputs from type i->j
double wmat[50][50][25] // wmat[i][j][k]==weight from type i->j for synapse k
double wd0[50][50][25] // wmat[i][j][k]==weight from type i->j for synapse k
double delm[50][50]//avg. delay from type i->j
double deld[50][50]//delay variance from type i->j
double conv[50][50]
double pmat[50][50]
double synloc[50][50]//location of synapses
double allcells[1],ecells[1],icells[1],setdviPT[1],dgcor[1],xpos[1],ypos[1],wseed[1],verbose[1]
double syty1[50][50] //synapse type lookup table
double syty2[50][50] //synapse type lookup table

public id,wmat,wd0,pmat,wire,div,delm,deld,conv,numc,setdviPT,allcells,ecells,icells,xpos,ypos,wire2col
public setwmat,setpmat,setdelmats,setsynloc,ix,ixe,cstim,defsynloc,wseed,verbose,syty1,syty2
public mknetnqss
objref this,ce,intf,cstim
objref cellsnq,connsnq,xcolconnsnq
strdef name
public ce,intf,allwtsoff,allwtson,inhibon,inhiboff,name
public ctt,cttr,stt
public cellsnq,connsnq,xcolconnsnq,mkcellsnq

proc setstim () {

//** allwtson -- turn all internal weights on
proc allwtson () { local i,j,s
  for i=0,CTYPi-1 for j=0,CTYPi-1 for s=0,STYPi-1 if(wmat[i][j][s]>0) wd0[i][j][s]=1

//** allwtsoff -- turn all internal weights off
proc allwtsoff () { local i,j,s
  for i=0,CTYPi-1 for j=0,CTYPi-1 for s=0,STYPi-1 wd0[i][j][s]=0

// ** inhiboff -- turn inhibition off -- only call after setting up network
proc inhiboff () { local i,from,to
  for i=0,allcells-1 if(ice(ce.o(i).type)) ce.o(i).prune(1) //prune all inhib cells
  for from=0,CTYPi-1 if(ice(from)) for to=0,CTYPi-1 wd0[from][to][GA]=0

// ** inhibon -- turn inhibition on
proc inhibon () { local i,from,to
  for i=0,allcells-1 if(ice(ce.o(i).type)) ce.o(i).prune(0) //unprune all inhib cells
  for from=0,CTYPi-1 if(ice(from)) for to=0,CTYPi-1 wd0[from][to][GA]=1

//** getNdv - gets # of outputs based on divergence settings, always returns at least 1
// $o1 = random object, $2 = avg. # of outputs
func getNdv () { local dv,N localobj rdm
  rdm=$o1 dv=$2
  if(setdviPT==PARETO) {
    rdm.avg=dv N=int(rdm.pick+0.5)
  } else if(setdviPT==NEGEXP) {
  } else if(setdviPT==UNIF) {
  } else if(setdviPT==NORM) {
  } else if(setdviPT==DVFIXED) {
    return dv
  return N

//** defsynloc - set default values in synloc
proc defsynloc () { 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

//** setsynloc - sets synapse location - not modifiable for now
proc setsynloc () { local i,sz,from,to,sy localobj nq,vfrom,vto,vloc
  if(numarg()==0) {
  } else {
    {nq=$o1 nq.verbose=0 nq.select("dist",0)}
    {vfrom=nq.getcol("from") vto=nq.getcol("to") vloc=nq.getcol("loc")}
    for i=0,sz-1 synloc[vfrom.x(i)][vto.x(i)]=vloc.x(i)
    {nq.tog("DB") nq.verbose=1}

//** prdiv - prints div info
proc prdiv () { local i,j,k
  for i=0,CTYPi-1 for j=0,CTYPi-1 for k=0,STYPi-1 if(wmat[i][j][k]>0) if(div[i][j]) {

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

//** clrpmat - sets pmat entries to 0
proc clrpmat () { local from,to
  for from=0,CTYPi-1 for to=0,CTYPi-1 pmat[from][to]=0 // clear it

//** setpmat(nqs with from,to,pij columns) - from=presynaptic type,
// to=postsynaptic type, pij=connection probability
proc setpmat () { local i,sz,from,to localobj nq,vfrom,vto,vpij
  {nq=$o1 nq.verbose=0 nq.select("dist",0)}
  clrpmat() // clear it
  {vfrom=nq.getcol("from") vto=nq.getcol("to") vpij=nq.getcol("pij")}
  for i=0,sz-1 pmat[vfrom.x(i)][vto.x(i)]=vpij.x(i)
  setdivmat() // set divergence vectors based on pmat
  {nq.tog("DB") nq.verbose=1}

//** setdivmat - sets div/conv based on values in pmat
proc setdivmat () { local from,to,cl
  for from=0,CTYPi-1 for to=0,CTYPi-1 if(pmat[from][to]) {
    div[from][to] =  ceilg(pmat[from][to]*numc[to])
    conv[from][to] = int(0.5 + pmat[from][to]*numc[from])

//** clrwmat - sets wmat entries to 0
proc clrwmat () { local from,to,sy
  for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=0,STYPi-1 wmat[from][to][sy]=0
//** clrwmd0 - sets wd0 entries to 0
proc clrwd0 () { local from,to,sy
  for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=0,STYPi-1 wd0[from][to][sy]=0

//** clrsyty - sets syty1,sty2 entries to -1
proc clrsyty () { local from,to
  for from=0,CTYPi-1 for to=0,CTYPi-1 syty1[from][to]=syty2[from][to]=-1

//** setwmat(nqs with from,to,sy,w columns) - from=presynaptic type,
// to=postsynaptic type, sy=synapse type, w=weight
proc setwmat () { local i,sz,from,to,sy localobj nqwm,vfrom,vto,vsy,vw
  {nqwm=$o1 nqwm.verbose=0 nqwm.select("dist",0) clrwmat() clrwd0()}
  {vfrom=nqwm.getcol("from") vto=nqwm.getcol("to")}
  {vsy=nqwm.getcol("sy") vw=nqwm.getcol("w") sz=vfrom.size}
  for i=0,sz-1 { 
    {sy=vsy.x(i) from=vfrom.x(i) to=vto.x(i)}
    {wmat[from][to][sy]=vw.x(i) wd0[from][to][sy]=1}
    if(sy==NM2) syty2[from][to]=sy else syty1[from][to]=sy
  {nqwm.tog("DB") nqwm.verbose=1}

//** setdelmats(nqs with from,to,delm,deld columns) - from=presynaptic type,
// to=postsynaptic type, sy=synapse type, delm=avg. delay, deld=variation in delay
proc setdelmats () {  local i,sz,from,to localobj nq,vfrom,vto,vdelm,vdeld
  {nq=$o1 nq.verbose=0 nq.select("dist",0)}
  {vfrom=nq.getcol("from") vto=nq.getcol("to") vdelm=nq.getcol("delm") vdeld=nq.getcol("deld")}
  for i=0,sz-1 {delm[vfrom.x(i)][vto.x(i)]=vdelm.x(i) deld[vfrom.x(i)][vto.x(i)]=vdeld.x(i)}
  {nq.tog("DB") nq.verbose=1}

//** clrnqss - clear nqs objects used for storing cell/connectivity info
proc clrnqss () {
  {nqsdel(cellsnq) cellsnq=nil}
  {nqsdel(connsnq) connsnq=nil}
  {nqsdel(xcolconnsnq)  xcolconnsnq=nil}

//** new COLUMN(ID,vnumc,[nqwiring,wiringseed,xpos,ypos,setdviPT,mknetnqss]])
proc init () {
  {id=$1 dgcor=0 setdviPT=DVFIXED wseed=1234 verbose=0}
  if(numarg()>=8) mknetnqss=$8 else mknetnqss=0
  double numc[CTYPi],ix[CTYPi],ixe[CTYPi],div[CTYPi][CTYPi]
  double wmat[CTYPi][CTYPi][STYPi],wd0[CTYPi][CTYPi][STYPi],delm[CTYPi][CTYPi],deld[CTYPi][CTYPi]
  double conv[CTYPi][CTYPi],pmat[CTYPi][CTYPi],synloc[CTYPi][CTYPi],syty1[CTYPi][CTYPi],syty2[CTYPi][CTYPi]
  {clrsyty() mkcells($o2) sprint(name,"C%d",id) clrnqss()}
  if(mknetnqss) mkcellsnq()  // make cellsnq and connsnq
  if(numarg()>2) {
    {setpmat($o3) setwmat($o3) setsynloc($o3) setdelmats($o3)}
    if(numarg()>=6) {xpos=$5 ypos=$6 sprint(name,"C%d_X%d_Y%d",id,xpos,ypos)} // position in 2D grid
    if(numarg()>=7) setdviPT=$7 // wiring scheme
    if(numarg()>3) wire(wseed=$4)

//** initdivrnd(seed) - initializes random object for wiring
obfunc initdivrnd () { local s localobj rd
  if(setdviPT==PARETO) rd=new rdmpareto(1,PARETOSH,s) else {rd=new Random() rd.ACG(s)}
  return rd

//** wire([seed]) - wires the column using div matrix
proc wire () { local x,y,z,ii,jj,a,del,prid,prty,poty,dv,dvt,lseed,dvrand,szo\
               localobj v1,v2,v3,v4,v5,v6,v7,vdvd,vidx,vdel,vdist,vwt1,vwt2,vtmp,opr,opo,st,l,rdm
  if(verbose) printf("wiring COLUMN %d",id)
  a=allocvecs(v1,v2,v3,v4,v5,v6,v7,vidx,vdel,vdvd,vtmp,vdist,vwt1,vwt2) z=0 l=new List() l.append(v1) l.append(v4)
  if (argtype(1)==0) lseed=$1 else lseed=1234
  if(dgcor) {  // div variability of 20% when dgcor==1
    dvrand=0.2 vdvd.setrnd(4,2*dvrand) vdvd.sub(dvrand) 
  rdm=initdivrnd(lseed)//initialize random # generator for use when dgcor==0
  // Create the connectivity NQS table.
  if(mknetnqss) {nqsdel(connsnq) connsnq=new NQS("id1","id2","del","dist","wt1","wt2")}
  for y=0,ce.count-1 { opr=ce.o(y)
    if(verbose) if(y%100==0)printf(".")
    prid=opr.id prty=opr.type
    for poty=0,CTYPi-1 if (numc[poty]!=0 && (dv=int(div[prty][poty]))>0) {      
      if(!dgcor) dv=getNdv(rdm,dv) else dv*=(1+vdvd.x[y])
      v3.setrnd(5,numc[poty]) v3.add(ix[poty]) // flag 5 provides integers
      if(prty==poty) {v1.cull(v3,prid) v3.copy(v1) v2.resize(v3.size)} // get rid of self connect      
      if( (cnt=v3.uniq(l,1)) > dv){ cnt = dv } // puts unsorted uniq vals into v4
      vrsz(cnt,v2,v4,v5,v6,v7) // v4 has poids
      v2.setrnd(4,2*deld[prty][poty]) // flag 4 provides doubles
      v2.abs() // no negative delays
      if(synloc[prty][poty]==DEND) v5.fill(1) else v5.fill(0)
      vidx.append(v4) vdel.append(v2) vdist.append(v5) // vidx.append(v1) will give sorted bug
      if(mknetnqss) {
        if(syty1[prty][poty]>=0) v6.fill(wmat[prty][poty][syty1[prty][poty]]) else v6.fill(0)
        if(syty2[prty][poty]>=0) v7.fill(wmat[prty][poty][syty2[prty][poty]]) else v7.fill(0)
        vwt1.append(v6) vwt2.append(v7)
    }//end poty
    if(vidx.size>0) {
      if(mknetnqss) for ii=0,vidx.size-1 connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
  }//end cell loop
  if(verbose) printf("\n")

//** wire2col(targetcol,nq,distance,ncl) - connect this COLUMN to another COLUMN
func wire2col () { local d,from,to,sz,dvm,dv,pij,dlym,dlyd,w1,w2,sy1,sy2,loc,i,idx,jdx,kdx,gid1,gid2,cnt,a\
                  localobj tcol,nq,ncl,vfrom,vto,vdelm,vdeld,vw,vsy,vpij,vloc,rdm,vidx,vdel,v1,v2,v3,v4,l,nc,prx,pox
  {tcol=$o1 nq=$o2 d=$3 ncl=$o4 nq.verbose=0}
  if(verbose) print "wiring COLUMN ", id, " to COLUMN ",tcol.id
  if(!nq.select("dist",d)) {nq.verbose=1 printf("wire2col WARNA: none found @ d=%d!\n",d) nq.tog("DB") return 0}
  a=allocvecs(vidx,vdel,v1,v2,v3,v4) l=new List() l.append(v1) l.append(v4)
  {vfrom=nq.getcol("from") vto=nq.getcol("to") vsy=nq.getcol("sy")}
  {vdelm=nq.getcol("delm") vdeld=nq.getcol("deld") vw=nq.getcol("w")}
  {vpij=nq.getcol("pij") vloc=nq.getcol("loc")}
  if(verbose) print "tcol.wseed = ",tcol.wseed, " vfrom.size = ",vfrom.size
  rdm=initdivrnd(tcol.wseed)//initialize random # generator
  if (mknetnqss && xcolconnsnq!=nil) xcolconnsnq=new NQS("id1","col2","id2","del","wt1","wt2")
  while(i<vfrom.size) {
    {from=vfrom.x(i) to=vto.x(i) dlym=vdelm.x(i) dlyd=vdeld.x(i)}
    {w1=vw.x(i) sy1=vsy.x(i) sy2=w2=-1 pij=vpij.x(i) loc=vloc.x(i)}
    dvm = ceilg(pij*tcol.numc[to]) // divergence
    if((sy1==AM || sy1==AM2) && i+1<vfrom.size) { // check for NMDA, which follows AMPA
      if(vfrom.x(i+1)==from && vto.x(i+1)==to && (vsy.x(i+1)==NM || vsy.x(i+1)==NM2)) {
        sy2=vsy.x(i+1) w2=vw.x(i+1)
    for idx=ix[from],ixe[from] { // go thru presynaptic cells
      hashseed_stats(ce.o(idx).gid,from,to,tcol.wseed) // why not call hashseed here?
      v3.setrnd(5,tcol.numc[to]) v3.add(tcol.ix[to])
      if(verbose && v3.size<1) print "dv=",dv," from ",CTYP.o(from).s, " to ",CTYP.o(to).s," numc=",tcol.numc[to]," ix=",tcol.ix[to]
      if((cnt=v3.uniq(l,1))>dv) cnt=dv 
      {vrsz(cnt,v2,v4) v2.setrnd(4,2*dlyd) v2.add(dlym-dlyd)}
      v2.abs() // no negative delays
      if(verbose>1) {
        print CTYP.o(from).s, " -> ", CTYP.o(to).s, ": pij ", pij, ", dvm ", dvm, ", dv ", dv, ", v4.size ", v4.size
        printf("v4: ") vlk(v4)
        printf("v2: ") vlk(v2)
      ce.o(idx).flag("out",1) // make sure NetCon events enabled from this cell
      for jdx=0,v4.size-1 {
        prx = ce.o(idx)
        pox = tcol.ce.o(v4.x(jdx))
        ncl.append(nc=new NetCon(prx,pox))
        if(sy2!=-1) nc.weight(sy2)=w2
        if(mknetnqss) {
          if(sy2!=-1) {
          } else {
    if(sy2!=-1) i+=2 else i+=1 // if used NMDA, skip it
  {nq.verbose=1 nq.tog("DB") dealloc(a) return 1}

//** setarrs(vector of num per type) - sets  up arrays/numc/cell counts
proc setarrs () { local ct,cnt,cl,ii localobj vnumc
  {vnumc=$o1 cnt=allcells=icells=ecells=0}
  for ct=0,CTYPi-1 if (vnumc.x(ct)>0) { numc[ct]=vnumc.x(ct)
    ix[ct]=cnt ixe[ct]=cnt+numc[ct]-1 cnt+=numc[ct]
    if(ice(ct))icells+=numc[ct] else ecells+=numc[ct]
  } else numc[ct]=0

//** mkcells(vector of num per type) - make the cells
proc mkcells () { local ct,idx,jdx,ty,nc,ic localobj vnumc,xo,st
  if(ce==nil) ce=new List()
  {vnumc=$o1 setarrs(vnumc)  }
  idx=0 // starting ID for cells in column - assumes all columns same size
  for ct=0,CTYPi-1 if(vnumc.x(ct)>0) { ic=ice(ct)
    for jdx=ix[ct],ixe[ct] {
      ce.append(xo=new INTF6(jdx,ct,ic,this.id))
  if(!ce.count) return  
  st=new String()
  if(verbose) print st.s
  {sprint(st.s,"%s.intf.flag(\"jcn\",1,1)",this) execute(st.s)}
  if(verbose) print " finished mkcells "

//** mkcellsnq() -- make NQS table for all of the cells in the column
proc mkcellsnq () { local ii localobj xo 
  cellsnq=new NQS("gid","id","col","ty","ice","xloc","yloc","zloc")
  for ii=0,ce.count - 1 {
    xo = ce.o(ii)

//*** ctt(&i) - iterate over cell types in increasing CTYPi order
iterator ctt () { local i
  for i=0,CTYPi-1 if(numc[i]>0) {
//*** cttr(&i) - iterate over cell types in decreasing CTYPi order
iterator cttr () { local i
  for(i=CTYPi-1;i>=0;i-=1) if(numc[i]>0) {
//*** stt(&i,&j,&k) - iterate over active synapses types in order
iterator stt () { local i,j,k
  for i=0,CTYPi-1 if(numc[i]>0) for j=0,CTYPi-1 if(numc[j]>0 && div[i][j]>0) {
    for k=0,STYPi-1 if(wmat[i][j][k]>0) {
      $&1=i $&2=j $&3=k

endtemplate COLUMN

Loading data, please wait...