Parallel network simulations with NEURON (Migliore et al 2006)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:64229
The NEURON simulation environment has been extended to support parallel network simulations. The performance of three published network models with very different spike patterns exhibits superlinear speedup on Beowulf clusters.
Reference:
1 . Migliore M, Cannia C, Lytton WW, Markram H, Hines ML (2006) Parallel network simulations with NEURON. J Comput Neurosci 21:119-29 [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:
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Methods;
Implementer(s): Hines, Michael [Michael.Hines at Yale.edu];
/
netmod
parscalebush
AMPA.mod *
arhRT03.mod *
cadecay.mod
cadRT03.mod
cah.mod
calRT03.mod
catRT03.mod *
GABAa.mod *
GABAb.mod *
intf.mod
k2RT03.mod *
kahpRT03.mod
kaRT03.mod *
kca.mod *
kcRT03.mod
kdr.mod
kdrp.mod
kdrRT03.mod *
kmRT03.mod *
misc.mod
myfit.mod
na.mod
nafRT03.mod *
nap.mod
napRT03.mod *
NMDA.mod *
nstim.mod
stats.mod
vecst.mod
batch.hoc
bg
bg_cvode.inc
boltz_cvode.inc
geneval_cvode.inc
geom.hoc
init.hoc
netcon.inc
network.hoc
nqsnet.hoc
nspike.dat
params.hoc
parnetwork.hoc
parnqsnet.hoc
perfrun.hoc
prebatch_.hoc
run.hoc
spkplt.hoc *
x_vs_nspike.hoc
                            
// $Id: parnqsnet.hoc,v 1.2 2006/02/11 13:34:17 hines Exp $

// load_file("nqsnet.hoc")
objref sp[5],nq[5]
datestr="06feb04"
scale=1 // redefined in network.hoc
double pmat[CTYPi][CTYPi][STYPi],numc[CTYPi],ix[CTYPi],wmat[CTYPi][CTYPi][STYPi]
obfunc cobj(){} // return cell given gid (defined in network.hoc)

//* routines
//** styp() sets synapse type based on presynaptic cell
func styp () { local pr,po
  pr=$1 po=$2
  if (pr==IN && po==IN) { return GA 
  } else if (pr==IN) { return IH
  } else if (pr==SU) { return EX
  } else if (pr==SM) { return AM
  } else printf("styp ERR %s->%s not classified",CTYP.object(pr).s,CTYP.object(po).s)
}

for ii=0,2 {sp[ii]=new NQS("PR","PO","SID","WT","DEL") sp[ii].clear(1e5*scale)}
batch_flag=1
for ii=0,2 sp[ii].mo(1)
batch_flag=0
// mkspmat() uses sp and sp[1], returns sp
proc mkspmat () { local a,sn,ii,jj,n,x,y,z,z1 localobj vr,vo
  x=$1 y=$2 z=$3
  rdm.ACG(903342)
  a=allocvecs(vr,vo)
  sn=pmat[x][y][z]*numc[x]*numc[y] // number of synapses
  for ii=0,1 sp[ii].clear
  sp.mo(2)
  sp.pad(1.5*sn) // leave room
  rdm.discunif(ix[x],ix[x]+numc[x]-1) // randomly chosen presynaptic cells
  PRv.setrand(rdm)
  rdm.discunif(ix[y],ix[y]+numc[y]-1)  // randomly chosen postsynaptic cells
  POv.setrand(rdm)
  if (z==EX) z1=AM else z1=z  
  sytyp(c[y],z1,vr) // will want to replace this since got rid of c[]
  SIDv.vfill(vr)
  rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2) // calc random weights
  WTv.setrand(rdm)   // set the weight column
  rdm.normal(2.2,0.1) // random delays
  DELv.setrand(rdm)    
  sp.elimrepeats("PR","PO","SID") // so don't have 2 connections between same cells
  sp.shuffle()  // randomize again
  sp.pad(sn) // reduce the length down to desired
  if (z==EX) { // now have to repeat everything, also will need for IX
    z1=NM
    sp[1].mo(2) // reset the PRv,POv,... pointers
    sp[1].cp(sp) // this will be the same as prior
    rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2)
    WTv.setrand(rdm)
    sytyp(c[y],z1,vo) // vo and vr represent list indices for cell synlists
    vo.sub(vr) // difference between the NMDA synnums and AMPA synnums
    if (! vo.ismono(0)){ printf("mkspmat ERRA\n") return }
    SIDv.add(vo.x[0])
    sp.grow(sp[1])
  }
  sp.mo(2) // reset the PRv,POv,... pointers
  WTv.negwrap() // take any neg weights and make them positive
  dealloc(a)
}

