newnetclamp.hoc

/************************************************************
For more information, consult ModelDoc.pdf
************************************************************/

loadstart = startsw()					// record the start time of the set up
/***********************************************************************************************
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/SynStore.hoc")}	// Contains the template that defines a 
													
{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

//{load_file("./setupfiles/check_for_files.hoc")}	// Checks that the necessary chosen files exist
												//  (datasets, connectivity, and stimulation)
												
default_var("gidOI", 0) // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("cellindOI", 0) // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("spkflag", 0) // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("stimflag", 0) // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("resultsfolder", "00001") // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("runname", "None") // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("origRunName", "CutSma_040_11") // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("celltype2use", "axoaxoniccell") // Run Comments  // -c 'strdef RunName ' 'RunName = "Woo"'
												
celltypeOI=cellindOI
//gidOI=0
i=0
ij=0
gid=0
/*strdef origRunName
origRunName = "CutSma_040_11"
PrintConnDetails=2*/
/***********************************************************************************************
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 path2ConnData
sprint(path2ConnData, "../networkclamp_results/%s/%s/conndata.dat", origRunName, resultsfolder)

strdef path2SynData
sprint(path2SynData, "../networkclamp_results/%s/%s/syndata.dat", origRunName, resultsfolder)

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

strdef tempFileStr						// Define a string reference to store the name of the
										//  current cell template file

proc loadCellTemplates(){local i		// Proc to load the template that defines (each) cell class

	//for i=0, numCellTypes-1 {			// Iterate over each cell type in cells2include (and art cells?)
	i = celltypeOI
		sprint(tempFileStr,"./cells/class_%s.hoc",celltype2use)	// Concatenate the
		//sprint(tempFileStr,"./cells/class_%s.hoc",cellType[i].technicalType)	// Concatenate the
																				//  path and file
																				
		load_file(tempFileStr)			// Load the file with the template that defines the class
										//  for each cell type
	//}
	if (stimflag==1) {
		load_file("./cells/class_sintrain.hoc")			// Load the file with the template that defines the class
		//tstop = 100
		//SimDuration = tstop
		print "Using PhasicData, SimDuration = ", SimDuration
	} else {
		load_file("./cells/class_ppvec.hoc")			// Load the file with the template that defines the class
	}
}	
loadCellTemplates()						// Run the newly defined proc

proc calcNetSize(){local i				// Calculate the final network size (after any cell death)
	//cellType[0].numCells = cellType[numCellTypes-1].cellEndGid - cellType[0].cellEndGid + 4*2 // just added now for spontburst
	//cellType[0].updateGidRange(0)

	totalCells = 0						// Initialize totalCells (which counts the number of 'real'
										//  cells) so we can add to it iteratively in the 'for' loop
										
	ncell = cellType[0].numCells		// Initialize ncell (which counts all 'real' and 'artificial'
										//  cells) so we can add to it iteratively in the 'for' loop
										
	for i=1,numCellTypes-1 {			// Run the following code for 'real' cell types only - need a different way of singling out real cells?
	
		if (cellType[i].layerflag==0) {	// For cell types in layer 0, which are susceptible to death
		
			cellType[i].numCells = int(cellType[i].numCells * ((100-PercentCellDeath)/100))
										// Calculate the number of cells surviving after cell loss
			
			if (cellType[i].numCells == 0) {cellType[i].numCells = 1}	// If all cells of one type
																		//  are killed, let 1 live
		}
		cellType[i].updateGidRange(cellType[i-1].cellEndGid+1)	// Update the gid range for each
																//  cell type
		
		totalCells = totalCells + cellType[i].numCells			// Update the total number of cells
																//   after sclerosis, not including
																//   artificial cells

		ncell = ncell + cellType[i].numCells 					// Update the total number of cells
																//   after sclerosis, including
																//   artificial cells
	}
}
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(scLongitudinalLength,scTransverseLength,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 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
				// if (i3 == $5) {iterator_statement}
				iterator_statement
				i1 += 1
				i2 += 1
		}
	}
}

iterator fastitr() {local i2, startgid	// Create iterator for use as a standard 'for' loop
										//  throughout given # cells usage:
										//  for fastitr(&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
				if (i3 == $5) {iterator_statement}
				// 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("cscript //NoLogo setupfiles/uuid.vbs", direx)} // pc
	} else {
		{system("uuidgen", direx)} // unix or mac
		strobj.left(direx, strobj.len(direx)-1)
	}
	UID = direx
}

strdef cmdstr, cmd
//{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
***********************************************************************************************/

