/***********************************************************************************************
I. LOAD LIBRARIES
***********************************************************************************************/
{load_file("nrngui.hoc")} // Standard definitions - NEURON library file
{load_file("netparmpi.hoc")} // Contains the template that defines the properties of
// the ParallelNetManager class, which is used to set
// up a network that runs on parallel processors
{load_file("./setupfiles/ranstream.hoc")} // Contains the template that defines a RandomStream
// class used to produce random numbers
// for the cell noise (what type of noise?)
{load_file("./setupfiles/CellCategoryInfo.hoc")} // Contains the template that defines a
// CellCategoryInfo class used to store
// celltype-specific parameters
{load_file("./setupfiles/defaultvar.hoc")} // Contains the proc definition for default_var proc
// that can't be changed at command line
{load_file("./setupfiles/parameters.hoc")} // Loads in operational and model parameters that can
// be changed at command line
{load_file("./setupfiles/set_other_parameters.hoc")}// Loads in operational and model parameters
/**********************/
default_var("BasedOnRunName","ReserveTrestles_01") // Name of simulation run that the netclamp is based on
default_var("CellOI","axoaxoniccell") // type of presynaptic cell
default_var("CellOItype",0) // type of presynaptic cell
default_var("gid",11) // gid of presynaptic cell
strdef cmdstr, cmd, RunName
sprint(RunName,"%s_%d", CellOI, gid)
/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
celsius=34
{load_file("./setupfiles/load_cell_category_info.hoc")} // Reads the 'cells2include.hoc' file and
// loads the information into one
// 'CellCategoryInfo' object for each cell
// type (bio or art cells?). Info includes
// number of cells, gid ranges, type name
strdef tempFileStr
sprint(tempFileStr,"./cells/class_%s.hoc",cellType[CellOItype].technicalType) // Concatenate the
// path and file
print "loading ", tempFileStr
load_file(tempFileStr) // Load the file with the template that defines the classs
sprint(tempFileStr,"./cells/class_%s.hoc","ppspont") // Concatenate the
print "loading ", tempFileStr
load_file(tempFileStr) // Load the file with the template that defines the class
totalCells = 1 // totalCells counts the number of 'real' cells
// Create numStims number of NetStims
strdef fn
objref f2
f2 = new File()
sprint(fn, "./netclamp/%s/%s_%d_NetStimConns.dat", BasedOnRunName, CellOI, gid)
f2.ropen(fn)
numStims = f2.scanvar
ncell = 1 + numStims // ncell counts all 'real' and 'artificial' cells
objref TMPcellType[numCellTypes+1]
for i=0,numCellTypes-1 {
if (i==CellOItype) {
cellType[CellOItype].numCells = 1
} else {
cellType[i].numCells = 0
}
if (i>0) {
cellType[i].updateGidRange(cellType[i-1].cellEndGid+1) // Update the gid range for each
} else {
cellType[i].updateGidRange(0) // Update the gid range for each
}
TMPcellType[i] = cellType[i]
}
TMPcellType[numCellTypes] = new CellCategoryInfo(numCellTypes) // Make one object for each cell type to store cell type info
TMPcellType[numCellTypes].setCellTypeParams("PatternStim", "ppspont", 1, numStims, 0, 1) // Set parameters
TMPcellType[numCellTypes].numCons = new Vector(numCellTypes,0)
numCellTypes = numCellTypes+1
objref cellType[numCellTypes]
for i=0,numCellTypes-1{
cellType[i] = TMPcellType[i]
print cellType[i].cellType_string, " ", cellType[i].numCells, "(|", cellType[i].cellStartGid, "| to |", cellType[i].cellEndGid, "|)"
}
objref TMPcellType
/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY AND WRITE OUT RUN RECEIPT
***********************************************************************************************/
objref pnm, pc, nc, nil
proc parallelizer() {
pnm = new ParallelNetManager(ncell) // Set up a parallel net manager for all the cells
pc = pnm.pc
pnm.round_robin() // Incorporate all processors - cells 0 through ncell-1
// are distributed throughout the hosts
// (cell 0 goes to host 0, cell 1 to host 1, etc)
}
parallelizer()
iterator pcitr() {local i2, startgid // Create iterator for use as a standard 'for' loop
// throughout given # cells usage:
// for pcitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
// it_start and it_end let you define range over
// which to iterate
// i1 is the index of the cell on the cell list for that host
// i2 is the index of that cell for that cell type on that host
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
i1 = numcycles+addcycle // the index into the cell # on this host.
i2 = 0 // the index of the cell in that cell type's list on that host
if (startgid<=$5) {
for (i3=startgid; i3 <= $5; i3 += pc.nhost) { // Just iterate through the cells on
// this host(this simple statement
// iterates through all the cells on
// this host and only these cells because
// the roundrobin call made earlier dealt
// the cells among the processors in an
// orderly manner (like a deck of cards)
$&1 = i1
$&2 = i2
$&3 = i3
iterator_statement
i1 += 1
i2 += 1
}
}
}
// Create one Real Cell
objref cells, ransynlist, ranstimlist
cells = new List()
ransynlist = new List()
ranstimlist = new List()
func xpos_algorithm() { // Arguments: 1-gid, 2-numCells, 3-startGid, 4-binNum, 5- binNum, 6-binSize; Return: x position of cell
return 0
}
func ypos_algorithm() { // Arguments: gid, numCells, startGid, binNum, binNum, binSize; Return: y position of cell
return 0
}
func zpos_algorithm() { // Arguments: 1-gid, 2-numCells, 3-startGid, 4-binNum, 5-binSize, cell layer Zo; Return: z position of cell
return 0
}
{load_file("./setupfiles/create_cells_pos.hoc")} // Creates each cell on its assigned host
// and sets its position using the algorithm
// defined above
objref cell
for r=cellType[numCellTypes-1].cellStartGid, cellType[numCellTypes-1].cellEndGid {
cell = pc.gid2cell(r)
cell.interval=1e9
cell.start=-1
}
// Add a recording to its soma potential and somatic ion channel currents, as well as a few other locations
strdef nickname, secstr
nickname="soma"
myi_flag=1
strdef cmdstr
strdef mname, tmpstr
objref strobj, mt, cell
objref mechstring[9], mechlength
mechlength = new Vector(1)
strobj = new StringFunctions()
objref cell, mt
cell = pc.gid2cell(cellType[CellOItype].cellStartGid)
secstr="soma"
sprint(cmdstr,"objref mytrace%s", nickname)
{execute1(cmdstr)}
{sprint(cmdstr,"mytrace%s = new Vector(tstop/dt)", nickname)}
{execute1(cmdstr)}
//{sprint(cmdstr,"mytrace%s.record(&cell.%s.v(0.5))", nickname, secstr)}
{sprint(cmdstr,"mytrace%s.record(&%s.v(0.5))", nickname, cell.myroot)}
{execute1(cmdstr) }
/* Get mechanisms from soma section "myroot" recorded */
if (myi_flag==1) {
mt = new MechanismType(0)
objref mechstring[mt.count()]
k = 0
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if( ismembrane(mname)) {
if (strcmp(mname,"capacitance")!=0 && strcmp(mname,"morphology")!=0 && strcmp(mname,"iconc_Ca")!=0 && strcmp(mname,"iconcCa")!=0 && strcmp(mname,"ccanl")!=0 && strcmp(mname,"cad")!=0 && strcmp(mname,"pas")!=0 && strcmp(mname,"vmax")!=0 ) { //
if (strobj.substr(mname,"_ion")==-1) {
//printf("myi_%s \n", mname)
sprint(tmpstr, "myi_%s", mname) // "cell.soma.%s(0.5)", tmpstr
mechstring[k] = new String()
mechstring[k].s = tmpstr
if (strcmp("myi_ch_Nav", mechstring[k].s)==0) {
print nickname, ": ", tmpstr}
{k = k+1}
}
}
}
}
{mechlength.x[0] = k}
if (strcmp("som", nickname)==0) {
print "num mechs: ", k}
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "{objref %s_%s}", nickname, mechstring[rt].s)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s = new Vector(%g)}", nickname, mechstring[rt].s, (tstop-tstart)/dt)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s.record(&%s.%s(0.5))}", nickname, mechstring[rt].s, cell.myroot, mechstring[rt].s)
execute1(cmdstr)
if (strcmp("myi_ch_Nav", mechstring[rt].s)==0) {
print "record: ", cmdstr}
}
}
// Connect NetStims to all CellOI input synapses
{load_file("./setupfiles/nc_append_functions.hoc")} // Defines nc_append and nc_appendo, which
{load_file("./setupfiles/load_cell_conns.hoc")} // Load in the cell connectivity info
pregid = 1
numConnTypes = f2.scanvar
postcellType = CellOItype
for r=0, numConnTypes-1 {
precellType = f2.scanvar
numSynTypes = f2.scanvar
synWeight = cellType[precellType].wgtConns.x[postcellType]
for SynNumber = 0, numSynTypes-1 {
print "pregid = ", pregid, ", precellType = ", precellType, ", synapse = ", SynNumber, ", weight = ", synWeight
nc_append(pregid, cellType[postcellType].cellStartGid, precellType, SynNumber, synWeight + (pregid+1)*1000, 0.5) // Make the connection
pregid += 1
cellType[postcellType].numCons.x[precellType] += 1
}
}
f2.close()
// Read in firing times of all net stims
strdef fn
objref f2
objref pattern_, tvec_, idvec_
f2 = new File()
sprint(fn, "./netclamp/%s/%s_%d_NetStimTimes.dat", BasedOnRunName, CellOI, gid)
f2.ropen(fn)
numspikes = f2.scanvar
pattern_ = new PatternStim()
tvec_ = new Vector(numspikes)
idvec_ = new Vector(numspikes)
for i=0, numspikes-1 {
tvec_.x[i] = f2.scanvar // spike time in ms
idvec_.x[i] = f2.scanvar // gid of NetStim to make fire
}
f2.close()
pattern_.fake_output = 1
pattern_.play(tvec_, idvec_)
// Run simulation & write out spike raster and recording results
proc init() { local dtsav, temp, secsav, secondordersav // Redefine the proc that initializes the
// simulation (why?)
dtsav = dt // Save the desired dt value to reset it after temporarily
// changing dt to run a quick presimulation to allow the
// network to reach steady state before we start 'recording'
secondordersav = secondorder // Save the desired secondorder value to reset it after
// temporarily changing secondorder to run a quick presimulation
// to allow the network to reach steady state before we start
// 'recording' its results
finitialize(v_init) // Call finitialize from within the custom init proc, just as the default
// init proc does. Note that finitialize will call the INITIAL block for
// all mechanisms and point processes inserted in the sections and set the
// initial voltage to v_init for all sections
t = -200 // Set the start time for (pre) simulation; -500 ms to allow the network to
// reach steady state before t = 0, when the real simulation begins
dt= 10 // Set dt to a large value so that the presimulation runs quickly
secondorder = 0 // Set secondorder to 0 to set the solver to the default fully implicit backward
// euler for numerical integration (see NEURON ref)
temp= cvode.active() // Check whether cvode, a type of solver (?) is on
if (temp!=0) {cvode.active(0)} // If cvode is on, turn it off while running the presimulation
while(t<-100) { fadvance() if (PrintTerminal>1) {print t}} // Run the presimulation from t = -500
// to t = -100 (why not 0?) to let the
// network and all its components reach
// steady state. Integrate all section
// equations over the interval dt,
// increment t by dt, and repeat until
// t at -100
if (temp!=0) {cvode.active(1)} // If cvode was on and then turned off, turn it back on now
t = tstart // Set t to the start time of the simulation
dt = dtsav // Reset dt to the specified value for the simulation
secondorder = secondordersav // Reset secondorder to the specified value for the simulation
if (cvode.active()){
cvode.re_init() // If cvode is active, initialize the integrator
} else {
fcurrent() // If cvode is not active, make all assigned variables
// (currents, conductances, etc) consistent with the
// values of the states
}
}
simnum=0
strdef cmdstr
sprint(cmdstr,"mkdir ./netclamp/%s/%s", BasedOnRunName,RunName)
print cmdstr
system(cmdstr)
proc spikeout() {local i, rank localobj f // Write out a spike raster (cell, spike time)
pc.barrier() // Wait for all ranks to get to this point
sprint(cmd,"./netclamp/%s/%s/spikeraster_%03.0f.dat", BasedOnRunName, RunName, simnum)
f = new File(cmd)
if (pc.id == 0) { // Write header to file 1 time only
f.wopen()
f.close()
}
for rank = 0, pc.nhost-1 { // For each processor, allow processor to append to file the spike times of its cells
if (rank == pc.id) { // Ensure that each processor runs once
f.aopen() // Open for appending to file
for i=0, pnm.idvec.size-1 {
f.printf("%.3f %d\n", pnm.spikevec.x[i], pnm.idvec.x[i]) // Print the spike time and spiking cell gid
}
f.close()
}
pc.barrier()
}
}
// also write out vector recordings
objref f, strobj
strobj = new StringFunctions()
strdef outfile, beforestr, afterstr
objref tgt
proc printtrace() {
{sprint(outfile, "./netclamp/%s/%s/%s__%d_%s_trace_%03.0f.dat", BasedOnRunName, RunName, CellOI, gid, secstr, simnum)}
{f = new File(outfile)}
{f.wopen()}
{f.printf("t\tv")} // write header
if (myi_flag==1) {
mt = new MechanismType(0)
objref mechstring[mt.count()]
k = 0
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if( ismembrane(mname)) {
if (strcmp(mname,"capacitance")!=0 && strcmp(mname,"morphology")!=0 && strcmp(mname,"iconc_Ca")!=0 && strcmp(mname,"iconcCa")!=0 && strcmp(mname,"ccanl")!=0 && strcmp(mname,"cad")!=0 && strcmp(mname,"pas")!=0 ) { //
if (strobj.substr(mname,"_ion")==-1) {
//printf("myi_%s \n", mname)
sprint(tmpstr, "myi_%s", mname) // "cell.soma.%s(0.5)", tmpstr
mechstring[k] = new String()
mechstring[k].s = tmpstr
if (strcmp("myi_ch_Nav", mechstring[k].s)==0) {
print nickname, ": ", tmpstr}
{k = k+1}
}
}
}
}
mechlength.x[0]=k
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "%s_%s", nickname, mechstring[rt].s)
if (strcmp("myi_ch_Nav", mechstring[rt].s)==0) {
print "Nav nickname: ", nickname}
if (strcmp("som", nickname)==0) {
print "cmdstr: ", cmdstr, " declared: ", name_declared(cmdstr) }
if (name_declared(cmdstr)>0) {
f.printf("\t%s", mechstring[rt].s)
}
}
}
{f.printf("\n")}
for i=0, (tstop-tstart)/dt-1 {
{sprint(cmdstr, "{f.printf(\"%%g\\t%%g\", i*dt, mytrace%s.x[i])}", nickname)}
{execute1(cmdstr)}
if (myi_flag==1) {
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "%s_%s", nickname, mechstring[rt].s)
if (name_declared(cmdstr)>0) {
sprint(cmdstr, "{f.printf(\"\\t%%g\", %s_%s.x[i])}", nickname, mechstring[rt].s)
execute1(cmdstr)
}
}
}
{f.printf("\n")}
}
{f.close()}
}
use_cache_efficient=1
get_spike_hist=0
use_bin_queue=0
use_spike_compression=0
if (use_spike_compression==1) {
maxstepval = 2.5
} else {
maxstepval = 0.5
}
proc rrun(){ // Run the network simulation and write out the results
pnm.want_all_spikes() // Record all spikes of all cells on this machine into the
// vectors pnm.spikevec (spiketimes) and pnm.idvec (gids)
local_minimum_delay = pc.set_maxstep(maxstepval) // Set every machine's max step size to minimum delay of
// all netcons created on pc using pc.gid_connect, but
// not larger than 10
stdinit() // Call the init fcn (which is redefined in this code) and
// then make other standard calls (to stop watch, etc)
pc.psolve(tstop) // Equivalent to calling cvode.solve(tstop) but for parallel NEURON;
// solve will be broken into steps determined by the result of
// set_maxstep
spikeout() // Write the spike times (in ms) and spiking cells (by gid) into a file called "spikeraster.dat"
printtrace()
}
walltime = startsw()
strdef cmd, cmdo
newtstop = tstop
warningflag=0
{cvode.cache_efficient(use_cache_efficient)} // always double check that this addition does not affect the spikeraster (via pointers in mod files, etc)
if (use_bin_queue==1) {
use_fixed_step_bin_queue = 1 // boolean
use_self_queue = 0 // boolean - this one may not be helpful for me, i think it's best for large numbers of artificial cells that receive large numbers of inputs
{mode = cvode.queue_mode(use_fixed_step_bin_queue, use_self_queue)}
}
if (maxstepval==1) {
nspike = 3 // compress spiketimes or not
gid_compress = 0 //only works if fewer than 256 cells on each proc
{nspike = pc.spike_compress(nspike, gid_compress)}
}
celsius = 34
simnum=0
numsims=10
lowerbound=0 // 5
rrun()
strdef f2n, f3n
objref f2, f3
f2 = new File()
f3 = new File()
sprint(f2n, "./netclamp/%s/%s_%d_NetStimConns.dat", BasedOnRunName, CellOI, gid)
sprint(f3n, "./netclamp/%s/%s_%d_key.dat", BasedOnRunName, CellOI, gid)
f3.wopen(f3n)
for simnum=1, numsims {
//nclist.remove_all
pregid=1
f2.ropen(f2n)
numStims = f2.scanvar
numConnTypes = f2.scanvar
postcellType = CellOItype
f3.printf("%d\t%f\n", simnum, 1 + (simnum-lowerbound)*.1)
for r=0, numConnTypes-1 {
precellType = f2.scanvar
numSynTypes = f2.scanvar
synWeight = cellType[precellType].wgtConns.x[postcellType]
if (strcmp(cellType[precellType].cellType_string,"ca3cell")==0 || strcmp(cellType[precellType].cellType_string,"eccell")==0) {
synWeight = synWeight*(1 + (simnum-lowerbound)*.1)
}
for SynNumber = 0, numSynTypes-1 {
nclist.o(pregid-1).weight = synWeight
print "pregid = ", pregid, ", precellType = ", precellType, ", synapse = ", SynNumber, ", weight = ", synWeight
//nc_append(pregid, cellType[postcellType].cellStartGid, precellType, SynNumber, synWeight + (pregid+1)*1000, 0.5) // Make the connection
pregid += 1
cellType[postcellType].numCons.x[precellType] += 1
}
}
f2.close()
rrun() // Run the network simulation (in proc rrun)
}
f3.close()
///////////////
{pc.runworker()} // Everything after this line is executed only by the host node
// The NEURON website describes this as "The master host returns immediately. Worker hosts
// start an infinite loop of requesting tasks for execution."
{pc.done()} // Sends the quit message to the worker processors, which then quit NEURON
quit() // Sends the quit message to the host processor, which then quits NEURON