// mkpomat() only make the array for a single postsynaptic cell
proc mkpomat () { local a,sn,ii,jj,n,x,y,yid,z,z1 localobj vr,vo
  x=$1 y=$2 z=$3 yid=$4
  // rdm.ACG(903342)
  // We want to have a postsynaptic gid specific stream but streams from
  // different gids should be statistically independent.
  // A stream is random but we want to be able to reproduce it at will
  // Assume this is the only place where this is required
  mcell_ran4_init(yid) // for independent sims add offset > i*allcells
  rdm.MCellRan4(1) //avoid 0 want to always start at known seed
  a=allocvecs(vr,vo)
  sn=pmat[x][y][z]*numc[x]*1 // number of synapses
  if (sn<1) { // treat sn<1 as a prob of connecting rather than a density
    if (rdm.uniform(0,1)<sn) sn=1 else sn=0 // flip coin
  } else if (sn<5) sn=round(sn) // possibly round it up
  for ii=0,1 sp[ii].clear
  if (sn==0) { dealloc(a) return }
  sp.mo(2)
  sp.pad(1.5*sn) // leave room
  rdm.discunif(ix[x],ix[x]+numc[x]-1) // randomly chosen presynaptic cells
  PRv.setrand(rdm)
  POv.fill(yid)
  if (z==EX) z1=AM else z1=z  
  sytyp(cobj(yid),z1,vr) // brief syncode.hoc routine translates synaptic types
  SIDv.vfill(vr)
  rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2) // calc random weights
  WTv.setrand(rdm)   // set the weight column
  rdm.normal(2.2,0.1) // random delays
  DELv.setrand(rdm)    
  sp.elimrepeats("PR","PO") // so don't have 2 connections between same cells
  sp.shuffle()  // randomize again
  sp.pad(sn) // reduce the length down to desired
  if (z==EX) { // now have to repeat everything, also will need for IX
    z1=NM
    sp[1].mo(2) // reset the PRv,POv,... pointers
    sp[1].cp(sp) // this will be the same as prior
    rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2)
    WTv.setrand(rdm)
    sytyp(cobj(yid),z1,vo) // vo and vr represent list indices for cell synlists
    vo.sub(vr) // difference between the NMDA synnums and AMPA synnums
    if (! vo.ismono(0)){ printf("mkspmat ERRA\n") return }
    SIDv.add(vo.x[0])
    sp.grow(sp[1])
  }
  sp.mo(2) // reset the PRv,POv,... pointers
  WTv.negwrap() // take any neg weights and make them positive
  dealloc(a)
}

