// Threshold spacetests using a binary search method.
//
// Program flow:
// Takes location range, Gaussian width for time distribution, and location increments.
// Also takes type of run: AMPA only, NMDA only, or both.
// Records number of synapses at which a single action potential is fired.
// Also may record voltage at soma (later)
// Updated to take a starting seed and to work with a cell template instead of the old ball and stick.
// Updated to use variable AMPA conductance.
strdef filename
proc ThreshSpace() { local AMPAc,UPPERLIM,resetHere,locRange,gaussTime,repStat,incrBy,found,branchLength,toggle,matInd,jvar,stSeed localobj search_ind,finMat,floovec,cell
// Runs the whole test.
// Inputs:
// $1: locRange is the range to be uniformly sampled in space
// $2: gaussTime is the width of the Gaussian to be sampled for timing
// $3: repStat is the number of repetitions for each location, to collect statistics in variation
// $4: incrBy is the distance in microns each trial is separated by.
// $5: branchLength is the number of spines on the tested branch.
// $6: toggle is the kind of synapse: 0 BOTH, 1 AMPA, 2 NMDA
// $s7: filename.
// $8: the starting seed
// $o9: the cell itself
// $10: AMPA conductance, with normal .0005 umho
UPPERLIM = 10000 // If the number of synapses goes above this, the trial moves on.
locRange = $1
gaussTime = $2
repStat = $3
incrBy = $4
branchLength = $5
toggle = $6
filename = $s7
stSeed = $8
cell = $o9
AMPAc = $10
// The binary search vector holding indices. Order: lower bound, current test number, upper bound.
search_ind = resetVec(200)
floovec = new Vector(1,((branchLength-locRange)/incrBy)+1)
// For the (test) case where incrBy is larger than locRange
if (incrBy >= locRange) {
floovec = new Vector(1,(branchLength/incrBy))
}
floovec.floor()
finMat = new Matrix(floovec.x[0]-1,repStat)
for r_ind = 0,repStat-1 {
seed = (stSeed+r_ind)*(r_ind+stSeed)
matInd = -1
for (l_ind=0;l_ind<(branchLength-locRange);l_ind=l_ind+incrBy) {
// For each run, search_ind starts from the previous threshold level
found = 0
matInd = matInd+1
if (search_ind.x[1] == -1) {
resetHere = UPPERLIM-1
} else {
resetHere = search_ind.x[1]
}
search_ind = resetVec(resetHere)
while (found != 1) {
// found variable is 1 when AP is fired at synnum n+1 but not at synnum n
found = testSyns(search_ind.x[1],locRange,gaussTime,toggle,seed,l_ind,cell,AMPAc)
// Sees if this synapse number works
print search_ind.x[1]
// Setting an upper limit for # synapses tested
if (search_ind.x[1] > UPPERLIM) {
found = 1
search_ind.x[1] = -1
}
search_ind = binSearch(search_ind,found)
// If synapse number doesn't work, changes the search index and upper limit
}
finMat.x[matInd][r_ind] = search_ind.x[1] + 1
// Saves final synapse number to the matrix
jvar = printf("Iter %d, loc %d: %d\n",r_ind,l_ind,search_ind.x[1])
}
}
saveSpace(finMat,filename)
}
/*------------------ MAIN FUNCTIONS ----------------------*/
func testSyns() { local AMPAc, search_ind, locRange, gaussTime, toggle, seed, l_ind, found localobj r2, tlist, tlist2, apc
// If an AP fires at search_ind and search_ind+1, returns 2.
// If no AP fires at either, returns 0.
// If AP fired at search_ind+1 and not search_ind (it's at threshold), returns 1.
//
// Inputs:
// $1: search_ind is how many synapses to use
// $2: locRange is the testing range
// $3: gaussTime is the width of the Gaussian for the time distribution
// $4: toggle is the usual synapse type ID
// $5: seed is for the norm dist for time and the unif dist for space.
// $6: l_ind is where the test is starting.
// $o7: the cell
// $8: AMPA conductance (normal .0005 umho)
search_ind = $1
locRange = $2
gaussTime = $3
toggle = $4
seed = $5
l_ind = $6
AMPAc = $8
r2 = new Vector(2,0)
// Runs for search_ind and search_ind+1
for ts_ind = 0,1 {
tlist = new List()
tlist2 = new List()
tlist = setSyns(search_ind+ts_ind,locRange,gaussTime,toggle,seed,l_ind,$o7,AMPAc)
if (toggle==0) {
tlist2 = setSyns(search_ind+ts_ind,locRange,gaussTime,2,seed,l_ind,$o7,AMPAc)
}
// Puts an APCount object at the axon
$o7.axon {
apc = new APCount(0.5)
apc.thresh = 0
apc.n = 0
}
//Takes tlist to adjust tstop.
runTest(tlist)
// Sees if an AP fired
if (apc.n > 0) {
if (ts_ind==0) {
return 2
} else {
r2.x[ts_ind] = 1
}
}
}
return r2.sum()
}
obfunc binSearch() {
// Adjusts the searchVec based on results from testSyns.
//
// Inputs:
// $o1: search_ind is a vector; the first element is the number of synapses being tested
// $2: found tells what to do with the search vector. 0 is none, 1 is threshold, 2 is both.
// If too low
if ($2 == 0) {
$o1.x[0] = $o1.x[1]
if ($o1.x[2] == $o1.x[1]) {
$o1.x[1] = $o1.x[1]*2
$o1.x[2] = $o1.x[1]
} else {
$o1.x[1] = (($o1.x[2]-$o1.x[1])/2)+$o1.x[1]
$o1.floor()
}
}
// If too high
if ($2 == 2) {
$o1.x[2] = $o1.x[1]
$o1.x[1] = (($o1.x[1]-$o1.x[0])/2)+$o1.x[0]
$o1.floor()
}
return $o1
}
proc saveSpace() { localobj mat, f
// Takes a matrix and saves it.
//
// Inputs:
// $o1: the matrix to be saved (rows are location, columns are trial)
// $s2: filename with the full directory
mat = $o1
f = new File()
f.wopen($s2)
mat.fprint(f, " %g")
f.close()
}
obfunc resetVec() { localobj vec
// Gives the index vector given some general starting point. It's a 3-element vector where the
// first element is the lower bound, middle is current index, last is upper bound. The function
// starts the upper bound at the current index value, since the search pattern works by doubling.
//
// Input:
// $1: the number of synapses to start out with.
vec = new Vector(3,0)
vec.x[1] = $1
vec.x[2] = $1
return vec
}
/*------------------ SUB FUNCTIONS -----------------------*/
// In testSyns()
obfunc setSyns() { local s_ind localobj tlist,r,normr,svec,tvec
// Creates a list of placed and timed synapses. Does not deal with BOTH case (creates AMPA for
// 0 and 1 toggle and NMDA for 2) because that's taken care of in testSyns().
//
// Inputs:
// $1: synapse number.
// $2: locRange, how widely synapses are distributed.
// $3: gaussTime, the width of the Gaussian used in sampling onset times.
// $4: toggle is explained above; slightly different for this function since it does one at a time.
// $5: seed
// $6: l_ind is the spine this test is starting on.
// $o7: the cell
// $8: AMPA conductance (normal .0005 umho), used only when toggle is 0 or 1
// Output:
// A single list of one kind of synapse; preset locations and firing times.
tlist = makeSyns($4,$1)
r = new Random($5)
r.uniform($6,$6+$2)
normr = new Random($5*10)
normr.normal(0,$3)
svec = new Vector($1)
tvec = new Vector($1)
svec.setrand(r)
svec.floor()
tvec.setrand(normr)
tvec.add((-1)*tvec.min()+2)
for s_ind = 0, $1-1 {
$o7.spine_head[svec.x[s_ind]] {
tlist.o(s_ind).loc(.5)
tlist.o(s_ind).onset() = tvec.x[s_ind]
if ($4!=2) {
tlist.o(s_ind).gmax() = $8
}
}
}
return tlist
}
// In testSyns()
proc runTest() { local i localobj vec
// Initializes cell and runs to tstop, which is 100 + last firing time.
// Input:
// $o1: tlist.
vec = giveTimes($o1)
init()
tstop = vec.max() + 100
while (t<tstop) {
fadvance()
}
}
// In setSyns()
obfunc makeSyns() { local i localobj tlist
// Takes kind of synapses and number and creates a list of them. 0 and 1 make AMPA, 2 makes NMDA.
//
// Inputs:
// $1: toggle
// $2: number of synapses.
// Output:
// tlist
tlist = new List()
if ($1!=2) {
for i=0, ($2-1) tlist.append(new syn_g())
} else {
for i=0, ($2-1) tlist.append(new nmda())
}
return tlist
}
// In runTest() and setSyns()
obfunc giveTimes() { local i localobj vec
// Takes synapse list and makes a vector of all onset times.
// Input:
// $o1 is the synapse list.
vec = new Vector()
for i = 0,$o1.count()-1 {
vec.append($o1.o(i).onset())
}
return vec
}