objref  memfile, topfile //, pccc
//{pccc = new ParallelContext()}
{zzz=0}
if (pc.id==0) {
	{sprint(cmd,"./networkclamp_results/%s/%s/memory.dat", runname, resultsfolder)}
	{memfile = new File(cmd)}
	{memfile.wopen()}

	{sprint(cmd,"./networkclamp_results/%s/%s/topoutput.dat", runname, resultsfolder)}
	{topfile = new File(cmd)}
	{topfile.wopen()}
}

strdef direx1, direx2, TopCommand

func mallinfo() {local m // this function is a wrapper for the system function mallinfo
        m = nrn_mallinfo(0)
        if (pc.id == 0) {

			if (PrintTerminal>1) {printf("Memory (host 0) - %s: %ld.  Since last call: %ld\n", $s2, m, m - $1)}
			memfile.printf("%s\t%f\t%ld\t%ld\n", $s2, startsw()-loadstart, m, m - $1)	// Print the spike time and spiking cell gid
			
			
			sprint(TopCommand,"top -p `pgrep %s | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1", TopProc)
			//print "TopCommand = ", TopCommand
			if (strcmp(TopProc,"")!=0) { system(TopCommand, direx1)}
			//system("top -p `pgrep special | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx1)
			//system("top -p `pgrep nrniv | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx2)
			//if (strcmp(direx1,direx2)<0) {
			//	topfile.printf("%s\t%s\n", $s2, direx2)	// Print the spike time and spiking cell gid
			//} else {
				topfile.printf("%s\t%s\n", $s2, direx1)	// Print the spike time and spiking cell gid
			//}
		}
        return m
}

{zzz = mallinfo(zzz, "Prior to creating cells")}

objref cells, ransynlist, ranstimlist, raninitlist
cells = new List()						
ransynlist = new List()
ranstimlist = new List()
raninitlist = new List()
													
{load_file("./setupfiles/clamp/netclamp_create_cells_pos.hoc")}	// Creates each cell on its assigned host
													//  and sets its position using the algorithm
													//  defined above
													
													
strdef cmd

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

oldtimeout = pc.timeout(0)

zzz = mallinfo(zzz, "After creating cells")

/***********************************************************************************************
V.	CONNECT THE MODEL CELLS AND CONNECT THE STIMULATION CELLS TO SOME MODEL CELLS
***********************************************************************************************/

celsius=34

{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)

{pc.barrier()}

{load_file("./setupfiles/clamp/list_randfaststim_connections.hoc")}			// Loads in a particular definition of the connectCells function,
							//  depending on the value of the Connectivity string

zzz = mallinfo(zzz, "After connecting cells")
{pc.barrier()}

if (stimflag==0) {
{load_file("./setupfiles/clamp/netclamp_stimulation.hoc")}			// Loads in a particular stimulation protocol, including the spike
							//  vectors for the artificial cells and their connections to the
							//  'real' cells (where do we define the # of art cells and create them?)
} else {
	if (stimflag==1) {
		{load_file("./setupfiles/clamp/netclamp_oscillation_stimulation.hoc")}
	}
}
zzz = mallinfo(zzz, "After defining stimulation")

{load_file("./setupfiles/write_conns_functions.hoc")}	// Defines two funcs (procs?): printNumConFile()
														//  writes the number of synapses between every
														//  combination of pre- and post-synaptic cell type
														//  and tracenet() lists the presynaptic cell type,
														//  postsynaptic cell type, and synapse type for
														//  every connection
														
if (PrintConnSummary==1) {printNumConFile()}	// Write "numcons.dat": a QUICK, SMALL summary file of the # of
										//  connections by pre and post cell types
