model-2.3.hoc

/************************************************************
Dentate Gyrus model from the Soltesz Lab

For the following details, consult the README.txt file:
Version History
Purpose of code
Input files
Output files
Code Outline

For even more information, consult Instruction Manual.pdf
************************************************************/

loadstart = startsw()					// record the start time of the set up
/***********************************************************************************************
I.  LOAD LIBRARIES & PARAMETERS
***********************************************************************************************/
{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

{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
													//  that can't be changed at command line


/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
{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 

{load_file("./setupfiles/load_cell_conns.hoc")}	// Load in the cell connectivity info

strdef tempFileStr						// Define string reference for the names of the cell template files
proc loadCellTemplates(){local i		// Define one template for each cell type in cells2include, plus perforant path cell(s)
	for i=0, numCellTypes-1 {
		sprint(tempFileStr,"./cells/class_%s.hoc",cellType[i].cellType_string)
		load_file(tempFileStr)			// Load the cell type's class template
	}
}	

loadCellTemplates()

proc calcNetSize(){local i				// Calculate final network size after sclerosis
	totalCells = 0						// Initialize totalCells so we can add to it iteratively in the for loop
	ncell = cellType[0].numCells		// Initialize ncell so we can add to it iteratively in the for loop
	for i=1,numCellTypes-1 {			// Run the following code for non ppstim cell types only
		cellType[i].updateGidRange(cellType[i-1].cellEndGid+1)	// Update the gid ranges for each cell type
		totalCells = totalCells + cellType[i].numCells			// Total number of cells after sclerosis, not including pp cells
		ncell = ncell + cellType[i].numCells 					// Total # cells including pp cells
	}
	random_stream_offset_= (totalCells*2+2)	// How far down the stream to start the next stream section -
											// 	must be higher than the length of each stream section, which
											// 	is the # of times this thing will be called per cell x random type
}
calcNetSize()

proc calcBinSize(){local NumGCells

	for i=0, numCellTypes-1 {		// Using the specified dimensions of the network (in um) and
									//  the total number of cells of each type, set the number
									//  of bins in X, Y, Z dimensions such that the cells will be
									//  evenly distributed throughout the allotted volume
									// just changed this so even the stim cells will be allotted, as now we have some
									// stimulation protocols that incorporate stim cell position
	
		cellType[i].setBins(LongitudinalLength,TransverseLength,LayerVector.x[cellType[i].layerflag])  
									// For the z length, use the height of the layer in which the
									// cell somata are found for this cell type
	}
}
calcBinSize()

/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY
***********************************************************************************************/
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
		}
	}
}


objref  strobj
strobj = new StringFunctions()
strdef direx
if (strcmp(UID,"0")==0 && pc.id==0) {
	type = unix_mac_pc() // 1 if unix, 2 if mac, 3 if mswin, or 4 if mac osx darwin
	if (type<3) {
		{system("uuidgen", direx)} // unix or mac
		strobj.left(direx, strobj.len(direx)-1)
	} else {
		{system("cscript //NoLogo setupfiles/uuid.vbs", direx)} // pc
	}
	UID = direx
}

{load_file("./setupfiles/save_run_info.hoc")}

loadtime = startsw() - loadstart		// Calculate the set up time (now - recorded start time) in seconds
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (set up)\n************\n", loadtime)}
createstart = startsw()					// Record the start time of the cell creation
/***********************************************************************************************
IV. CREATE, UNIQUELY ID, AND POSITION CELLS
***********************************************************************************************/
strdef cmd
objref cells, ranconlist, ransynlist
cells = new List()						
ranconlist = new List()
ransynlist = new List()

{load_file("./setupfiles/position_functions.hoc")}	// Defines the algorithm used to set the
													//  positions of the cells
													
{load_file("./setupfiles/create_cells_pos.hoc")}	// Creates each cell on its assigned host
													//  and sets its position using the algorithm
													//  defined above

strdef cmd
{load_file("./setupfiles/write_positions.hoc")}		// Defines the proc to write the position

if (PrintCellPositions>0) {posout()}	// Calls the proc that writes the position of each cell to the

createtime = startsw() - createstart	// Calculate time taken to create the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (created cells)\n************\n", createtime)}
connectstart = startsw()				// Grab start time of cell connection

objref ic
cellType[3].CellList[int(cellType[3].numCells/2)].soma ic = new IClamp(0.5)
ic.del = 10
ic.dur = 4000
ic.amp = .05

/***********************************************************************************************
V.	CONNECT THE CELLS AND CONNECT THE PERFORANT PATH TO SOME CELLS
***********************************************************************************************/


{load_file("./setupfiles/nc_append_functions.hoc")}	// Defines nc_append and nc_appendo, which 
													//  are procs (?) needed to create the netcon
													//  object associated with each connection
													//  (synapse)

strdef cmd
{sprint(cmd,"./connections/%s_connections.hoc", Connectivity)}
{load_file(cmd)}

