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]
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
parbulbNet
README *
cadecay.mod *
flushf.mod *
kA.mod *
kca.mod *
kfasttab.mod *
kM.mod *
kslowtab.mod *
lcafixed.mod *
nafast.mod *
nagran.mod *
nmdanet.mod *
bulb.hoc
calcisilag.hoc *
ddi_baseline.gnu *
ddi_baseline.ses *
experiment_ddi_baseline.hoc *
experiment_odour_baseline.hoc *
granule.tem *
init.hoc *
input.hoc *
input1 *
mathslib.hoc *
mitral.tem *
modstat
mosinit.hoc *
odour_baseline.gnu *
odour_baseline.ses *
par_batch1.hoc
par_bulb.hoc
par_calcisilag.hoc
par_experiment_ddi_baseline.hoc
par_granule.tem
par_init.hoc
par_input.hoc
par_mitral.tem
par_netpar.hoc
par_notes
parameters_ddi_baseline.hoc *
parameters_odour_baseline.hoc *
screenshot.png *
tabchannels.dat *
tabchannels.hoc *
test1.sh
                            
// calcisilag.hoc
// Olfactory bulb network model: calculate inter-spike interval
//                             : and lag time statistics
// Andrew Davison, The Babraham Institute, 2000.

/* The following procedures are defined for writing results to file
 * fileroot is the filename root. A suffix will be added to
 * this, e.g. fileroot.synch for print_si().
 * i,j are the mitral cell indices
 *
 * print_smooth_hist(variance,fileroot)
 * print_gran_smooth_hist(variance,fileroot)
 * print_spiketimes(i,j,fileroot)
 * print_raster(fileroot)
 * print_gran_raster(fileroot)
 * print_isis(i,j,fileroot)
 * print_isi_stats(fileroot)
 * print_lags(i,j,fileroot) 	
 * print_si(fileroot)
 */

// Variables used in this file
strdef tstr
objref work, work2, outputarray, isi, lags, hist, tmpspike
work 	= new Vector()
work2 	= new Vector()
outputarray = new Matrix(nmitx,nmity)
isi 	= new Vector()
lags 	= new Vector()
tmpspike = new Vector()

// Procedures for processing spike times -------------------------------

proc calc_isis() { local i,j,k,n // 3 args - indices of mitral cell, transient time
  if ($1 > nmitx || $2 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i = int($1)
  j = int($2)
  isi.resize(0)
  n = mit[i][j].spiketimes.size()
  if (n > 1) {
    for k = 1,n-1 {
      if (mit[i][j].spiketimes.x[k-1] > $3) {
        isi.append(mit[i][j].spiketimes.x[k]-mit[i][j].spiketimes.x[k-1])
      }
    }
  }
}

proc calc_gran_isis() { local i,j,k,n // 3 args - indices of granule cell, transient time
  if ($1 > ngranx || $2 > ngrany) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i = int($1)
  j = int($2)
  isi.resize(0)
  n = gran[i][j].spiketimes.size()
  if (n > 1) {
    for k = 1,n-1 {
      if (gran[i][j].spiketimes.x[k-1] > $3) {
        isi.append(gran[i][j].spiketimes.x[k]-gran[i][j].spiketimes.x[k-1])
      }
    }
  }
}

func minisi() { local i,j,min	// find shortest mean ISI
  min = 1e6
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      if (isi.size() > 0) {
        if (isi.mean() < min) { min = isi.mean() }
      }
    }
  }
  return min
}
      
  
proc rate_array() { local i,j
  for i = 0,nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      if (isi.size() > 0) {
        outputarray.x[i][j] = 1000/isi.mean()
      } else {
        outputarray.x[i][j] = 0
      }
    }
  }  
  print "Max: ",arraymax(outputarray)
  outputarray.muls(1/arraymax(outputarray))
} 

proc calc_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time
    if ($1 > nmitx || $2 > nmity || $3 > nmitx || $4 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (mit[i2][j2].spiketimes.size > 0) {
    for k = 1,mit[i1][j1].spiketimes.size()-2 {
      if (mit[i1][j1].spiketimes.x[k] > $5) {
        work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k]
        isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k]
        if (min > isiprev/2 && min < isinext/2) {
          lags.append(min)
        }
      }
    }
  }
}

proc calc_phase_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time
    if ($1 > nmitx || $2 > nmity || $3 > nmitx || $4 > nmity) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (mit[i2][j2].spiketimes.size > 0) {
    for k = 1,mit[i1][j1].spiketimes.size()-2 {
      if (mit[i1][j1].spiketimes.x[k] > $5) {
        work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k]
        isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k]
        if (min > isiprev/2 && min < isinext/2) {
          if (min < 0) {
            lags.append(min/isiprev)
          } else {
            lags.append(min/isinext)
          }
        }
      }
    }
  }
}

