/************************************************************
Current version:
electrophys-1.0
separate stim protocols MC 2010-10-19
************************************************************/
loadstart = startsw() // record the start time of the set up
/***********************************************************************************************
I. LOAD LIBRARIES
***********************************************************************************************/
{load_file("nrngui.hoc")} // Standard definitions - NEURON library file
{load_file("stdlib.hoc")} // Standard library, used by the default_var procedure
{load_file("setupfiles/defaultvar.hoc")} // Contains the proc definition for default_var proc
default_var("relpath","/home/casem/repos/ca1")
default_var("toolpath","/home/casem/matlab/work/RunOrganizer/tools")
default_var("resultspath","cellclamp_results")
default_var("sl","/ ")
objref strobj, fAll, fSyn
strobj = new StringFunctions()
sllength = strobj.len(sl)
strobj.left(sl,sllength-1)
strdef pathstr
{load_file("netparmpi.hoc")} // Contains the template that defines the properties of the
// ParallelNetManager class
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "ranstream.hoc")}
{load_file(pathstr)} // Defines a RandomStream class used to produce random numbers
// for the cell noise
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "CellCategoryInfo.hoc")}
{load_file(pathstr)}// Defines a CellCategoryInfo class used to store
// celltype-specific parameters
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "SynStore.hoc")}
{load_file(pathstr)} // Contains the template that defines a holder of
// synapse type info
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "parameters.hoc")}
{load_file(pathstr)} // Loads in operational and model parameters
//sprint(pathstr,"%s%s", relpath, "setupfiles/loadbal.hoc")
//{load_file(pathstr)} // Loads in a LoadBalance class
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "set_other_parameters.hoc")}
{load_file(pathstr)}// Loads in operational and model parameters
// that can't be changed at command line
{sprint(pathstr,"%s%s%s%s%s%s%s", relpath, sl, "setupfiles", sl, "clamp", sl, "class_pptype.hoc")}
{load_file(pathstr)}
strdef mversion
{mversion = "electrophys"}
/***********************************************************************************************
V? SINGLE CHANNEL RECORDINGS
***********************************************************************************************/
{default_var("onecell",0)} // controls whether to perform single cell recordings
{default_var("cellmorph",0)} // controls whether to plot cell shape
{default_var("pair",0)} // controls whether to perform paired cell recordings
{default_var("mech",0)} // controls whether to perform ion channel recordings
{default_var("conn",100)} // for paired recording, use these synapse weights
{default_var("mysyn",120)} // for paired recording, use these synapse weights
{default_var("cellnum",100)} // ?
{default_var("cellmechs",0)} // ?
{default_var("numeapair",1)} // for paired recording, run this number of pairs and average them
objref cl, f, f2, myvec, mygvec
strdef myfile, channel, chan, conductstr, revpotstr, cmdstr
strdef pathstr, clamppath
//{clamppath = "tools/clamp/"}
{sprint(clamppath,"%s%s%s%sclamp%s", relpath, sl, resultspath, sl, sl)}
{f2 = new File()}
//if (mech==1) {
{sprint(pathstr,"%s%ssetupfiles%sclamp%smechclamp.dat", relpath, sl, sl, sl)}
{f2.ropen(pathstr)}
// mechact = f2.scanvar // 16.8
// mechiv = f2.scanvar // 16.8
somadim = f2.scanvar // 16.8
Ratmp = f2.scanvar // 210
cmtmp = f2.scanvar // 1 // 32*10^(-5) ? // for the p-type channel, 32 pF
Catmp = f2.scanvar // 1 // 32*10^(-5) ? // for the p-type channel, 32 pF
celsius = f2.scanvar // 37
mydt = f2.scanvar // 37
//insert pas
//}
if (mech==1) { // for some reason this creation of soma has to happen in a separate if statement from setting its dimensions
create soma
}
if (mech==1) {
soma {L=somadim diam=somadim}
Ra = Ratmp
cm = cmtmp
/*** New Calcium Stuff ***/
insert iconc_Ca
catau_iconc_Ca = 10
caiinf_iconc_Ca = Catmp //5.e-6
/*** End New Calcium Stuff ***/
stepto_stepsize = f2.scanvar
stepto_numsteps = f2.scanvar
stepto_amplitude1 = f2.scanvar
stepto_amplitude2 = f2.scanvar
stepto_duration1 = f2.scanvar
stepto_duration2 = f2.scanvar
hold_numsteps = f2.scanvar
hold_amplitude1 = f2.scanvar
hold_amplitude2 = f2.scanvar
hold_stepsize = f2.scanvar
hold_duration1 = f2.scanvar
hold_duration2 = f2.scanvar
while (f2.eof()==0) {
f2.scanstr(chan)
{sprint(channel,"ch_%s",chan)}
density = f2.scanvar
// f2.scanstr(conductstr)
conductstr = "gmax"
{f2.scanstr(revpotstr)}
revpot = f2.scanvar
{sprint(cmdstr,"soma {insert %s}", channel)}
{execute1(cmdstr)}
print cmdstr
{sprint(cmdstr,"{%s_%s = %f}", conductstr, channel, density)}
{execute1(cmdstr)}
{sprint(cmdstr,"{%s = %f}", revpotstr, revpot)}
{execute1(cmdstr)}
soma cl = new SEClamp(0.5)
for i = 0, stepto_numsteps-1 { // 16 for p
cl.amp1 = stepto_amplitude1 // mV
cl.dur1 = stepto_duration1
cl.dur2 = stepto_duration2
cl.amp2 = stepto_amplitude2 + i*stepto_stepsize // up to 10 for act/inact, up to 50 for IV curve, -110 for p
tstop = stepto_duration1 + stepto_duration2
dt = mydt
myvec = new Vector(tstop/dt)
mygvec = new Vector(tstop/dt)
{sprint(cmdstr,"{myvec.record(&soma.myi_%s(0.5))}", channel)}
{execute1(cmdstr)}
{sprint(cmdstr,"{mygvec.record(&soma.g_%s(0.5))}", channel)}
{execute1(cmdstr)}
//finitialize(v_init)
stdinit()
run()
// print
f = new File()
sprint(myfile,"%s%s%sstepto_%g.dat", clamppath, chan, sl, cl.amp2)
f.wopen(myfile)
f.printf("t\ti\tg\n")
for j=0, tstop/dt - 1 {
{f.printf("%g\t%f\t%f\n", j*dt, myvec.x[j], mygvec.x[j])}
}
{f.close}
}
for i = 0, hold_numsteps-1 { // 13 for p-type
{cl.amp1 = hold_amplitude1 + i*hold_stepsize} // mV, P-type: -80 + i*10
{cl.dur1 = hold_duration1}
{cl.dur2 = hold_duration2}
{cl.amp2 = hold_amplitude2} // up to 10 for act/inact, up to 50 for IV curve
{tstop = hold_duration1 + hold_duration2}
{dt = mydt}
{myvec = new Vector(tstop/dt)}
{mygvec = new Vector(tstop/dt)}
{sprint(cmdstr,"{myvec.record(&soma.myi_%s(0.5))}", channel)}
{execute1(cmdstr)}
{sprint(cmdstr,"{mygvec.record(&soma.g_%s(0.5))}", channel)}
{execute1(cmdstr)}
{stdinit()}
{run()}
// print
//sprint(dircmd, "mkdir %s%s", , clamppath, chan)
//{system(dircmd, direx)}
{f = new File()}
{sprint(myfile,"%s%s%shold_%g.dat", clamppath, chan, sl, cl.amp1)}
{f.wopen(myfile)}
{f.printf("t\ti\tg\n")}
for j=0, tstop/dt - 1 {
{f.printf("%g\t%f\t%f\n", j*dt, myvec.x[j], mygvec.x[j])}
}
{f.close}
}
// delete soma
{sprint(cmdstr,"soma {uninsert %s}", channel)}
{execute1(cmdstr)}
}
}
{f2.close()}
if ((onecell+pair+cellmorph+cellmechs)==0) {
quit()
}
/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
objref dentateZLength
{tstart = 0} // Start time of simulation
{tstop = mytstop} //1000+2*50 //
{dt = mydt} // Integration interval for fadvance (see NEURON reference)
{myi_flag=1}
strdef newstr
objref f2, f2c // Define object reference for the cells2include file
//objref cellnumvar, celltypestring[1], ReplaceType[1], OrigType[1], cellType[1] // Define placeholder objects, redefine with correct size in fcn
//double cellLayerflag[1], numCellTypes[1] // Define placeholder doubles, redefine with correct size in fcn
//proc loadCellCategoryInfo() {local i, startpos // Load celltype info into a CellCategoryInfo object, 1 per cell type
// f2 = new File()
// sprint(pathstr,"tools/clamp/icell.dat")
// f2.ropen(pathstr)
// numCellTypes = f2.scanvar // # cell types, including 1 for pp cells
// cellnumvar = new Vector(numCellTypes)
// objref celltypestring[numCellTypes], cellType[numCellTypes] // Define variables to temporarily hold data scanned from file
// double cellLayerflag[numCellTypes]
//
// startpos = 1 // 0 is reserved for the pptype which excites all the pairs
// numcellea = 2
//
// // ppspont
// //celltypestring[0]= new String()
// //celltypestring[0].s = "pptype" // Scan in the cell name
// //cellnumvar.x[0]=numCellTypes-1 // make one presynaptic cell per cell type
//
// //cellLayerflag[i]=0 // Scan the layer flag (hilar=2, granular=1, molecular=0), where hilar cells
//
//
// for i=0, numCellTypes-1 {
// celltypestring[i]= new String()
// f2.scanstr(celltypestring[i].s) // Scan in the cell name
// cellnumvar.x[i]=numcellea // One cell for the single cell simulation, one for the pairs
// cellLayerflag[i]=0 // Scan the layer flag (hilar=2, granular=1, molecular=0), where hilar cells
// // are subject to sclerosis
// }
// f2.close()
//
// f2c = new File()
// sprint(cmdstr, "%sconnections/cellnumbers_%g.dat", relpath, cellnum)
// f2c.ropen(cmdstr) // Open the celltype
// numCellTypesLook = f2c.scanvar // Scan the first line, which contains a number giving the
// objref ReplaceType[numCellTypesLook], OrigType[numCellTypesLook]
// for i=0, numCellTypesLook-1 {
// ReplaceType[i]= new String()
// OrigType[i] = new String()
// f2c.scanstr(ReplaceType[i].s) // Scan in the cell name
// f2c.scanstr(OrigType[i].s)
// tmpvar=f2c.scanvar
// tmpvar=f2c.scanvar
// }
// f2c.close()
//
// for i=0, numCellTypes-1 {
// cellType[i] = new CellCategoryInfo(i) // Make one object for each cell type to store cell type info
//
// for j=0, numCellTypesLook-1 {
// if (strcmp(celltypestring[i].s,OrigType[j].s)==0) {
// newstr = ReplaceType[j].s
// }
// }
//
// newstr = celltypestring[i].s
// if (strcmp(celltypestring[i].s,"poolosyncell")==0 || strcmp(celltypestring[i].s,"sprustoncell")==0) {
// newstr = "pyramidalcell"
// }
// if (strcmp(celltypestring[i].s,"ppca3sin")==0) {
// newstr = "ca3cell"
// }
// if (strcmp(celltypestring[i].s,"ppecsin")==0) {
// newstr = "eccell"
// }
//
// print i, ": ", newstr, " (", celltypestring[i].s, ")"
// cellType[i].setCellTypeParams(newstr, celltypestring[i].s, startpos, numcellea, 0, 0) // Set parameters
// cellType[i].numCons = new Vector(numCellTypes,0)
// startpos = startpos + cellType[i].numCells // just make one of each cell
// }
//}
{sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "load_cell_category_info.hoc")}
{load_file(pathstr)}
//loadCellCategoryInfo()
strdef outfile
ConnData = conn
sprint(outfile, "%s%ssetupfiles%sload_cell_conns.hoc", relpath, sl, sl)
{load_file(outfile)} // Load in the cell connectivity info (weight)
SynData = mysyn
sprint(outfile, "%s%ssetupfiles%sload_cell_syns.hoc", relpath, sl, sl)
{load_file(outfile)} // Load in the cell connectivity info (weight)
strdef cmdstr
f2 = new File()
sprint(pathstr,"%s%ssetupfiles%sclamp%sicell.dat", relpath, sl, sl, sl)
f2.ropen(pathstr)
numCellidxes = f2.scanvar // # cell types, including 1 for pp cells
if (numCellidxes>0) {
sprint(cmdstr,"objref cellRectypestring[numCellidxes]")
} else {
sprint(cmdstr,"objref cellRectypestring[1]")
}
execute(cmdstr)
if (numCellidxes>0) {
for i=0, numCellidxes-1 {
cellRectypestring[i]= new String()
f2.scanstr(cellRectypestring[i].s) // Scan in the cell name
}
}
f2.close()
print "pair = ", pair
objref cellPost[1]
strdef pre_type
if (pair>0) {
f2 = new File()
{sprint(pathstr,"%s%ssetupfiles%sclamp%spairtypes.txt", relpath, sl, sl, sl)}
f2.ropen(pathstr)
numpairs = f2.scanvar
if (numCellidxes>0) {
sprint(cmdstr,"objref cellPost[numpairs]")
} else {
sprint(cmdstr,"objref cellPost[1]")
}
execute(cmdstr)
if (numpairs>0) {
objref cellPost[numpairs]
for r = 0, numpairs - 1 {
f2.scanstr(pre_type)
cellPost[r] = new String()
f2.scanstr(cellPost[r].s)
}
}
f2.close()
}
startpos = 1 // 0 is reserved for the pptype which excites all the pairs
numcellea = 2
objref mycellvecs
mycellvecs = new Vector(numCellidxes)
idxes=0
for i=0, numCellTypes-1 {
// cellType[i].numCells = numcellea
// cellType[i].updateGidRange(startpos)
// startpos = startpos + cellType[i].numCells // just make one of each cell
myfl=0
for j=0, numCellidxes-1 {
if (strcmp(cellType[i].cellType_string,cellRectypestring[j].s)==0) {
mycellvecs.x[idxes] = i
idxes = idxes+1
myfl=1
}
}
if (myfl==1) {
cellType[i].numCells = numcellea
} else {
if (pair>0) {
for j=0, numpairs-1 {
if (strcmp(cellType[i].cellType_string,cellPost[j].s)==0) {
myfl=1
}
}
if (myfl==1) {
cellType[i].numCells = numcellea
} else {
cellType[i].numCells = 0
}
} else {
cellType[i].numCells = 0
}
}
cellType[i].updateGidRange(startpos)
startpos = startpos + cellType[i].numCells // just make one of each cell
}
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)
ncell = 1
totalCells = 0
for i=0, numCellTypes-1 {
sprint(tempFileStr,"%s%scells%sclass_%s.hoc", relpath, sl, sl, cellType[i].technicalType)
load_file(tempFileStr) // Load the cell type's class template
print "we'll make ", cellType[i].numCells, " ", cellType[i].cellType_string, "s"
totalCells = numCellTypes + cellType[i].numCells
ncell = ncell + cellType[i].numCells
}
for i=0, numCellTypes-1 {
cellType[i].layerflag=cellType[i].cellStartGid+1
}
random_stream_offset_= (totalCells*2+2)
}
loadCellTemplates()
/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY
***********************************************************************************************/
loadstart = startsw()
strdef iteratorflag
objref pnm, pc, nil, nc
print "ncell: ", ncell
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
}
parallelizer()
iterator pcitr() {local i2, startgid, startgididx, endgididx // Create iterator for use as a standard 'for' loop throughout given # cells
// usage:
// for pcitr(&i, &reli, &gid, it_start, it_end) {do stuff} // i = index of this cell in entire cell list on this host
// reli = relative index: index of this cell in celltype-specific list on this host
// gid = gid of this cell
// it_start and it_end let you define range over which to iterate
i1 = 0
// eventually switch this pcitr to a vector as well so that it is the same as above.
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
if (startgid<=$5) {
for (i2=startgid; i2 <= $5; i2 += pc.nhost) { // Just iterate through the cells on this host
// (this works because of the roundrobin call made earlier)
$&1 = numcycles+addcycle+i1
$&2 = i1
$&3 = i2
iterator_statement
i1 += 1
}
}
}
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()
ransynlist.append(new RandomStream(1, 0)) // placeholder randomstream for the pptype cell, which doesn't need it
objref precell
{precell = new pptype(0,0)}
cells.append(precell)
{precell.connect_pre(nil, nc)} // Create an empty connection for use by the spike detector
{pc.cell(0, nc)} // Associate the cell with its gid and its spike generation location
chdir("cells/")
//print "current dir: ", getcwd()
proc createCells(){ local i, ij, gotil, startat, reli, si, pci, cellind, cnum, runresult, gid // Create cells and assign a GID to each cell
for cellind=0, numCellTypes-1 {
if ($1==1) {
startat = cellType[cellind].cellStartGid
gotil = cellType[cellind].cellStartGid
} else {
startat = cellType[cellind].cellStartGid //+1
gotil = cellType[cellind].cellEndGid
}
for pcitr(&i, &ij, &gid, startat, gotil) {// use the pciter over all cells of this type
print "type: ", cellType[cellind].technicalType, " ", startat, " ", gotil
if (pc.gid_exists(gid)) {
sprint(cmd, "cellType[%g].CellList[%g]=new %s(%g,%g,%g)", cellind, ij, cellType[cellind].technicalType, gid, i, cellind) //+cellType[cellind].cellStartGid) // why add the startgid to the gid?
//print cmd
{runresult=execute1(cmd)} // This command was written as a string so
// the cell object doesn't have to be hard coded
cells.append(cellType[cellind].CellList[ij]) // Append each cell to cells list
ransynlist.append(new RandomStream(1, gid)) // Create a new random number generator for each cell,
{cellType[cellind].CellList[ij].connect_pre(nil, nc)} // Create an empty connection for use by the spike detector
{pc.cell(gid, nc)} // Associate the cell with its gid and its spike generation location
if (cellind>-1 && cellType[cellind].CellList[ij].is_art==0) { // For non ppstim cells, assign position, initialize synapse cid and sid
for si=0, cellType[cellind].CellList[ij].pre_list.count-1 { // Iterate over each pre cell type's synapse list
for j=0, cellType[cellind].CellList[ij].pre_list.o(si).count-1 { // Iterate through each synapse in the list
{cellType[cellind].CellList[ij].pre_list.o(si).o(j).cid=gid} // Set the cell id for each synapse
// Note: Parameters added to Syn2Gid mechanism
}
}
if ((cnum%int(cellType[cellind].numCells/10+1) == 0) && (PrintTerminal>1)) {print cellType[cellind].cellType_string, ": ", reli}
cellType[cellind].CellList[ij].position(gid*5, 0, 0)
}
}
}
}
nc = nil // Then clear the reference to the netcon object, which should destroy the netcon (because all refs would have been removed)
if (PrintTerminal>0) {print "Host ", pc.id, " created cells."}
}
createCells(2)
chdir("../")
objref cell
if (pc.gid_exists(1)) {
cell = pc.gid2cell(1)}
strdef mname, tmpstr
objref strobj
objref mechstring[9], mechlength
proc print_mechs(){local i, k localobj mt, cell // this code courtesy of Jose Ambros-Ingerson via the NEURON forum
mechlength = new Vector(1)
strobj = new StringFunctions()
cell = pc.gid2cell($1)
access cell.soma {
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_CaZ")!=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
{k = k+1}
}
}
}
}
}
//{printf("mt count = %g, mechlength0 = %g mechstring = %s\n", mt.count(), k-1, mechstring[k-1].s)}
{mechlength.x[0] = k}
//{printf("'bout to leave. mechstring count = %g \n", mechstring.count)}
//return mechstring
}
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
/***********************************************************************************************
V? SINGLE CELL DATA
***********************************************************************************************/
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
}
frecord_init() // see email from ted - fadvance() increments the recorder, so we need to fix the index it ends up at
}
proc rrun(){ // Run the network simulation
//pnm.want_all_spikes()
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)
pc.psolve(tstop) // Equivalent to calling cvode.solve(tstop), for parallel NEURON, where solve will be broken into steps determined by the result of set_maxstep
}
objref f2, f, stimAMPvector, cell, strobj //, mechstring[5]
strobj = new StringFunctions()
//local i, sec, duration
strdef mname, headstr, outfile, cmdstr
objref mt, cell, strobj
strobj = new StringFunctions()
if (cellmechs==1) {
for idxes=0, numCellidxes-1 {
cellind = mycellvecs.x[idxes]
gid=cellType[cellind].cellStartGid
if (pc.gid_exists(gid)) {
// add recording locations to cell:
cell = pc.gid2cell(gid)
if (cell.is_art==0) {
sprint(outfile, "%s%s%s%scelldata_%s.dat", relpath, sl, resultspath, sl, cellType[cellind].cellType_string)
f = new File(outfile)
f.wopen()
access cell.soma // access cell.soma[0]
distance() // set the origin for measuring the distance
sprint(headstr,"secname\tpos\tx\ty\tz\tdist\tdiam\tcm\tRa")
mt = new MechanismType(0)
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if (strobj.substr(mname,"ch_")>-1) {
sprint(headstr,"%s\tgmax_%s", headstr, mname)
}
}
sprint(headstr,"%s\n", headstr)
f.printf(headstr)
forsec cell.all { // forsec cell.all // forall
for (x) {
xdist = distance(x)
sprint(headstr,"%s\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g", secname(), x, x3d(0), y3d(0), z3d(0), xdist, diam, cm, Ra)
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if (strobj.substr(mname,"ch_")>-1) {
g=0
e=0
if ( ismembrane(mname)) {
sprint(cmdstr, "g = gmax_%s", mname)
execute1(cmdstr)
}
sprint(headstr,"%s\t%g", headstr, g)
}
}
sprint(headstr,"%s\n", headstr)
f.printf(headstr)
}
}
f.close()
}
}
}
}
objref zz, fcell
tmp1unused=0
tmp2unused=1
strdef cmdstr
if (onecell==1) {
// single cell props
sprint(cmdstr,"%s%s%s%scellproperties.dat", relpath, sl, resultspath, sl)
fcell = new File(cmdstr)
fcell.wopen()
fcell.printf("CellType\tInputResistance\n")
for idxes=0, numCellidxes-1 {
cellind = mycellvecs.x[idxes]
gid=cellType[cellind].cellStartGid
if (pc.gid_exists(gid)) {
// add recording locations to cell:
strdef sname, cmdstr, outfile
cell = pc.gid2cell(gid)
if (cell.is_art==0) {
stdinit()
zz = new Impedance() // this little chunk of code modified from Ted & Michael's MRF tutorial
sprint(cmdstr,"%s zz.loc(0.5) ", cell.myroot)
execute(cmdstr)
zz.compute(0) // DC input R
sprint(cmdstr,"%s { rn = zz.input(0.5) }", cell.myroot) // Input resistance in MegaOhms
execute(cmdstr)
fcell.printf("%s\t%f\n",cellType[cellind].cellType_string, rn)
}
}
}
// what else would we like to see? area of cell (and possibly cell parts)
// somatic conductance
// somatic capacitance: print out the specific cap (cm) at the soma, as well as its area, then compute cap
// area: secname for (x,0) print x, area(x)
// also take into account the "well-clamped" computed capacitance: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3273682/ (taylor 2012, J comput neurosci)
// axial resistance:
// dendritic conductance:
// dendritic capacitance: print out the specific cap (cm) at various locations along the dendrite, as well as its area, then compute cap
// time constant: after reaching RMP, 0.15 nA negative current pulse of 100 ms duration. How long does it take to reach 63% of total steady state hyperpolarization? (Isokawa, 1997) http://www.sciencedirect.com/science/article/pii/S1385299X96000165
// .15 nA may be too large if Ih is involved, try <40 pA (http://www.jneurosci.org/content/27/46/12440.long)
// specific membrane resistance
// keep in mind space clamp issues: http://www.ncbi.nlm.nih.gov/pubmed/18184885
// issues to be aware of in measuring tau and cap: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2775376/
// how are these different properties related?
fcell.close()
// iClamp stuff
strdef tempString
f2 = new File()
sprint(pathstr,"%s%s%s%s%s", relpath, sl, "setupfiles", sl, "iclamp.dat")
f2.ropen(pathstr)
numiclamps = f2.scanvar
stimAMPvector = new Vector(numiclamps)
for i = 0, numiclamps-1 {
stimAMPvector.x(i) = f2.scanvar
}
f2.close()
for i = 0, numiclamps-1 {
for idxes=0, numCellidxes-1 {
cellind = mycellvecs.x[idxes]
gid=cellType[cellind].cellStartGid
if (pc.gid_exists(gid)) {
// add recording locations to cell:
strdef sname, cmdstr, outfile
cell = pc.gid2cell(gid)
if (cell.is_art==0) {
cell.soma { //forsec cell.all { // add recording location to middle of each section
sname = secname()
{index = strobj.substr(sname, ".")+1}
{strobj.right(sname, index)}
{index = strobj.substr(sname, "[")}
if (index>0) {strobj.left(sname, index)}
sprint(cmdstr, "{objref trace_%s}", cellType[cellind].cellType_string)
execute1(cmdstr)
sprint(cmdstr, "{trace_%s = new Vector(%g)}", cellType[cellind].cellType_string, (tstop-tstart)/dt)
execute1(cmdstr)
sprint(cmdstr, "{trace_%s.record(&cell.%s.v(0.5))}", cellType[cellind].cellType_string, sname)
execute1(cmdstr)
if (myi_flag==1) {
print_mechs(gid)
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "{objref %s_%s}", cellType[cellind].cellType_string, mechstring[rt].s)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s = new Vector(%g)}", cellType[cellind].cellType_string, mechstring[rt].s, (tstop-tstart)/dt)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s.record(&cell.%s.%s(0.5))}", cellType[cellind].cellType_string, mechstring[rt].s, sname, mechstring[rt].s)
execute1(cmdstr)
}
}
}
// add an iClamp (to be activated at various time points first)
sprint(cmdstr, "objref stim%g", cellind)
execute1(cmdstr)
sprint(cmdstr, "cell.soma stim%g = new IClamp(0.5)", cellind)
execute1(cmdstr)
sprint(cmdstr, "stim%g.dur = duration", cellind) //4000
execute1(cmdstr)
sprint(cmdstr, "stim%g.del = starttime", cellind) //0.01
execute1(cmdstr)
sprint(cmdstr, "stim%g.amp = stimAMPvector.x(i)", cellind) // can't be over 36.7 and must be greater than -22.85
execute1(cmdstr)
}
}
}
// run model
rrun() // Run the network simulation
//for cellind=0, numCellTypes-1 {
for idxes=0, numCellidxes-1 {
cellind = mycellvecs.x[idxes]
gid=cellType[cellind].cellStartGid
if (pc.gid_exists(gid)) {
// add recording locations to cell:
strdef sname, cmdstr, outfile
cell = pc.gid2cell(gid)
if (cell.is_art==0) {
// write out recordings
cell.soma { //forsec cell.all {
sname = secname()
index = strobj.substr(sname, ".")+1
strobj.right(sname, index)
index = strobj.substr(sname, "[")
if (index>0) {strobj.left(sname, index)}
sprint(outfile, "%s%s%s%strace_%s.%s(0.5).%g.dat", relpath, sl, resultspath, sl, cellType[cellind].cellType_string, sname, i)
f = new File(outfile)
f.wopen()
print_mechs(gid) // mechstring =
f.printf("t\tv")
if (myi_flag==1) {
for rt = 0, mechlength.x[0]-1 {
f.printf("\t%s", mechstring[rt].s)
}
}
f.printf("\n")
for j=0, (tstop-tstart)/dt-1 {
sprint(cmdstr, "{f.printf(\"%%g\\t%%g\", j*dt, trace_%s.x[j])}", cellType[cellind].cellType_string)
execute1(cmdstr)
if (myi_flag==1) {
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "{f.printf(\"\\t%%g\", %s_%s.x[j])}", cellType[cellind].cellType_string, mechstring[rt].s)
execute1(cmdstr)
}
}
{f.printf("\n")}
}
f.close()
}
sprint(cmdstr, "{objref stim%g}", cellind)
execute1(cmdstr)
}
}
}
} // end of iclamp
}
/***********************************************************************************************
V? PAIRED RECORDINGS
***********************************************************************************************/
objref f2, f3, syn
strdef tempString, prestr, poststr, spname
objref cell, nc[1], nil
objref ss, nobj
strdef pre_type, post_type, cmdstr1
objref synvec, fpair
if ((pair==2) || (pair==4)) {
print "2 or 4 pair = ", pair
f2 = new File()
{sprint(pathstr,"%s%ssetupfiles%sclamp%spairclamp.dat", relpath, sl, sl, sl)}
f2.ropen(pathstr)
current = f2.scanvar
holding = f2.scanvar
pretime = f2.scanvar
posttime = f2.scanvar
pairrevpot = f2.scanvar
pairrevflag = f2.scanvar
tstop = posttime // pretime + duration + posttime
f2.close
for i=0, numCellTypes-1 {
sprint(cmdstr, "%s_idx = %g", cellType[i].cellType_string, i)
{execute1(cmdstr)}
print cmdstr
}
f2 = new File()
{sprint(pathstr,"%s%ssetupfiles%sclamp%spairtypes.txt", relpath, sl, sl, sl)}
f2.ropen(pathstr)
numpairs = f2.scanvar
for r = 0, numpairs - 1 {
f2.scanstr(pre_type)
f2.scanstr(post_type)
sprint(cmdstr, "%s_idx", pre_type)
sprint(cmdstr1, "%s_idx", post_type)
if ((name_declared(cmdstr) > 0) && (name_declared(cmdstr1) > 0)) {
sprint(cmdstr, "precellType = %s_idx", pre_type)
{execute1(cmdstr)}
sprint(cmdstr, "postcellType = %s_idx", post_type)
{execute1(cmdstr)}
/*********** PRECELL ************/
precell.pp.number = 1
precell.pp.start = pretime
prestr=cellType[precellType].cellType_string
print "pre: ", prestr, " ", cellType[precellType].technicalType
/*********** POSTCELL ************/
jgid=cellType[postcellType].cellEndGid
poststr=cellType[postcellType].cellType_string
print "post: ", poststr, " ", cellType[postcellType].technicalType, " jgid: ", jgid
// add a voltage clamp and watch to the postcell at the soma
cell = pc.gid2cell(jgid)
sprint(sname,"o%s_%s", prestr, poststr)
sprint(cmdstr, "{objref vcl%s}", sname)
execute1(cmdstr)
if (pair==2) {
sprint(cmdstr, "{cell.soma vcl%s = new SEClamp(0.5)}", sname)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.dur1 = %g}", sname, tstop)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.amp1 = %g}", sname, holding)
execute1(cmdstr)
} else {
sprint(cmdstr, "{cell.soma vcl%s = new IClamp(0.5)}", sname)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.del = %g}", sname, 0)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.dur = %g}", sname, tstop)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.amp = %g}", sname, current)
execute1(cmdstr)
}
sprint(cmdstr, "{objref %s}", sname)
execute1(cmdstr)
sprint(cmdstr, "{%s = new Vector(%g)}", sname, (tstop-tstart)/dt)
execute1(cmdstr)
if (pair==2) {
sprint(cmdstr, "{%s.record(&vcl%s.i)}", sname, sname)
execute1(cmdstr)
} else {
sprint(cmdstr, "{%s.record(&cell.soma.v(0.5))}", sname)
execute1(cmdstr)
}
synWeight = cellType[precellType].wgtConns.x[postcellType]
numSyns = cellType[precellType].numSyns.x[postcellType]
conDelay = 0
synvec = new Vector(numSyns)
print prestr, " -> ", poststr, ": ", synWeight
if (synWeight>0) {
numSynTypes = cell.pre_list.o(precellType).count
sprint(outfile,"%s%s%s%s%s.%s.sids.dat", relpath, sl, resultspath, sl, prestr, poststr)
fpair = new File(outfile)
fpair.wopen()
fpair.printf("%g\n2\n", numeapair)
for p = 0, numeapair-1 {
objref nc[numSyns]
ransynlist.object(cell.randi).r.discunif(0,numSynTypes-1) // Create a uniform random INTEGER variable over the range specified (0 to # synapse types-1),
for q = 0, numSyns-1 {
synvec.x[q] = ransynlist.object(cell.randi).repick
if (pairrevflag==1) { //(strcmp(prestr,"ivycell")==0 || strcmp(prestr,"ngfcell")==0) {
if (strcmp(prestr,"ngfcell")==0) {
cell.pre_list.o(precellType).o(synvec.x[q]).ea = pairrevpot
cell.pre_list.o(precellType).o(synvec.x[q]).eb = pairrevpot
} else {
cell.pre_list.o(precellType).o(synvec.x[q]).e = pairrevpot
}
// cell.pre_list.o(precellType).o(synvec.x[q]).ea = pairrevpot
// cell.pre_list.o(precellType).o(synvec.x[q]).eb = -91 //pairrevpot
//} else {
// cell.pre_list.o(precellType).o(synvec.x[q]).e = pairrevpot
}
//print "p: ", p, " nw: ", cell.pre_list.o(precellType).o(synvec.x[q]).nw
nc[q] = pc.gid_connect(0, cell.pre_list.o(precellType).o(synvec.x[q])) // Connect the soma of the pre cell to the synapse on the post cell
nc[q].weight = synWeight // Set a synaptic connection weight
nc[q].delay = conDelay // Set a delay time (axonal delay)
}
rrun()
nc = nil
// print out watches
sprint(outfile,"%s%s%s%s%s.%s.%g.trace.dat", relpath, sl, resultspath, sl, prestr, poststr, p)
f = new File(outfile)
f.wopen()
f.printf("t\tpre\tpost\n")
for j=0, (tstop-tstart)/dt-1 {
sprint(cmdstr, "{f.printf(\"%%g\\t0\\t%%g\\n\", j*dt, %s.x[j])}", sname)
execute1(cmdstr)
}
f.close()
for q=0, numSyns-1 {
fpair.printf("%g\t%g\n", p, cell.pre_list.o(precellType).o(synvec.x[q]).sid)
}
}
fpair.close()
}
// clear postsyn clamp and watch
sprint(sname,"o%s_%s", prestr, poststr)
sprint(cmdstr, "%s = nil", sname)
execute1(cmdstr)
sprint(cmdstr, "vcl%s = nil", sname)
execute1(cmdstr)
}
}
f2.close
}
if ((pair==1) || (pair==3)) {
print "1 or 3 pair = ", pair, " numCellTypes = ", numCellTypes
f2 = new File()
{sprint(pathstr,"%s%ssetupfiles%sclamp%spairclamp.dat", relpath, sl, sl, sl)}
f2.ropen(pathstr)
current = f2.scanvar
holding = f2.scanvar
pretime = f2.scanvar
posttime = f2.scanvar
pairrevpot = f2.scanvar
pairrevflag = f2.scanvar
tstop = posttime // pretime + duration + posttime
f2.close
ConnData = conn
sprint(outfile, "%s%ssetupfiles%sload_cell_conns.hoc", relpath, sl, sl)
{load_file(outfile)} // Load in the cell connectivity info (weight)
for precellType= 0, numCellTypes - 1 {
igid=0
prestr=cellType[precellType].cellType_string
print prestr
precell.pp.number = 1
precell.pp.start = pretime
for postcellType= 0, numCellTypes - 1 {
synWeight = cellType[precellType].wgtConns.x[postcellType]
if (synWeight>0) {
jgid=cellType[postcellType].cellEndGid
poststr=cellType[postcellType].cellType_string
objref fpair
sprint(outfile,"%s%s%s%s%s.%s.sids.dat", relpath, sl, resultspath, sl, prestr, poststr)
fpair = new File(outfile)
fpair.wopen()
fpair.printf("%g\n2\n", numeapair)
// add a voltage clamp and watch to the postcell at the soma
cell = pc.gid2cell(jgid)
sprint(sname,"o%s_%s", prestr, poststr)
sprint(cmdstr, "{objref vcl%s}", sname)
execute1(cmdstr)
if (pair==1) {
sprint(cmdstr, "{cell.soma vcl%s = new SEClamp(0.5)}", sname)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.dur1 = %g}", sname, tstop)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.amp1 = %g}", sname, holding)
execute1(cmdstr)
} else {
sprint(cmdstr, "{cell.soma vcl%s = new IClamp(0.5)}", sname)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.del = %g}", sname, 0)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.dur = %g}", sname, tstop)
execute1(cmdstr)
sprint(cmdstr, "{vcl%s.amp = %g}", sname, current)
execute1(cmdstr)
}
for p = 0, numeapair-1 {
sprint(cmdstr, "{objref %s}", sname)
execute1(cmdstr)
sprint(cmdstr, "{%s = new Vector(%g)}", sname, (tstop-tstart)/dt)
execute1(cmdstr)
if (pair==1) {
sprint(cmdstr, "{%s.record(&vcl%s.i)}", sname, sname)
execute1(cmdstr)
} else {
sprint(cmdstr, "{%s.record(&cell.soma.v(0.5))}", sname)
execute1(cmdstr)
}
conDelay = 0
numSyns = cellType[precellType].numSyns.x[postcellType]
objref synvec, nc[numSyns]
synvec = new Vector(numSyns)
numSynTypes = cell.pre_list.o(precellType).count
objref nc[numSyns]
ransynlist.object(cell.randi).r.discunif(0,numSynTypes-1) // Create a uniform random INTEGER variable over the range specified (0 to # synapse types-1),
for q = 0, numSyns-1 {
synvec.x[q] = ransynlist.object(cell.randi).repick
if (pairrevflag==1) { //(strcmp(prestr,"ivycell")==0 || strcmp(prestr,"ngfcell")==0) {
if (strcmp(prestr,"ngfcell")==0) {
cell.pre_list.o(precellType).o(synvec.x[q]).ea = pairrevpot
cell.pre_list.o(precellType).o(synvec.x[q]).eb = pairrevpot
} else {
cell.pre_list.o(precellType).o(synvec.x[q]).e = pairrevpot
}
// cell.pre_list.o(precellType).o(synvec.x[q]).ea = pairrevpot
// cell.pre_list.o(precellType).o(synvec.x[q]).eb = -91 //pairrevpot
//} else {
// cell.pre_list.o(precellType).o(synvec.x[q]).e = pairrevpot
}
nc[q] = pc.gid_connect(0, cell.pre_list.o(precellType).o(synvec.x[q])) // Connect the soma of the pre cell to the synapse on the post cell
nc[q].weight = synWeight // Set a synaptic connection weight
nc[q].delay = conDelay // Set a delay time (axonal delay)
}
rrun()
nc = nil
for q=0, synvec.size()-1 {
fpair.printf("%g\t%g\n", p, cell.pre_list.o(precellType).o(synvec.x[q]).sid)
}
// print out watches
sprint(outfile,"%s%s%s%s%s.%s.%g.trace.dat", relpath, sl, resultspath, sl, prestr, poststr, p)
f = new File(outfile)
f.wopen()
f.printf("t\tpre\tpost\n")
for j=0, (tstop-tstart)/dt-1 {
sprint(cmdstr, "{f.printf(\"%%g\\t0\\t%%g\\n\", j*dt, %s.x[j])}", sname)
execute1(cmdstr)
}
f.close()
}
// clear postsyn clamp and watch
sprint(sname,"o%s_%s", prestr, poststr)
sprint(cmdstr, "%s = nil", sname)
execute1(cmdstr)
sprint(cmdstr, "vcl%s = nil", sname)
execute1(cmdstr)
fpair.close()
}
}
}
}
/***********************************************************************************************
V? MORPHOLOGY
***********************************************************************************************/
create xScale, yScale, zScale
proc anatscale() {
if ($4>0) { // if length arg is <= 0 then do nothing
xScale {
pt3dclear()
pt3dadd($1, $2, $3, 2)
pt3dadd($1+$4, $2, $3, 2)
}
yScale {
pt3dclear()
pt3dadd($1, $2, $3, 2)
pt3dadd($1, $2+$4, $3, 2)
}
zScale {
pt3dclear()
pt3dadd($1, $2, $3, 2)
pt3dadd($1, $2, $3+$4, 2)
}
}
}
objref showthese, ss, nobj
strdef pathstr
double vec[4]
showthese = new SectionList()
xScale {showthese.append()}
yScale {showthese.append()}
zScale {showthese.append()}
if (cellmorph==1) {
print "in morph!"
for idxes=0, numCellidxes-1 {
cellind = mycellvecs.x[idxes]
gid=cellType[cellind].cellStartGid
if (pc.gid_exists(gid)) {
// add recording locations to cell:
cell = pc.gid2cell(gid)
if (cell.is_art==0) {
showthese.append(cell.all)
ss = new Shape(showthese)
ss.color_list(showthese, 1)
ss.show(0)
ss.exec_menu("View = plot")
ss.size(&vec[0])
xpt = (vec[1]-vec[0])*.1+vec[0]
ypt = (vec[3]-vec[2])*.1+vec[2]
print "xpt: ", xpt, ", ypt: ", ypt
print "x: ", vec[0], " to ", vec[1], "; y: ", vec[2], " to ", vec[3]
anatscale(0,0,0,100) // xyz coords of origin, and length
//anatscale(vec[0]-10,vec[2]-10,0,100) // xyz coords of origin, and length
//ss.label(0.10,0.10,"100 um")
sprint(pathstr,"%s%s%s%s%s_morph.eps", relpath, sl, resultspath, sl, cellType[cellind].cellType_string)
ss.printfile(pathstr)
//ss = nobj
}
}
}
}
{pc.runworker()} // Everything after this line is executed only by the host node
// "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