/*										
if (PrintConnDetails==2) {allCellConnsParallel()}		// Write "connections.dat": a VERY SLOW, LARGE file of all cell
										//	connections (pre and post gids, synapse type, host on which
										//  (target?) cell exists)
	if (PrintConnDetails==1) {allCellConnsParallel()}	// Write the file of all cell connections (pre and post gids, synapse type,
										//  host on which cell exists), "connections.dat"

	if ((PrintConnDetails==0) && (PrintVoltage==2)) {recordedCellConnsParallel()}	// If detailed connections are not being written out, at least write them out for cells that are being traced
*/
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)}

/***********************************************************************************************
VI.	INITIALIZE AND RUN NETWORK, OUTPUT RESULT FILES
***********************************************************************************************/

objref mytrace, myvrec, cell
objref lfpkmatrix

lfpkmatrix = new Vector()

cell = pc.gid2cell(gidOI)

/*********************
inserted stuff right below here:
************/


/* 
Computes xyz coords of nodes in a model cell  whose topology & geometry are defined by pt3d data.
Code by Ted Carnevale.
*/

proc interpxyz() { local ii, nn, kk, xr, nsegs localobj xx, yy, zz, xint, yint, zint, ll, range
    
    // get the data for the section
    nn = $1
    nsegs = $2
    xx = $o3
    yy = $o4
    zz = $o5
    ll = $o6
    xint = $o7
    yint = $o8
    zint = $o9
    
    // to use Vector class's .interpolate() 
    // must first scale the independent variable
    // i.e. normalize length along centroid
    ll.div(ll.x[nn-1])
    
    // initialize the destination "independent" vector
    range = new Vector(nsegs+2)
    range.indgen(1/nsegs)
    range.sub(1/(2*nsegs))
    range.x[0]=0
    range.x[nsegs+1]=1
    
    // length contains the normalized distances of the pt3d points 
    // along the centroid of the section.  These are spaced at 
    // irregular intervals.
    // range contains the normalized distances of the nodes along the 
    // centroid of the section.  These are spaced at regular intervals.
    // Ready to interpolate.
    
    xint.interpolate(range, ll, xx)
    yint.interpolate(range, ll, yy)
    zint.interpolate(range, ll, zz)
}


proc setup_npole_lfp() { local i, j, k, m, n, r, h, s, avglfp, typei, celltype, ex, ey, ez, sx, sy, sz, lfptype, ii, nn localobj myr, dists, lst, xx, yy, zz, ll, xint, yint, zint strdef ks
    
    rho = 333.0 // extracellular resistivity, [ohm cm]
    fdst = 0.1 // percent of distant cells to include in the computation
    
    ex = mypoint.x[0] //$1
    ey = mypoint.x[1] //$2
    ez = mypoint.x[2] //$3
    
    //printf ("host %d: entering setup_npole_lfp\n", pc.id)
    
    // Vector for storing longitudinal and perpendicular distances
    // between recording electrode and compartments
    dists = new Vector(2)
     
	m = 0
	n = 0

	// Relative to the recording electrode position

	forsec cell.all { 
		insert extracellular
		n = n + nseg
	}

	lfpkmatrix = new Vector(n+1)
		//printf ("host %d: npole_lfp: lfpkmatrix.nrow = %d lfpkmatrix.ncol = %d\n", pc.id, lfpkmatrix.nrow, lfpkmatrix.ncol)

	// Iterates over each compartment of the cell
	j = 0
	forsec cell.all { 
		if (ismembrane("extracellular")) {
			nn = n3d()
					
			xx = new Vector(nn)
			yy = new Vector(nn)
			zz = new Vector(nn)
			ll = new Vector(nn)
				
			for ii = 0,nn-1 {
				xx.x[ii] = x3d(ii)
				yy.x[ii] = y3d(ii)
				zz.x[ii] = z3d(ii)
				ll.x[ii] = arc3d(ii)
			}
				
			xint = new Vector(nseg+2)
			yint = new Vector(nseg+2)
			zint = new Vector(nseg+2)

			interpxyz(nn,nseg,xx,yy,zz,ll,xint,yint,zint)
			
			j = 0
			for (x,0) {
				sx = xint.x[j]
				sy = yint.x[j]
				sz = zint.x[j]
				
				// l = L/nseg is compartment length 
				// r is the perpendicular distance from the electrode to a line through the compartment
				// h is longitudinal distance along this line from the electrode to one end of the compartment
				// s = l + h is longitudinal distance to the other end of the compartment
				l = L/nseg
				r = sqrt((ex-sx)*(ex-sx) + (ey-sy)*(ey-sy) + (ez-sz)*(ez-sz))
				h = l/2
				s = l + h
				k = 0.0001 * area(x) * (rho / (4 * PI * l)) * abs(log((sqrt(h*h + r*r) - h) / (sqrt(s*s + r*r) - s)))
				sprint(ks,"%g",k)
				if ((strcmp(ks,"nan") == 0) || (strcmp(ks,"-nan") == 0)) {
					k = 0
				}
				lfpkmatrix.x[j] = k
					
				j = j + 1
			}
		}
	}
} 
setup_npole_lfp()

