1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
// $Id: initbatpar.hoc,v 1.7 2011/06/10 01:46:17 ted Exp $
// Batch implementation suitable for serial or parallel execution
// for stimulus current i iterating over a range of values
// run a simulation and report spike frequency
// save i & corresponding f to a file
// optionally plot fi curve
// wrapping code in { } prevents printing undesired 0s etc. to terminal
{ load_file("stdgui.hoc") } // load standard library but don't open NMM toolbar
objref pc // create ParallelContext instance here
pc = new ParallelContext() // so it exists whenever we need it
///// Simulation parameters
TSTOP = 500 // ms, more than long enough for 15 spikes at ISI = 25 ms
AMP0 = 0.1 // nA--minimum stimulus
D_AMP = 0.02 // nA--stim increment between runs
NRUNS = 30
///// Model specification
{ load_file("cell.hoc") }
///// Instrumentation
// experimental manipulations
objref stim
soma stim = new IClamp(0.5)
stim.del = 1 // ms
stim.dur = 1e9
stim.amp = 0.1 // nA
// $1 is run #
// returns corresponding stimulus amplitude
func fstimamp() {
return AMP0 + $1*D_AMP
}
// sets stimulus amplitude to appropriate value for this run
proc setparams() {
stim.amp = fstimamp($1)
}
// data recording and analysis
objref nc, spvec, nil // to record spike times
// count only those spikes that get to distal end of dend
dend nc = new NetCon(&v(1), nil)
nc.threshold = -10 // mV
spvec = new Vector()
{ nc.record(spvec) }
NSETTLE = 5 // ignore the first NSETTLE ISI (allow freq to stablize)
NINVL = 10 // # ISI from which frequency will be calculated
NMIN = NSETTLE + NINVL // ignore recordings with fewer than this # of ISIs
freq = 0
proc postproc() { local nspikes, t1, t2
freq = 0
nspikes = spvec.size()
if (nspikes > NMIN) {
t2 = spvec.x(nspikes-1) // last spike
t1 = spvec.x(nspikes-1-NINVL) // NINVL prior to last spike
freq = NINVL*1e3/(t2-t1) // t1 and t2 are separated by NINVL ISIs
}
}
///// Simulation control
tstop = TSTOP
func fi() { // set params, execute a simulation, analyze and return results
setparams($1) // set parameters for this run
run()
postproc() // analyze the data
return freq
}
// batch control
trun = startsw()
{ pc.runworker() } // start execute loop on the workers
// code beyond this point is executed only by the master
// the master must now post jobs to the bulletin board
objref svec, fvec
proc batchrun() { local ii, tmp
svec = new Vector()
fvec = new Vector()
for ii = 0, $1-1 pc.submit("fi", ii) // post all jobs to bulletin board
// retrieve results from bulletin board
while (pc.working) { // if a result is ready, get it; if not, pick a job and do it
fvec.append(pc.retval()) // get frequency
pc.unpack(&tmp)
svec.append(tmp) // get job number
printf(".") // indicate progress
}
}
batchrun(NRUNS)
{ pc.done() }
// all simulations have been completed
// and the workers have been released
// but the boss still has things to do
///// Reporting of results
fvec = fvec.index(svec.sortindex()) // rearrange fvec according to svec sortindex
{
svec.sort()
// but svec contains job numbers, not actual stimulus amplitudes
svec.apply("fstimamp") // convert job numbers to stimulus amplitudes
}
// now svec holds the current amplitudes in increasing order
// and fvec holds the corresponding firing frequencies
proc saveresults() { local ii localobj file
file = new File($s1)
file.wopen()
file.printf("label:%s\n%d\n", "f", fvec.size)
for ii=0, fvec.size-1 {
file.printf("%g\t%g\n", svec.x[ii], fvec.x[ii])
}
file.close()
}
saveresults("fi.dat")
print " "
print startsw() - trun, "seconds" // report run time
quit()
|
