-
Notifications
You must be signed in to change notification settings - Fork 1
/
init_exp.hoc
180 lines (162 loc) · 6.94 KB
/
init_exp.hoc
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
//
// Entry point for the model. There are many options and switches here, some defined in
// included hoc files. See startup_*.sh scripts in the same directory for examples of
// how to run this file. Most parameters are specified on the command-line.
//
VERBOSE=0
DEBUG=0
if (name_declared("verbose")==5) {
if (verbose >= 1) {
printf("*** verbose enabled!\n")
VERBOSE = 1
}
if (verbose >= 2) {
printf("*** debug enabled!\n")
DEBUG = 1
}
}
/*
* Before anything, figure out if we're running Richy (Cell 1) or
* Starfish (Cell 2).
*/
if (name_declared("cell") != 5) {
printf("*** FATAL ERROR: no cell type! ABORTING...\n")
quit()
}
if (cell != 1 && cell != 2) {
printf("*** FATAL ERROR: cell type must be 1 or 2! ABORTING...\n")
quit()
}
load_file("nrngui.hoc")
if (cell == 1) {
load_file("initactivesdRichy.hoc")
} else {
load_file("initactivesdStarfish.hoc")
}
load_file("rn.hoc")
load_file("iprotocol.hoc")
load_file("params.hoc")
load_file("nsegfuncs.hoc")
load_file("param_active_cond.hoc")
load_file("active_mods.hoc") // for changing gh on the fly across model
load_file("ihold_funcs.hoc") // for refitting ihold after changing gh (above)
load_file("dendIh_funcs.hoc")
v_init = -74
cvode_active(1)
cvode.atol(0.005)
print "DC input resistance of cell is ", rn(), "megohms"
// For some reason this needs to be after calling rn()
geom_nseg() // need this as Cm changed from loading fitted ses file
// Define maximum conductance densities as per command-line args
err = 0
if (set_gbars()) { err = 1 }
if (err) {
printf("*** FATAL ERROR. ABORTING... \n")
quit()
}
// Fit holding current, only if ihold not specified as commandline argument
tstop=1000 // only run for initial period
objref recv_hold
if (name_declared("ihold") == 5) {
holdclamp.amp = ihold
} else {
injclamp_del = injclamp.del // probably 180, see in iprotocol.hoc
injclamp.del=1000 // so that we can do the fitting properly, then set it back afterwards
recv_hold = new Vector()
recv_hold.record(&soma[0].v(0.5)) // first need to find spike, so use variable dt from cvode
target_Vm = -60 - 13.8 // == -73.8mV; NB: -13.8mV for JP correction
//leeway = 2 // leeway in fitting holding current; start off at 2mV
leeway = 5 // leeway in fitting holding current; start off at 5mV
//target_error = 0.05 // mV
//target_error = 0.1 // mV - changed it as 0.05 wouldn't always converge
target_error = 0.15 // mV - changed it as 0.1 wouldn't always converge
// start from -15pA, we will then either go more positive (likely), or more negative, as needed
start_amp = -0.015
holdclamp.amp = start_amp // holdclamp defined in iprotocol.hoc
stopspike_amp = start_amp // amount of amp to stop spike with
amp_increment = 0.001 // amount to increment amp, at first
found_spike = 0
found_hold = 0
while (!found_hold) {
printf("*** FITTING holding clamp amp: trying %g nA\n", holdclamp.amp)
run()
/*
Two phases: 1) Find amp to stop spiking at, using variable dt; 2) Find amp
to get avg_Vm to target_Vm, using fixed dt. We need both phases because of
Vector.record() functioning. If we just have one phase, and set variable dt,
then we won't get accurate avg_Vm in the last fifth of initial period as there
will be too few data points, so won't converge properly. If we have fixed dt,
then likely won't catch the spike in the case where the spike threshold is
lower than target_Vm. In this case, the fitting procedure will run indefinitely
because the spiking will throw off avg_Vm.
So our solution are these two phases. The first is to definitively find
the spike threshold (rather, amount of holding current needed to stop spiking).
Then we start the second phase with fixed dt in recv_hold.record() vector. This
is where we look for the holding current to keep avg_Vm close to target_Vm.
*But* if in the course of this second phase we notice that our holding clamp
amp is getting to be smaller (i.e., higher nA value) than the current needed to
stop the spike, then we know that this is one of the model cells where spiking
occurs even before we've reached target_Vm from below. In this case we stop
with the holding clamp amp set to stopspike_amp.
*/
if (!found_spike) {
if (recv_hold.max() >= -10) {
stopspike_amp = holdclamp.amp - amp_increment
printf("*** GENERATED SPIKE!!! FOUND holding clamp amp to stop spike at %g nA\n", stopspike_amp)
found_spike = 1
} else if (holdclamp.amp >= 0.015) {
printf("*** CELL NOT INTRINSICALLY SPIKING; reached upper bound on check. Continuing with holding clamp fitting...\n")
stopspike_amp = holdclamp.amp
found_spike = 1
} else {
holdclamp.amp = holdclamp.amp + amp_increment
}
if (found_spike) {
// now re-init vector and use fixed dt in order to get good avg_Vm calculation
recv_hold = new Vector()
recv_hold.record(&soma[0].v(0.5),20) // this will give us 50 timepoints in recv_hold
holdclamp.amp = start_amp // reset holdclamp.amp too
}
} else {
if (holdclamp.amp > stopspike_amp) {
holdclamp.amp = holdclamp.amp - amp_increment // just to output what we actually tested for
if (want_output) {
write_err_and_quit()
} else {
printf("*** FATAL ERROR: holding current fitting algorithm found inappropriate model. ABORTING...\n")
quit()
}
} else {
// Get average of somatic Vm in last 200ms of initial period, so last 5th of 1s-long trace
end = recv_hold.size() - 1
start = (end) * 4/5 + 1
avg_Vm = recv_hold.sum(start, end) / (end - start + 1)
printf("*** AVG Vm %g mV\n", avg_Vm)
if (avg_Vm >= target_Vm - leeway && avg_Vm <= target_Vm + leeway) {
if (abs(abs(avg_Vm) - abs(target_Vm)) <= target_error) {
// we found our holding current amp!
printf("*** FOUND holding clamp amp at %g nA\n", holdclamp.amp)
found_hold = 1
} else {
leeway = leeway * 2/3
amp_increment = amp_increment / 2
printf("*** GOING TO FINER resolution of amp_increment: %g and leeway: %g avg_Vm: %g target_Vm: %g\n", amp_increment, leeway, avg_Vm, target_Vm)
}
} else if (avg_Vm < target_Vm - leeway) {
// still hyperpolarized wrt target Vm, so inject *less* hyperpolarizing current
holdclamp.amp = holdclamp.amp + amp_increment
} else if (avg_Vm > target_Vm - leeway) {
// we are depolarized wrt target Vm, so inject *more* hyperpolarizing current
holdclamp.amp = holdclamp.amp - amp_increment
} else {
// umm, something broke
printf("*** FATAL ERROR: impossible control flow in holding current code; avg_Vm %g holdclamp.amp %g", avg_Vm, holdclamp.amp)
quit()
}
}
}
}
injclamp.del = injclamp_del
}
tstop=1300
run()