/************** done inserting *********/


func npole_lfp() { local vlfp
	vlfp = 0 // Initialize the LFP variable
	forsec cell.all { 
		if (ismembrane("extracellular")) {
			j = 0
			for (x,0) {
				vlfp = vlfp + (i_membrane(x) * lfpkmatrix.x[j])
				j = j + 1
			}
		}
	}  //  vlfp = cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)  
   return vlfp
}

objref lfplist, lfptracelist, vec
lfplist = new List()
lfptracelist = new List()

lfp_dt = 0.5 // sampling interval for LFP

objref mytraceidxlist, lfptracelist
proc sample_lfp() { local i, avglfp, startgid, endgid, gid, celltype localobj vec, lst
    
   avglfp = npole_lfp(MaxEDist) // x, y, z position of electrode in microns
    
    if (pc.id == 0) {
        vec = new Vector()
        vec.append(t, avglfp)
        lfplist.append(vec.c)
    }
    cvode.event(t + lfp_dt, "sample_lfp()") // execute sample_lfp lfp_dt ms from now
}



strdef filestr
sprint(filestr, "./networkclamp_results/%s/%s/allsyns.dat", origRunName, resultsfolder)
objref fAll, fSyn
fAll = new File(filestr)
fAll.wopen()
sprint(filestr, "./networkclamp_results/%s/%s/synlist.dat", origRunName, resultsfolder)
fSyn = new File(filestr)
fSyn.wopen()

strdef cmdstr, tempFileStr
objref newSecRef


if (pc.gid_exists(gidOI)) {
	cell = pc.gid2cell(gidOI)
	print "cell is located at: ", cell.x, ", ", cell.y, ", ", cell.z
	mytrace = new Vector((tstop-tstart)/dt)
	mytrace.record(&cell.soma.v(0.5))

	ist = 0
	ien = 0
	myi = 0
	access cell.soma[0]
	{distance()}

	for precellType = 0, numCellTypes-1 {
		ist = myi
		for r=0, cellType[celltypeOI].SynList[precellType].count()-1 {
			sprint(cmdstr,"newSecRef = cell.%s",cellType[celltypeOI].SynList[precellType].object(r).SecRefStr)
			execute(cmdstr) // sets newSecRef
						
			forsec newSecRef {		
				for (x,0) {
					//y=0
					execute(cellType[celltypeOI].SynList[precellType].object(r).CondStr)
					 if (y==1) {
						fSyn.printf("%s %s %g %s\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, myi, secname())
						myi = myi + 1
					}
				}
			}
		}
		ien = myi - 1
		if (ien>=ist) {fAll.printf("%s %s %g %g\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, ist, ien)}
	}		

	fAll.close()
	fSyn.close()
	
	/*forsec cell.all {
		insert extracellular
		insert xtra
	}

	myvrec = new Vector((tstop-tstart)/dt)
	myvrec.record(&lfp_xtra)*/
} else {
	print "couldnt make cell ", gidOI
}
/* 
load_file("clamp/interpxyz.hoc")	// only interpolates sections that have extracellular
load_file("clamp/setpointers.hoc")	// automatically calls grindaway() in interpxyz.hoc
load_file("clamp/calcrxc.hoc")	// computes transfer r between segments and recording electrodes
*/



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>2) {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
	}
	frecord_init() // see email from ted - fadvance() increments the recorder, so we need to fix the index it ends up at
}


{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)
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 = 10
}														
															
