Electrostimulation to reduce synaptic scaling driven progression of Alzheimers (Rowan et al. 2014)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:154096
"... As cells die and synapses lose their drive, remaining cells suffer an initial decrease in activity. Neuronal homeostatic synaptic scaling then provides a feedback mechanism to restore activity. ... The scaling mechanism increases the firing rates of remaining cells in the network to compensate for decreases in network activity. However, this effect can itself become a pathology, ... Here, we present a mechanistic explanation of how directed brain stimulation might be expected to slow AD progression based on computational simulations in a 470-neuron biomimetic model of a neocortical column. ... "
Reference:
1 . Rowan MS, Neymotin SA, Lytton WW (2014) Electrostimulation to reduce synaptic scaling driven progression of Alzheimer's disease. Front Comput Neurosci 8:39 [PubMed]
Citations  Citation Browser
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; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; Python;
Model Concept(s): Long-term Synaptic Plasticity; Aging/Alzheimer`s; Deep brain stimulation; Homeostasis;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org]; Rowan, Mark [m.s.rowan at cs.bham.ac.uk];
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; Glutamate;
/
RowanEtAl2014
batchscripts
mod
README
alz.hoc
alzinfo.m
autotune.hoc *
basestdp.hoc *
batch.hoc *
batch2.hoc *
batchcommon
checkirreg.hoc *
clusterrun.sh
col.dot *
col.hoc *
comppowspec.hoc *
condisconcellfig.hoc *
condisconpowfig.hoc *
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
e2hubsdisconpow.hoc *
e2incconpow.hoc *
filtutils.hoc *
flexinput.hoc
geom.hoc *
graphplug.hoc *
grvec.hoc *
infot.hoc *
init.hoc *
labels.hoc *
load.hoc *
local.hoc *
makepopspikenq.hoc *
matfftpowplug.hoc *
matpmtmplug.hoc *
matpmtmsubpopplug.hoc *
matspecplug.hoc *
mosinit.hoc
network.hoc *
nload.hoc *
nqpplug.hoc *
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc *
params.hoc
plot.py
plotavg.py
plotbatch.sh
plotbatchcluster.sh
plotdeletions.py
plotntes.py
powchgtest.hoc *
pyhoc.py
python.hoc *
pywrap.hoc *
ratlfp.dat *
redE2.hoc *
run.hoc
runsim.sh
setup.hoc *
shufmua.hoc *
sim.hoc
simctrl.hoc *
spkts.hoc *
stats.hoc *
syncode.hoc *
vsampenplug.hoc *
writedata.hoc
xgetargs.hoc *
                            
// $Id: network.hoc,v 1.205 2012/01/20 02:50:20 samn Exp $

//* 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())
declare("autotune",0,"useplast",1) // 
declare("randsy",-1) // whether to randomize wts (iff non-zero). iff > 0 use uniform distrib, if < 0 use normal
{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

//* swire & related column variables
declare("colside",30) // column diameter in micrometers
declare("swire",3) // whether to use 'spatial' wiring, 2 means use swirecut, 3 means use swirecutfl
declare("checkers",0) // whether to arrange cells in checkerboard pattern
declare("cxinc",3,"cyinc",3,"crad",1)
declare("EEsq",(colside/2)^2,"EIsq",(colside/2)^2,"IEsq",(colside/1)^2,"IIsq",(colside/2)^2) // params for swirecut
double dels[CTYPi][CTYPi]//[STYPi] //stdev of delays
double delv[CTYPi][CTYPi]//variance of delays
double vcond[CTYPi][CTYPi]//[STYPi] //conduction velocities
declare("layerzvar",25) // range of z location of cells in a layer (in micrometers)
declare("colxydist",50) // distance between adjacent columns (x,y)
declare("vlayerz",new Vector(CTYPi))
declare("pseed",4321)
declare("dvcond","d[2]") //intra layer conduction delays
declare("maxsfall",0.001) // max fall-off in prob of connection @ opposite side of column (used by swire)
declare("slambda",colside/3) // spatial length constant for probability of connections, used in swirecut
dvcond[0] = 1   //if presynaptic cell is excit, conduction delay == 1.0 m/s
dvcond[1] = 0.4 //if presynaptic cell is inhib, conduction delay == 0.4 m/s
declare("lambda",30,"hval",0)
declare("mknetnqss",1) // whether COLUMN should make connsnq,cellsnq
declare("DopeRL",0) // whether to use dopamine-RL learning

if(autotune) {
  seadsetting_INTF6 = 2
  wsetting_INTF6 = 1 // use sywv during sim
} else if(useplast) {
  wsetting_INTF6 = 1 // use sywv during sim
  seadsetting_INTF6 = 3 // use plasticity
} else {
  seadsetting_INTF6 = 2
  wsetting_INTF6 = 0
}

if(DopeRL) {
  DOPE_INTF6=1
  EDOPE_INTF6 = IDOPE_INTF6 = 1
  ESTDP_INTF6 = ISTDP_INTF6 = 0  
}

declare("wnqstr","") // intracol wiring NQS (single column)

//* 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) // 2 (takes a lot longer!)
declare("pmatscale",1/scale) // scale for pmat - allows keeping it fixed while changing # of cells in network
declare("wmatscale",1) // scale for wmat - called after setwmat
declare("pmatEE",1,"pmatEI",1,"pmatIE",1,"pmatII",1) // scaling for pmat
declare("onelyr",0) // whether to keep network as a single layer (L2/3)

declare("fc",1)
// declare("EEGain",fc*1,"EIGain",fc*1,"IEGain",fc*1,"IIGain",fc*1)
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)

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("cc12vals",0) // use pmat,wmat,del #s from /u/samn/papers/gif/cc12:936_table1.gif
declare("autogain",0) // automatically set EE,EI,IE,II gains based on whether using cc12vals

if(autogain) {
  print "WARNING: autogain is on!"
  if(cc12vals) {
    EEGain = fc*1
    EIGain = fc*1
    IEGain = fc*1
    IIGain = fc*1
  } else {
    EEGain=4*15/11
    EIGain=15
    IEGain=4*15/11
    IIGain=4*15/11
  }
}

sprint(tstr,"d[%d][%d]",CTYPi,CTYPi)
if(!(i=name_declared("delmscale"))) {
  declare("delmscale",tstr) // scale delm values by this #
}
if(!i) {
  for i=0,CTYPi-1 for j=0,CTYPi-1 delmscale[i][j]=1
}

//* 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(onelyr) {

      cpercol[E2]  = 150 * scale
      cpercol[I2L] =  13 * scale
      cpercol[I2]  =  25 * scale

    } else 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,ii,jj,kk,ct
  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

  if(cc12vals) {

    for case(&ii,E5R,E5B) for case(&jj,E5R,E5B) pmat[ii][jj][0] = 1/11
    pmat[E2][E2][0] = (1/4 + 1/10) / 2
    pmat[E4][E4][0] = 1 / 5.7    
    for case(&ct,E5R,E5B) pmat[E2][ct][0] = 1 / 1.8
    for case(&ct,E5R,E5B) pmat[ct][E2][0] = 1 / 29
    pmat[E4][E2][0] = 1 / 3.6    
    for case(&ct,E5R,E5B) for case(&jj,I5,I5L) pmat[ct][jj][0] = 1 / 10.4
    for case(&ct,I5,I5L) for case(&jj,E5R,E5B) pmat[ct][jj][0] = 1 / 8
    //div by 2 are when cat/rat values differ, so take average
    for case(&ct,I2,I2L) pmat[E2][ct][0] = (1/5 + 1/4) / 2
    for case(&ct,I2,I2L) pmat[ct][E2][0] = (1/6.3 + 1/3.6) / 2
    for case(&ct,I4,I4L) pmat[E4][ct][0] = 1 / 5
    for case(&ct,I4,I4L) pmat[ct][E4][0] = 1 / 10
    for case(&ct,I2,I2L) pmat[E4][ct][0] = 1 / 10    
    for case(&ct,I4,I4L) pmat[E2][ct][0] = (1/12 + 1/5.3) / 2
    for case(&ct,I4,I4L) pmat[ct][E2][0] = (1/2 + 1/3.7) / 2       
    for case(&ct,I2,I2L) for case(&jj,I2,I2L) pmat[ct][jj][0] = (1/4 + 1) / 2
    for case(&ct,I4,I4L) for case(&jj,I4,I4L) pmat[ct][jj][0] = 1/2
    for case(&ct,I5,I5L) for case(&jj,I5,I5L) pmat[ct][jj][0] = 1 / 1.7
    for case(&ct,I4,I4L) for case(&jj,I2,I2L) pmat[ct][jj][0] = 1 / 2 // 1:1 but keeping @ 50%
    // for L6, no data from that paper, so guessing & using KMJ values
    pmat[E6][E5B][0]=0.028 * 3
    pmat[E6][E5R][0]=0.006 * 3
    pmat[E6][E6][0]=0.028 * 3
    pmat[E6][I6L][0]=0.51
    pmat[E6][I6][0]=0.43
    pmat[E6][I6][1]=0.14
    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    

  } else {
    pmat[E2][E2][0]=0.187 
    pmat[E2][E2][1]=0//0.14
    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[,EE,EI,IE,II]) - multiply values in pmat by fctr
proc scalepmat () { local fctr,EE,EI,IE,II,from,to,cl
  fctr=$1 
  for from=0,CTYPi-1 for to=0,CTYPi-1 for cl=0,1 pmat[from][to][cl] *= fctr
  if(numarg()==5) {
    EE=$2 EI=$3 IE=$4 II=$5
    for from=0,CTYPi-1 for to=0,CTYPi-1 for cl=0,1 {
      if(ice(from)) {
        if(ice(to)) {
          pmat[from][to][cl] *= II
        } else {
          pmat[from][to][cl] *= IE
        }
      } else {
        if(ice(to)) {
          pmat[from][to][cl] *= EI
        } else {
          pmat[from][to][cl] *= EE
        }
      }
    }
  }
}

//* 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,ct
  if(cc12vals) {

    for case(&ii,E5R,E5B) for case(&jj,E5R,E5B) delm[ii][jj] = (1.8+1.5)
    for case(&ii,E5R,E5B) for case(&jj,E5R,E5B) deld[ii][jj] = (1.4+0.15)

    delm[E2][E2] = (1.86+2.5)/2 + (1.5+1.7)/2
    deld[E2][E2] = (0.8+1.3)/2 + (0.3+0.79)/2

    delm[E4][E4] = 2 + .6
    deld[E4][E4] = 1.5 // unknown

    for case(&ct,E5R,E5B) delm[E2][ct] = 2 + 1.4
    for case(&ct,E5R,E5B) deld[E2][ct] = 0.9 + 0.3

    for case(&ct,E5R,E5B) delm[ct][E2] = 2.6 //unknown
    for case(&ct,E5R,E5B) deld[ct][E2] = 1.5 //unknown

    delm[E4][E2] = 1.9 + 0.9 // using spiny stellate values, since others close anyway
    deld[E4][E2] = 0.4 + 0.24

    for case(&ct,E5R,E5B) for case(&jj,I5,I5L) delm[ct][jj] = 0.6 + 1.4
    for case(&ct,E5R,E5B) for case(&jj,I5,I5L) deld[ct][jj] = 1 //unknown

    for case(&jj,E5R,E5B) delm[I5][jj] = 4.2 + 1.1 - 2 // -2 since FS faster
    for case(&jj,E5R,E5B) deld[I5][jj] = 2.8 + 0.2 - 1 // -1 since FS faster

    for case(&jj,E5R,E5B) delm[I5L][jj] = 4.2 + 1.1 // LTS slightly slower
    for case(&jj,E5R,E5B) deld[I5L][jj] = 2.8 + 0.2

    for case(&ct,I2,I2L) delm[E2][ct] = (1.2+2.2)/2 + (1.3+0.95)/2
    for case(&ct,I2,I2L) deld[E2][ct] = (.9+1.4)/2 + (.3+.2)/2

    delm[I2][E2] = 3.9 + 1.8 -2 // -2 since FS faster
    deld[I2][E2] = 1.5 + 0.8 -1 // -1 since FS faster

    delm[I2L][E2] = 3.9 + 1.8
    deld[I2L][E2] = 1.5 + 0.8

    for case(&ct,I4,I4L) delm[E4][ct] = 0.7 + 1.2
    for case(&ct,I4,I4L) deld[E4][ct] = 1 //unknown

    delm[I4][E4] = 3.7 + 1.2 -2 // -2 since FS faster
    deld[I4][E4] = 0.3 

    delm[I4L][E4] = 3.7 + 1.2
    deld[I4L][E4] = 0.3

    for case(&ct,I2,I2L) delm[E4][ct] = 0.72 + 0.6
    for case(&ct,I2,I2L) deld[E4][ct] = 0.04 + 0.1

    for case(&ct,I4,I4L) delm[E2][ct] = 0.85 + 1.0
    for case(&ct,I4,I4L) deld[E2][ct] = 0.2 + 0.2

    delm[I4][E2] = (2.8+3+3.5+4.4)/4 + (1.1+1.0+0.8+1.9)/4 -2 // -2 since FS faster
    deld[I4][E2] = (1.6 + 1.1)/2 -0.5 // -0.5 since FS faster

    delm[I4L][E2] = (2.8+3+3.5+4.4)/4 + (1.1+1.0+0.8+1.9)/4 
    deld[I4L][E2] = (1.6 + 1.1)/2 

    for case(&jj,I2,I2L) delm[I2][jj] = (2+2.4+2.3+3.9)/4 + (0.7+1.5+1.0+0.8)/4 -2 // -2 since FS faster
    for case(&jj,I2,I2L) deld[I2][jj] = (1.9+0.8)/2 -0.5 // -0.5 since FS faster

    for case(&jj,I2,I2L) delm[I2L][jj] = (2+2.4+2.3+3.9)/4 + (0.7+1.5+1.0+0.8)/4 
    for case(&jj,I2,I2L) deld[I2L][jj] = (1.9+0.8)/2 

    for case(&jj,I4,I4L) delm[I4][jj] = (2.7+2.9)/2 + (0.8+1.1)/2 -2 // -2 since FS faster
    for case(&jj,I4,I4L) deld[I4][jj] = 1.35 -0.5 // -0.5 since FS faster

    for case(&jj,I4,I4L) delm[I4L][jj] = (2.7+2.9)/2 + (0.8+1.1)/2 
    for case(&jj,I4,I4L) deld[I4L][jj] = 1.35 

    //rest unknown

    for case(&jj,I5,I5L) delm[I5][jj] = (2.7+2.9)/2 + (0.8+1.1)/2 -2 // -2 since FS faster
    for case(&jj,I5,I5L) deld[I5][jj] = 1.35 -0.5 // -0.5 since FS faster

    for case(&jj,I5,I5L) delm[I5L][jj] = (2.7+2.9)/2 + (0.8+1.1)/2  
    for case(&jj,I5,I5L) deld[I5L][jj] = 1.35 


    for case(&jj,I2,I2L) delm[I4][jj] = (2.7+2.9)/2 + (0.8+1.1)/2 -2 // -2 since FS faster
    for case(&jj,I2,I2L) deld[I4][jj] = 1.35 -0.5 // -0.5 since FS faster

    for case(&jj,I2,I2L) delm[I4L][jj] = (2.7+2.9)/2 + (0.8+1.1)/2
    for case(&jj,I2,I2L) deld[I4L][jj] = 1.35 

    // for L6, no data from that paper, so guessing & using KMJ values
    delm[E6][E5B]= 3.8
    deld[E6][E5B]= 1.6

    delm[E6][E5R]= 3.8
    deld[E6][E5R]= 1.6

    delm[E6][E6]= 3.8
    deld[E6][E6]= 1.6

    delm[E6][I6L]= 2.825
    deld[E6][I6L]= 1.4

    delm[E6][I6]= 2.825
    deld[E6][I6]= 1.4

    for case(&ct,E2,E5B,E5R,E6) {
      delm[I6L][ct]= 3.9 + 1.8
      deld[I6L][ct]= 1.5 + 0.8
    }

    for case(&ct,I2,I5,I6L,I6) {
      delm[I6L][ct]= (2+2.4+2.3+3.9)/4 + (0.7+1.5+1.0+0.8)/4 
      deld[I6L][ct]= (1.9+0.8)/2 
    }

    delm[I6][E6]= 3.9 + 1.8 -2 // -2 since FS faster
    deld[I6][E6]= 1.5 + 0.8 -1 // -1 since FS faster

    for case(&ct,I6L,I6) {
      delm[I6][ct]= (2+2.4+2.3+3.9)/4 + (0.7+1.5+1.0+0.8)/4 -2 // -2 since FS faster
      deld[I6][ct]= (1.9+0.8)/2 -0.5 // -0.5 since FS faster
    }

  } else {
    for from=0,CTYPi-1 for to=0,CTYPi-1 {
      if(synloc[from][to]==DEND) {
        delm[from][to]=4 * delmscale[from][to]
        deld[from][to]=1
      } else {
        delm[from][to]=2.0 * delmscale[from][to]
        deld[from][to]=0.2
      }
    }
  }
}

//* weight params
//** delay all 2+/-0.02 within column for now
proc setwmat () { local ii,jj,from,to,sy,gn,c,fc,ct
  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

  if(cc12vals) {

    for case(&ii,E5R,E5B) for case(&jj,E5R,E5B) wmat[ii][jj][AM2][0] = 1.7
    wmat[E2][E2][AM2][0] = (1.7 + 1.4) / 2
    wmat[E4][E4][AM2][0] = 1.1
    for case(&ct,E5R,E5B) wmat[E2][ct][AM2][0] = 1.4
    for case(&ct,E5R,E5B) wmat[ct][E2][AM2][0] = 0.3
    wmat[E4][E2][AM2][0] = (1.1+3.3+5.9)/3
    for case(&ct,E5R,E5B) for case(&jj,I5,I5L) wmat[ct][jj][AM2][0] = 0.9
    for case(&ct,E5R,E5B) wmat[I5][ct][GA][0] = 1.23
    for case(&ct,E5R,E5B) wmat[I5L][ct][GA2][0] = 1.23
    //div by 2 are when cat/rat values differ, so take average 
    for case(&ct,I2,I2L) wmat[E2][ct][AM2][0] = (1.9+3.1)/2
    wmat[I2][E2][GA][0] = 0.65
    wmat[I2L][E2][GA2][0] = 0.65

    for case(&ct,I4,I4L) wmat[E4][ct][AM2][0] = 3.7
    wmat[I4][E4][GA][0] = 0.85
    wmat[I4L][E4][GA2][0] = 0.85

    for case(&ct,I2,I2L) wmat[E4][ct][AM2][0] = 1.2

    for case(&ct,I4,I4L) wmat[E2][ct][AM2][0] = 1.3

    wmat[I4][E2][GA][0] = (1.75 + (.43+1)/2)/2
    wmat[I4L][E2][GA2][0] = (1.75 + (.43+1)/2)/2

    for case(&ct,I2,I2L) wmat[I2][ct][GA][0] = (2 + .7 + .8 + 1.3)/4
    for case(&ct,I2,I2L) wmat[I2L][ct][GA2][0] = (2 + .7 + .8 + 1.3)/4

    for case(&ct,I4,I4L) wmat[I4][ct][GA][0] = (.4+2.7)/2
    for case(&ct,I4,I4L) wmat[I4L][ct][GA2][0] = (.4+2.7)/2

    for case(&ct,I5,I5L) wmat[I5][ct][GA][0] = (.4+2.7)/2
    for case(&ct,I5,I5L) wmat[I5L][ct][GA2][0] = (.4+2.7)/2

    for case(&ct,I2,I2L) wmat[I4][ct][GA][0] = (.4+2.7)/2
    for case(&ct,I2,I2L) wmat[I4L][ct][GA2][0] = (.4+2.7)/2

    // for L6, no data from that paper, so guessing & using KMJ values -- appear similar range here
    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[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

  } else {
    fc = 1
    //*** neocx -> neocx
    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 * fc
    wmat[E2][E5R][AM2][0]=0.93 * fc
    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 * fc
    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 * fc
    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 
  }
}

//* scalewmat(fctr) - multiply values in wmat by fctr
proc scalewmat () { local fctr,from,to,sy,c
  fctr=$1
  for from=0,CTYPi-1 for to=0,CTYPi-1 {      
    for sy=0,STYPi-1 for c=0,colr-1 wmat[from][to][sy][c] *= fctr
  }
}

// %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

//** 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,1) // 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,1) // 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,1) // 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
    }    
  }
  ce.o(0).finishdvir()
  {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 localobj nq
  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,mknetnqss,1))
    col[id]=gcol[y][x]
    col[id].verbose=verbose_INTF6
    //if(strlen(wnqstr)>0) { // strlen is an undefined function in this particular model
    //  nq=new NQS(wnqstr)
    //  col[id].wirenq(nq)
    //  nqsdel(nq)
    //} else if(swire>0) {
    if(swire>0) {
      col[id].setcellpos(vlayerz,layerzvar,pseed*(id+1),colside)
      if(swire==1) {
        col[id].swire(col[id].wseed,maxsfall)
      } else if(swire==3) {
        swirecutfl(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda) // uses intf6::floc for speedup
      }
    } else col[id].wire(col[id].wseed)
    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
  from=$1 to=$2 sy=$3 w=$4
  for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) {
    if(xo.pre.type==from && xo.syn.type==to) xo.weight(sy)=w
  }
}

//* conprob
func conprob () { // calculate connection probability using global slambda and hval
  if ($1<=slambda) return hval else return hval*exp(1-$1/slambda) // probability of connect
}

//** swirecutfl (col,EEsq,EIsq,IEsq,IIsq[,seed,slambda]) - spatial wiring: wires the column using pmat and cell positions
//                   (wiring probability effected by distance btwn cells)
// seed is random # seed
// slambda specifies length-constant for spatially-dependent fall-off in wiring probability
// at one slambda away, probability of connect is value in pmat
proc swirecutfl () { local x,y,z,ii,jj,a,del,prid,poid,prty,poty,dv,lseed,h,prob,slambda,dsq,dist,slambdasq,\
                  EEsq,EIsq,IEsq,IIsq,ic1,ic2,pdx,n,nn,maxd,rsz\
               localobj col,ce,vidx,vdel,vdist,vwt1,vwt2,vtmp,vpre,opr,opo,st,rdm,vprob,oq,ol,mmaxd
  oq=new NQS("id1","id2","del","dist","wt1","wt2","dist0","prob") oq.verbose=0 ol=new List() mmaxd=new Matrix(2,2)
  a=allocvecs(vidx,vdel,vtmp,vdist,vwt1,vwt2,vprob,vpre) z=0
  {ol.append(vpre) ol.append(vidx) ol.append(vdel) ol.append(vdist) ol.append(vwt1) ol.append(vwt2)}
  col=$o1 ce=col.ce EEsq=$2 EIsq=$3 IEsq=$4 IIsq=$5
  if (argtype(6)==0) lseed=$6 else lseed=1234  
  if(numarg()>6) slambda=$7 else slambda=10
  if(slambda<=0){
    printf("swirecut WARN: invalid slambda=%g, setting slambda to %g\n",colside/3)
    slambda=colside/3
  }
  slambdasq = slambda^2 // using squared distance
  vrsz(1e4,vidx,vdel,vdist,vtmp,vpre)
  rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator
  rdm.uniform(0,1) rsz=ce.count^2 vprob.resize(rsz) vprob.setrand(rdm) pdx=0
  {mmaxd.x(0,0)=sqrt(EEsq) mmaxd.x(0,1)=sqrt(EIsq) mmaxd.x(1,0)=sqrt(IEsq) mmaxd.x(1,1)=sqrt(IIsq)}
  // Create the connectivity NQS table.
  if(col.mknetnqss) {nqsdel(col.connsnq) col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")}
  for y=0,ce.count-1 { opr=ce.o(y)
    vrsz(0,vidx,vdel,vdist,vwt1,vwt2)
    prid=opr.id prty=opr.type ic1=ice(opr.type)
    for poty=0,CTYPi-1 if (col.numc[poty]!=0 && (h=hval=col.pmat[prty][poty])>0) {
      ic2 = ice(poty)
      oq.clear()
      if ((n=opr.floc(opr.xloc,opr.yloc,1e9,oq.getcol("id2"),oq.getcol("dist"),mmaxd.x(ic1,ic2),poty))==0) {
        if(verbose) printf("swirecutfl: No possibilities for connections found from %d\n",opr.id)
      } else {
        oq.getcol("dist0").copy(oq.getcol("dist"))
        oq.pad(n)
        oq.getcol("dist0").apply("conprob") // turn distances into probabilities
        if (pdx+n>=rsz) { rdm.uniform(0,1) vprob.setrand(rdm) pdx=0 }
        oq.getcol("prob").copy(vprob,pdx,pdx+n-1)  // pick out some random probabilities to compare to
        pdx+=n 
        oq.calc("<dist0>.sub(<prob>)") 
        if ((nn=oq.select("dist0",">",0,"id2","!=",prid))==0) { // looks for dist_prob > prob
          if(verbose) printf("swirecutfl: No connections found from %d out of %d assayed\n",opr.id,n)
        } else {
          oq.cpout()
          //oq.fill("id1",prid)
          rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty])
          oq.getcol("del").setrand(rdm)
          if (col.syty1[prty][poty]>=0) oq.getcol("wt1").fill(col.wmat[prty][poty][col.syty1[prty][poty]]) else oq.getcol("wt1").fill(0)
          if (col.syty2[prty][poty]>=0) oq.getcol("wt2").fill(col.wmat[prty][poty][col.syty2[prty][poty]]) else oq.getcol("wt2").fill(0)
          {vwt1.append(oq.getcol("wt1")) vwt2.append(oq.getcol("wt2"))}
          {vidx.append(oq.getcol("id2")) vdel.append(oq.getcol("del"))}
          vdist.resize(vdist.size+nn)
          if(synloc[prty][poty]==DEND) vdist.fill(1,vdist.size-nn,vdist.size-1) else vdist.fill(0,vdist.size-nn,vdist.size-1)
        }
      }
    }
    if(vidx.size>0) {
      if(wsetting_INTF6==1) opr.setdvi(vidx,vdel,vdist,1,vwt1,vwt2) else opr.setdvi(vidx,vdel,vdist,1)
      if(col.mknetnqss){vpre.resize(vidx.size) vpre.fill(prid) col.connsnq.append(ol)}
//        for ii=0,vidx.size-1 col.connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
    }
  }
  col.ce.o(0).finishdvir()
  dealloc(a)
  nqsdel(oq)
  if(verbose) printf("\n")
}

//* randsywts(col) - randomize synaptic weights
proc randsywts () { local w,i,j,k,sy,sead,a localobj xo,vidx,vw1,vw2,col,ce,rdm
  col=$o1 ce=col.ce
  if(numarg()>1) sead=$2 else sead=1234
  rdm=new Random() rdm.ACG(sead)
  a=allocvecs(vidx,vw1,vw2)
  for ltr(xo,ce) {
    {vrsz(xo.getdvi(vidx),vw1,vw2) vw1.fill(0) vw2.fill(0)}
    for vtr(&i,vidx,&j) {
      if(synloc[xo.type][ce.o(i).type]==DEND) {
        if(ice(xo.type)) sy=GA2 else sy=AM2
      } else {
        if(ice(xo.type)) sy=GA else sy=AM
      }
      w = wmat[xo.type][ce.o(i).type][sy+0][0]
      if(randsy > 0) {
        vw1.x(j)=MAXxy(0,rdm.uniform((1-randsy)*w,(1+randsy)*w))
      } else vw1.x(j)=MAXxy(0,rdm.normal(w,w/(-randsy)))
      if((w=wmat[xo.type][ce.o(i).type][sy+1][0])>0)vw2.x(j)=vw1.x(j)*NMAMR// only NM2 for wt2
    }
    xo.setsywv(vw1,vw2)
  }
  dealloc(a)
}

//* drcelllocs - draw location of cells
proc drcelllocs () { local i,gvt,ty
  gvt=gvmarkflag
  {gvmarkflag=1 col[CDX].cellsnq.marksym="O"}
  if(numarg()>0) {
    if(col[CDX].cellsnq.select("ty",$1)) col[CDX].cellsnq.gr("yloc","xloc",0,1,3)
  } else for i=0,1 if(col[CDX].cellsnq.select("ice",i)) col[CDX].cellsnq.gr("yloc","xloc",0,i+2,3)
  gvmarkflag=gvt
}

//* drcelldiv(cellid) - draw outputs from a cell
proc drcelldiv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt localobj cnq,nq
  gvt=gvmarkflag gvmarkflag=1
  cnq=col[CDX].connsnq
  nq=col[CDX].cellsnq
  id1=$1 nq.verbose=0
  nq.select(-1,"id",id1)
  x1=nq.v[5].x(nq.ind.x(0))
  y1=nq.v[6].x(nq.ind.x(0))
  ic1=nq.v[4].x(nq.ind.x(0))
  if(ic1) clr1=3 else clr1=2
  nq.select("id",id1)
  nq.gr("yloc","xloc",0,clr1,15)
  cnq.select(-1,"id1",id1)
  cnq.verbose=0
  for vtr(&i,cnq.ind) {
    id2 = cnq.v[1].x(i)
    if(nq.select("id",id2)) {
      x2=nq.fetch("xloc")
      y2=nq.fetch("yloc")
      ic2=nq.fetch("ice")
      if(ic2) clr2=3 else clr2=2
      nq.gr("yloc","xloc",0,clr2,8)    
      drline(x1,y1,x2,y2,g,clr1,7)
    }
  }
  cnq.verbose=nq.verbose=1
  gvmarkflag=gvt
}

//* drcellconv(cellid) - draw inputs to a cell
proc drcellconv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt localobj cnq,nq
  gvt=gvmarkflag gvmarkflag=1
  cnq=col[CDX].connsnq
  nq=col[CDX].cellsnq
  id2=$1 nq.verbose=0
  nq.select(-1,"id",id2)
  x2=nq.v[5].x(nq.ind.x(0))
  y2=nq.v[6].x(nq.ind.x(0))
  ic2=nq.v[4].x(nq.ind.x(0))
  if(ic2) clr2=3 else clr2=2
  nq.select("id",id2)
  nq.gr("yloc","xloc",0,clr2,15)
  cnq.select(-1,"id2",id2) 
  cnq.verbose=0
  for vtr(&i,cnq.ind) {
    id1 = cnq.v[0].x(i)
    if(nq.select("id",id1)) {
      x1=nq.fetch("xloc")
      y1=nq.fetch("yloc")
      ic1=nq.fetch("ice")
      if(ic1) clr1=3 else clr1=2
      nq.gr("yloc","xloc",0,clr1,8)    
      drline(x1,y1,x2,y2,g,clr1,7)
    }
  }
  cnq.verbose=nq.verbose=1
  gvmarkflag=gvt
}

//* mkwgnq(col) - return nqs with weight info
obfunc mkwgnq () { local a,i,j,id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist\
                localobj mycel,vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6,wgnq,col
  a=allocvecs(vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6) col=$o1
  {nqsdel(wgnq) wgnq=new NQS("id1","id2","wg","ice","ty1","ty2","wt1","wt2","del","dist")}
  for i=0,col.allcells-1 {
    mycel=col.ce.o(i)
    mycel.getdvi(myv,myv1,myv2,myv3,myv4,myv5,myv6)
    id1=i
    ic=ice(mycel.type)
    ty1=mycel.type
    {vrsz(mycel.getdvi,vwt1,vwt2) vwt1.fill(0) vwt2.fill(0)}
    if(wsetting_INTF6==1) mycel.getsywv(vwt1,vwt2)
    for j=0,myv6.size-1 {
      id2=myv.x(j)
      ty2=col.ce.o(id2).type
      wg = myv6.x(j)
      wt1 = wg * vwt1.x(j)
      if(ic) wt2 = vwt2.x(j) else wt2 = wg * vwt2.x(j)
      del = myv1.x(j)
      dist = myv5.x(j)
      wgnq.append(id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist)
    }
  }
  dealloc(a)
  return wgnq
}

//* 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()
scalepmat(pmatscale,pmatEE,pmatEI,pmatIE,pmatII)

//** setup synapse locations,delays,wmat
setsynloc()
setdelmats()
setwmat() // new KMJ version
if(wmatscale!=1) scalewmat(wmatscale)

scrsz=50*1e3
double scr[scrsz]

//** make cells, columns, wire columns
mkcolnqs()
mkcols()
wirecols(1)
if(randsy) randsywts(col,randsy)