proc calc_gran_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of granule cells, transient time
    if ($1 > ngranx || $2 > ngrany || $3 > ngranx || $4 > ngrany) {
    print "Sorry - index out of range. Please try again."
    return
  }
  i1 = int($1)
  j1 = int($2)
  i2 = int($3)
  j2 = int($4)
  lags.resize(0)
  // for each spiketime in cell 1, find closest spike in cell 2
  // Note: first and last spikes ignored since can't calculate previous ISI
  if (gran[i2][j2].spiketimes.size > 0) {
    for k = 1,gran[i1][j1].spiketimes.size()-2 {
      if (gran[i1][j1].spiketimes.x[k] > $5) {
        work = gran[i2][j2].spiketimes.c.add(-gran[i1][j1].spiketimes.x[k])
        minidx = work.c.abs.min_ind()
        min = work.x[minidx]
        isiprev = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k-1]
        isinext = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k+1]
        if (min < isiprev/2 && min > isinext/2) {
          lags.append(min)
        }
      }
    }
  }
}

proc time_hist() { // 1 arg - time step
  work.resize(0)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      work.append(mit[i][j].spiketimes)
    }
  }
  hist = work.histogram(0,tstop,$1)
  hist.printf("%d\n")
}

func synch_index() { local i1,j1,i2,j2,n
  synchindex = 0
  n = 0
  for i1 = 0, nmitx-1 {
    for j1 = 0, nmity-1 {
      if (mit[i1][j1].spiketimes.size() > 0) {
        for i2 = 0, nmitx-1 {
          for j2 = 0, nmity-1 {
            if (i1 != i2 || j1 != j2) {
              calc_phase_lags(i1,j1,i2,j2,ttrans)
              n += lags.size()
              synchindex += lags.reduce("abs",0)
            }
          }
        }
      }
    }
  }
  if (n > 0) {
    return synchindex/n
  } else {
    return 1e6
  }
}


func phaselock_index() { local n,i1,j1,i2,j2
  synchindex = 0
  n = 0
  for i1 = 0, nmitx-1 {
    for j1 = 0, nmity-1 {
      if (mit[i1][j1].spiketimes.size() > 0) {
        for i2 = 0, nmitx-1 {
          for j2 = 0, nmity-1 {
            if (i1 != i2 || j1 != j2) {
              calc_phase_lags(i1,j1,i2,j2,ttrans)
              if (lags.size() > 1) {
                synchindex += lags.var()
                n += 1
              }
            }
          }
        }
      }
    }
  }
  if (n > 0) {
    synchindex = sqrt(synchindex/n)
    return synchindex
  } else {
    return 1e6
  }
}

// Procedures for writing out data --------------------------------

proc calc_smooth_hist() { // 1 arg - variance
  work.resize(0)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      work.append(mit[i][j].spiketimes)
    }
  }
  hist = work.sumgauss(0,tstop,1,$1)
}

proc print_smooth_hist() { // 2 args - variance, filename root
  calc_smooth_hist($1)
  sprint(filename,"%s.smhist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Mitral cell smoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
/*
  work.resize(0)
  hist.remove(0,ttrans)
  work.spctrm(hist)
  sprint(filename,"%s.pow",$s2)
  outfile.wopen(filename)
  outfile.printf("# Power spectrum of Mitral cell smoothed histogram\n")
  work.printf(outfile,"%9.5f\n")
  outfile.close()
  */
}

proc calc_gran_smooth_hist() { // 1 arg - variance
  work.resize(0)
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
      work.append(gran[i][j].spiketimes)
    }
  }
  hist = work.sumgauss(0,tstop,1,$1)
}

proc print_gran_smooth_hist() { // 2 args - variance, filename root
  calc_gran_smooth_hist($1)
  sprint(filename,"%s.gran.smhist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Granule cell smoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
  //work.resize(0)
  //hist.remove(0,ttrans)
  //work.spctrm(hist)
  //sprint(filename,"%s.gran.pow",$s2)
  //outfile.wopen(filename)
  //outfile.printf("# Power spectrum of Granule cell smoothed histogram\n")
  //work.printf(outfile,"%9.5f\n")
  //outfile.close()
}

proc print_gran_hist() { // 2 args - binsize, filename root
  work.resize(0)
  pnm.pc.context("post_raster(ngranx, ngrany, grangid)\n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
	take_raster(i, j, grangid)
      work.append(tmpspike)
    }
  }
  hist = work.histogram(0,tstop,$1)
  sprint(filename,"%s.gran.hist",$s2)
  outfile.wopen(filename)
  outfile.printf("# Granule cell unsmoothed histogram\n")
  hist.printf(outfile,"%8.3f\n")
  outfile.close()
}