{sprint(cmd,"./stimulation/%s_stimulation.hoc", Stimulation)}
{load_file(cmd)}

{load_file("./setupfiles/write_conns_functions.hoc")}
if (PrintConnSummary==1) {printNumConFile()}	// Write a summary file of all cell connections by pre and post cell types "numcons.dat"
if (PrintConnDetails==2) {tracenet()}	// Write the file of all cell connections (pre and post gids, synapse type, host on which cell exists), "connections.dat"


connecttime = startsw() - connectstart			// calculate time taken to connect the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (connected cells)\n************\n", connecttime)}
runstart = startsw()							// grab start time of the simulation
/***********************************************************************************************
VI.	INITIALIZE AND RUN NETWORK, OUTPUT RESULT FILES
***********************************************************************************************/
if (PrintVoltage==1) {load_file("./setupfiles/recordvoltageauto.hoc")}	// If we are going to record some voltages, then set up the recordings

proc init() { local dtsav, temp, secsav, secondordersav	// initialize the simulation
	dtsav = dt						// Save desired dt value to reset after temporarily changing dt
	secondordersav = secondorder	// Save desired secondorder value to reset after temporarily changing secondorder

	finitialize(v_init)	// Call finitialize (since we are replacing the default init proc that calls this)
						// 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 = -500			// Set the start time for (pre) simulation; -500 to prepare network in advance of start at 0
	dt= 10				// Set dt to large value
	secondorder = 0		// Set secondorder to 0 to set the default fully implicit backward euler for numerical integration (see NEURON ref)
		
	temp= cvode.active()
	if (temp!=0) {cvode.active(0)}	// If cvode is on, turn off temporarily to do large fixed step
	// Now, do a large pre run from t = -500 to t = -100 to set the network 'settle' and all components to reach steady state
	while(t<-100) { fadvance() if (PrintTerminal>1) {print t}}	// Integrate all section equations over the interval dt. increment t by dt
															//	and repeat till t at -100
	if (temp!=0) {cvode.active(1)}	// If cvode was on and then turned off, turn it back on now
	
	t = tstart 						// Start time of the simulation
	dt = dtsav						// Reset dt to specified value
	secondorder = secondordersav	// Reset secondorder to specified value
	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
	}
}


{load_file("./setupfiles/write_results_functions.hoc")}		// Defines the procs spikeout() and timeout()
															//  which write out the spikeraster (spikeout)
															//  and the time taken to run each section of
															//  the code (timeout)
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)
												
	pc.set_maxstep(10)							// 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)

	runstart = startsw()							// grab start time of the simulation
	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

	runtime = startsw() - runstart				// Calculate runtime of simulation
												// Print a time summary message to screen
	comptime = pc.step_time
	avgcomp = pc.allreduce(comptime, 1)/pc.nhost
	maxcomp = pc.allreduce(comptime, 2)
	if (maxcomp>0) {
		if (pc.id == 0) { printf("load_balance = %g\n", avgcomp/maxcomp)}
		if (pc.id == 0) { printf("exchange_time = %g\n",  runtime - maxcomp) }
	} else {
		if (pc.id == 0) { printf("no load balance info available\nno spike exchange info available\n")}
	}
	
	if (pc.id == 0) {printf("****\nTIME SUMMARY for host 0\nset up in %g seconds\ncreated cells in %g seconds\nconnected cells in %g seconds\nran simulation in %g seconds\n************\nThis run is called: %s\n************\n", loadtime, createtime, connecttime, runtime, RunName)}
		
	spikeout()	// Write the spike times (in ms) and spiking cells (by gid) into a file called "spikeraster.dat"
	
	timeout()	// Write out a file of run times for each code section

	if (PrintVoltage==1) {voltageout()}	// Write out any voltages that were recorded during simulation
	//if (PrintConnDetails==1) {tracenetfast()}	// Same as tracenet below but each processor writes its own file and they get concatenated if CatFlag==1
	if (PrintConnDetails==1) {tracenet()}	// Write one file of all cell connections (pre and post gids, synapse type,
										//  host on which cell exists), "connections.dat"

	//if ((PrintConnDetails==0) && (PrintVoltage==1)) {tracecellsfast()}	// Same as tracecells below, but each processor writes its own file and they get concat if CatFlag==1
	if ((PrintConnDetails==0) && (PrintVoltage==1)) {tracecells()}	// If detailed connections are not being written out, at least write them out for cells that are being traced
	
	sumnumout(runtime)	// Write out a file with the total number of cells, spikes, and connections (including contributions from artificial cells)

        if (PrintNeuroML2==1) {
            {load_file("./setupfiles/write_neuroml2.hoc")}
            printNeuroML2Network()  // Write out the structure of the network in NeuroML 2 in  network.nml
        }
}
		
rrun()	// Run the network simulation
//{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