proc rrun(){									// Run the network simulation and write out the results

	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
	if (pc.id==0 && use_spike_compression==1) {print "Host ", pc.id, " has local_minimum_delay=", local_minimum_delay}


	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
	writestart = startsw()
	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 = %.3f\n", avgcomp/maxcomp)}
		if (pc.id == 0) { printf("exchange_time = %.2f ms\n",  runtime - maxcomp) }
	} else {
		if (pc.id == 0) { printf("no load balance info available\nno spike exchange info available\n")}
	}
	
	/*
	spikeoutfast()
	
	TimeOutSerial()	// Write out a file of run times for each code section
	
	if (PrintVoltage>0) {voltageout()}	// Write out any voltages that were recorded during simulation
	if (PrintConnDetails==1) {allCellConnsParallel()}	// Write the file of all cell connections (pre and post gids, synapse type,
										//  host on which cell exists), "connections.dat"

	if ((PrintConnDetails==0) && (PrintVoltage==1)) {recordedCellConnsParallel()}	// If detailed connections are not being written out, at least write them out for cells that are being traced
	
	SummaryOut(runtime)	// Write out a file with the total number of cells, spikes, and connections (including contributions from artificial cells)
	*/
	
	writetime = startsw() - writestart				// Calculate write time of program
	totaltime = startsw() - loadstart
	allottedtime = JobHours*3600
	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\nwrote results in %g seconds\nTOTAL TIME   = %g seconds\nALLOTTED TIME = %g seconds\n************\nThis run is called: %s %s\n************\n", loadtime, createtime, connecttime, runtime, writetime, totaltime, allottedtime, origRunName, resultsfolder)}
}

objref fih
// ComputeNpoleLFP=1 ComputeDipoleLFP=0
if (ComputeNpoleLFP > 0) {
	// execute sample_lfp() at t = 0,
	// right after the mechanism INITIAL blocks have been executed
	fih = new FInitializeHandler("sample_lfp()")
}

objref fihw
fihw = new FInitializeHandler(2, "midbal()")
walltime = startsw()
strdef cmd, cmdo
newtstop = tstop
warningflag=0
proc midbal() {local wt, thisstep
	wt = startsw()
	if (t>0) {
		thisstep = wt - walltime
		simleft = tstop - t
		compleft = JobHours*3600 - (startsw() - loadstart)
		//print "t=",t,", tstop=",tstop,", simleft=",simleft,", jobtime=",JobHours*3600,", timeused=",startsw() - loadstart,", compleft=", compleft 
		if (warningflag==0 && (simleft/StepBy*thisstep+EstWriteTime)>compleft && pc.id == 0) {
			newtstop = int((compleft-EstWriteTime)/thisstep)*StepBy + t
			print "Not enough time to complete ", tstop,  " ms simulation, simulation will likely stop around ", newtstop, " ms"
			warningflag=1
		}
		if (pc.id == 0) { printf("%g ms interval at t=%g ms was %g s\n", StepBy, t, thisstep) }
		
		if ((thisstep+EstWriteTime)>compleft && t<tstop) { // not enough time for another step, end now
			{pc.barrier()}
			sprint(cmd,"cat results/%s/runreceipt.txt | sed -e 's/^SimDuration = [^ ]*/SimDuration = %g;/' > x ; mv x results/%s/runreceipt.txt", RunName, t, RunName)
			if (pc.id == 0) { {system(cmd,cmdo)} }
			if (pc.id == 0) { print "simulation stopped early at t=", t, " ms"}
			tstop = t
			stoprun=1
		}
	}
	walltime = wt
	cvode.event(t+StepBy, "midbal(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 (use_spike_compression==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)}
}

objref spkhist
if (get_spike_hist==1) {
	spkhist = new Vector(pc.nhost)
	if (pc.id==0) {
		pc.max_histogram(spkhist)
	}
}

zzz = mallinfo(zzz, "Before running simulation")

{rrun()}	// Run the network simulation (in proc rrun)

zzz = mallinfo(zzz, "After running simulation")

if (pc.id==0) {
	{memfile.close()}
	{topfile.close()}
}

objref f4
strdef cmd
if (pc.id==0 && get_spike_hist==1) {
	sprint(cmd,"./results/%s/spkhist.dat", RunName)
	f4 = new File(cmd)
	f4.wopen()

			// Open for appending to file
	for i=0, pc.nhost-1 {
		f4.printf("%d\t%d\n", i, spkhist.x[i])	// Print the spike time and spiking cell gid
	}
	f4.close()
}

strdef outfile
objref f
if (pc.gid_exists(0)) {	// If cell exists on this machine
	sprint(outfile, "./networkclamp_results/%s/%s/mytrace_%d_soma.dat", runname, resultsfolder, gidOI)
	//ssprint(outfile, "./results/%s/mytrace_%d_soma.dat", RunName, gidOI)
	f = new File(outfile)
	f.wopen()
	f.printf("t\tv\n")
	for i=0, (tstop-tstart)/dt-1 {
		f.printf("%g\t%g\n", i*dt, mytrace.x[i])
	}
	f.close()
/*
	sprint(outfile, "../networkclamp_results/%s/%s/myvrec_%d_soma.dat", RunName, resultsfolder, gidOI)
	f = new File(outfile)
	f.wopen()
	f.printf("t\tv\n")
	for i=0, (tstop-tstart)/dt-1 {
		f.printf("%g\t%f\n", i*dt, myvrec.x[i])
	}
	f.close()*/

}

proc SpikeOutSerial() {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,"./networkclamp_results/%s/%s/spikeraster.dat", runname, resultsfolder)
	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()
	}
}
SpikeOutSerial()