proc print_spiketimes() { // 3 args - indices of mitral cell plus filename root
  if (numarg() == 3) {
    sprint(filename,"%s_%d_%d.isi",$s3,$1,$2)
    outfile.wopen(filename)
    outfile.printf("# Spiketimes for mitral cell [%d][%d]",$1,$2)
    mit[$1][$2].spiketimes.printf(outfile,"%10.3f")
  }
  if (numarg() == 1) {
    sprint(filename,"%s.isi",$s1)
    outfile.wopen(filename)
    for i = 0, nmitx-1 {
      for j = 0, nmity-1 {
        outfile.printf("[%d][%d]",i,j)
        mit[i][j].spiketimes.printf(outfile,"%10.3f")
      }
    }
  }
  outfile.close()
}

proc post_raster() {local i, j
  for i = 0, $1-1 {
    for j = 0, $2-1 {
	if (pnm.pc.gid_exists($o3.x[i][j])) {
		sprint(tstr, "%d %d spikes", i, j)		
		if ($o3 == mitgid) {
			pnm.pc.post(tstr, mit[i][j].spiketimes)
		}else{
			pnm.pc.post(tstr, gran[i][j].spiketimes)
		}
	}
    }
  }
}

proc take_raster() {
	if (pnm.pc.gid_exists($o3.x[$1][$2])) {
		if ($o3 == mitgid) {
			tmpspike = mit[$1][$2].spiketimes.c
		}else{
			tmpspike = gran[$1][$2].spiketimes.c
		}
	}else{
		sprint(tstr, "%d %d spikes", $1, $2)
		pnm.pc.take(tstr, tmpspike)
	}
}

proc print_raster() { // 1 arg - filename root

  sprint(filename,"%s.ras",$s1)
  outfile.wopen(filename)
  outfile.printf("# Mitral cell raster plot\n")
  pnm.pc.context("post_raster(nmitx, nmity, mitgid)\n")
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
	take_raster(i, j, mitgid)
      for k = 0, tmpspike.size()-1 {
        outfile.printf("%d %d %d %10.3f\n",i,j,i*nmity+j,tmpspike.x[k])
      }
    }
  }
  outfile.close()
}

proc print_gran_raster() { // 1 arg - filename root
  sprint(filename,"%s.gran.ras",$s1)
  outfile.wopen(filename)
  outfile.printf("# Granule cell raster plot\n")
  pnm.pc.context("post_raster(ngranx, ngrany, grangid)\n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
	take_raster(i, j, grangid)
      for k = 0, tmpspike.size()-1 {
        outfile.printf("%d %d %d %10.3f\n",i,j,i*ngrany+j,tmpspike.x[k])
      }
    }
  }
  outfile.close()
}

proc print_isis() { // 3 args - indices of mitral cell plus filename root
  calc_isis($1,$2,ttrans)
  sprint(filename,"%s_%d_%d.isi",$s3,$1,$2)
  outfile.wopen(filename)
  outfile.printf("# Interspike intervals for mitral cell [%d][%d]",$1,$2)
  isi.printf(outfile,"%10.3f")
  outfile.close()
}

proc print_isi_stats() { // 1 arg - filename root
  sprint(filename,"%s.stats",$s1)
  outfile.wopen(filename)
  outfile.printf("#Interspike interval statistics for mitral cells\n")
  outfile.printf("# i  j   n   mean   median   stdev  \n")
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_isis(i,j,ttrans)
      outfile.printf("%3d%3d%4d",i,j,isi.size())
      if (isi.size() > 0) {
        outfile.printf("%8.2f%8.2f",isi.mean(),isi.median())
        if (isi.size() > 1) {
          outfile.printf("%8.2f\n",isi.stdev())
        } else {
          printf("\n")
        }
      } else { outfile.printf("\n") }
    }
  }
  outfile.printf("#Interspike interval statistics for granule cells\n")
  outfile.printf("# i  j   n   mean   median   stdev  \n")
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
      calc_gran_isis(i,j,ttrans)
      outfile.printf("%3d%3d%4d",i,j,isi.size())
      if (isi.size() > 0) {
        outfile.printf("%8.2f%8.2f",isi.mean(),isi.median())
        if (isi.size() > 1) {
          outfile.printf("%8.2f\n",isi.stdev())
        }
      } else { outfile.printf("\n") }
    }
  }
  outfile.close()
}

proc print_lags() { local i,j // 3 args - indices of mitral cell + filename root
  sprint(filename,"%s_%d_%d.lags",$s3,$1,$2)
  outfile.wopen(filename)
  outfile.printf("# Lag times for mitral cell [%d][%d]\n",$1,$2)
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      calc_lags($1,$2,i,j,ttrans)
      outfile.printf("[%d,%d]",i,j)
      lags.printf(outfile,"%10.3f")
    }
  }
  outfile.close()
}

proc print_si() { // 1 arg - fileroot
  print "Calculating synchronization indices"
  sprint(filename,"%s.synch",$s1)
  outfile.wopen(filename)
  outfile.printf("Synchronization index: %10.3f\n",synch_index())
  outfile.printf("Phase-locking index:   %10.3f\n",phaselock_index())
  outfile.close()
}

Loading data, please wait...