netclamp.hoc

/***********************************************************************************************
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