proc mkall () {  local x,y,z localobj fname
  fname=new String()
  for case(&y,FP,TP,B5) { // FP,TP,B5
    sp[2].clear
    for case(&x,SM,FP,TP,B5) for case(&z,AM,GA,GB,EX,E2) {
      if (pmat[x][y][z]!=0) {
        mkspmat(x,y,z) 
        sp[2].grow(sp)
      }
    }
    sprint(fname.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
    sp[2].sv(fname.s)
  }
}

// nqs2txt()
proc nqs2txt () { local ii,last,lnum,n,fl1 localobj f,fn
  nq=new NQS("GID","LINE","NUML","TELL")
  fn=new String()
  fl1=lnum=n=0
  f=new File()
  sprint(tstr,"data/%sS%02d.syns",datestr,scale)
  f.wopen(tstr)
  for case(&y,FP) { // FP,TP,B5 need proper order
    sprint(fn.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
    sp.rd(fn.s)
    sp.sort("SID") sp.sort("PR") sp.sort("PO")
    if (fl1==0) {  
      last=sp.v[1].x[0]
      nq.append(last,lnum+1,-1,f.tell)
      fl1=1
    }
    for ii=0,sp.size(1)-1 {
      PR=sp.v[0].x[ii] PO=sp.v[1].x[ii] SID=sp.v[2].x[ii] WT=sp.v[3].x[ii] DEL=sp.v[4].x[ii]
      if (PO==last) {
        lnum+=1 n+=1
        f.printf("%d %d %d %g %g\n",PO,PR,SID,WT,DEL)
      } else {
        nq.set("NUML",-1,n)
        n=0
        if (PO==last+1) {
          nq.append(PO,lnum+1,-1,f.tell+1)
        } else { 
          printf("%d is not a target\n")
          nq.append(PO,lnum+1,0,f.tell+1)
        }
        last=PO
        ii-=1 // do this one again
      } 
    }
    nq.set("NUML",-1,n)
  }
  sprint(tstr,"data/%sS%02d.index",datestr,scale)
  f.wopen(tstr)
  for ii=0,nq.size(1)-1 {
    f.printf("%d %d %d %d\n",nq.v[0].x[ii],nq.v[1].x[ii],nq.v[2].x[ii],nq.v[3].x[ii])
  }
  f.close
}

// nqssplit() // divide into multiple smaller files
proc nqssplit () { local ii,last,lnum,n,fl1 localobj f,fn
  fn=new String()
  fl1=lnum=n=0
  f=new File()
  for case(&y,FP,TP,B5) {
    sprint(fn.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
    sp.rd(fn.s)
    sp.sort("PO")
    for (ii=ix[y];ii<ix[y]+numc[y];ii+=200) {
      sp.select("PO","[]",ii,ii+200-1)
      sprint(tstr,"data/S%dzip/%s_%d-%d_S%02d.nqs",scale,datestr,ii,ii+200-1,scale)
      sp.out.sv(tstr)
    }
  }
}

proc ncl2mat () { local ii,wt,del,pc,pr,po,prid,poid,ps,sy,syn localobj Xo
  if (!isassigned(sp)) {
    sp=new NQS("PRID","POID","PR","PO","SYNID","DEL","WT0","THRESH") }
  for ii=0,ncl.count-1 { 
    Xo=ncl.o(ii)
    pc=id(Xo.pre) prid=id(Xo.precell) sy=id(Xo.syn) poid=id(Xo.postcell)
    pr=ojtnum(Xo.precell) po=objnum(Xo.postcell)
    if (!isojt(Xo.postcell,nil)) syn=Xo.postcell.synlist.index(Xo.syn)
    if (poid==-1 && sy==-1) {  // nil targ
    } else if (poid==-1) { // ACELL targ -- not implemented
    } else if (prid==-1) { // ACELL src
      sp.append(pc,poid,ojtnum(Xo.pre),po,syn,Xo.delay,Xo.weight,Xo.threshold)
    } else {
      sp.append(prid,poid,pr,          po,syn,Xo.delay,Xo.weight,Xo.threshold)
    }
  }
}

//* fill everything in and save it
// mksfc()     // synfire chain
// ncl2mat()
// mkall()

// go to global id's -- can leave fp alone since starts at 0
proc rel2glbl () {
  sp.select("PRID",TP) sp.spr("<PR>.add(idtp)") sp.delect()
  sp.select("POID",TP) sp.spr("<PO>.add(idtp)") sp.delect()
  sp.select("PRID",B5) sp.spr("<PR>.add(idb5)") sp.delect()
  sp.select("POID",B5) sp.spr("<PO>.add(idb5)") sp.delect()
  sp.select("PRID",SM) sp.spr("<PR>.add(idstim)") sp.delect()
}