proc allCellConnsSerial() {local i, rank, srcid localobj tgt, f	// Write out the connections list, including synapse type
	pc.barrier()									// Wait for all ranks to get to this point
	sprint(cmd,"./networkclamp_results/%s/%s/connections.dat", runname, resultsfolder)
	f = new File(cmd)
	if (pc.id == 0) { 								// Write header to file 1 time only
		f.wopen()
		f.printf("source\ttarget\tsynapse\n")
		f.close()
	}
	for rank = 0, pc.nhost-1 {				// For each processor, allow processor to append to file its connection info
		if (rank == pc.id) {				// Ensure that each processor runs once
			f.aopen() 						// Open for appending to file
			for i=0, nclist.count -1 {		// For each connection in the list
				srcid = nclist.o(i).srcgid()	// Get the gid of the source cell
				tgt = nclist.o(i).syn			// Get a reference to the target cell
				f.printf("%d\t%d\t%d\n", srcid, tgt.cid, tgt.sid)	// Print source gid, target gid, synapse id
			}
			f.close()
		}		
		pc.barrier()
	}
}
allCellConnsSerial()

proc writeCells() {
	pc.barrier()					// Wait for all ranks to get to this point
	sprint(cmd,"./networkclamp_results/%s/%s/celltype.dat", runname, resultsfolder)
	f = new File(cmd)
	if (pc.id == 0) { 				// Write header to file 1 time only
		f.wopen()
		f.printf("celltype\ttypeIndex\trangeStart\trangeEnd\n")
		for i=0, numCellTypes-1 {
			f.printf("%s\t%s\t%d\t%d\t%d\n", cellType[i].cellType_string, cellType[i].technicalType, i, cellType[i].cellStartGid, cellType[i].cellEndGid)
		}
		f.close()
	}
}
writeCells()


proc writeLFP() { localobj f, lfpvec, lfptrace
	sprint(cmd,"./networkclamp_results/%s/%s/lfp.dat", runname, resultsfolder)
    f = new File(cmd)
    f.wopen()
    
    // Open for appending to file
    for i=0, lfplist.count()-1 {
        lfpvec = lfplist.o(i)
		f.printf("%g\t%g\n", lfpvec.x[0], lfpvec.x[1])	// Prints time and average LFP value
    }
    f.close()
 }
writeLFP()

{pc.barrier}
{quit()}

/*
// Old quit method
{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
*/