-
Notifications
You must be signed in to change notification settings - Fork 2
/
integrator.h
354 lines (329 loc) · 10.3 KB
/
integrator.h
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
//general parameters
#define PI 3.1415926
#define NCELL 8 //number STN, GPe, or GPi cells
#define STNGPE NCELL*2 //sum of STN and GPe cells
#define NEQN 6 //number of ODEs per cell
#define NTCCELL 2 //number of TC cells
#define NTCEQN 3 //number of ODEs per TC cell
#define P_cm 1.0
//STN parameters
#define P_phin 0.75
#define P_phih 0.75
#define P_phica 0.75 //0.75 in Terman's code, not available from his paper.
#define P_eps 3.75e-5 //5.0e-5 in Terman's code
#define P_kca 22.5
#define P_phir 0.2
#define P_alpha 1.0
#define P_thetag 30.0
#define P_beta 0.05
#define P_gGS 0.9
#define P_vGS -100.0
#define P_gl 2.25
#define P_vl -60.0
#define P_gk 45.0
#define P_vk -80.0
#define P_gna 37.5
#define P_vna 55.0
#define P_gt 0.5
#define P_vca 140.0
#define P_gca 0.5
#define P_gahp 9.0
#define P_k1 15.0
#define P_thetam -30.0
#define P_sigmam 15.0
#define P_thetaa -63.0
#define P_sigmaa 7.8
#define P_thetas -39.0
#define P_sigmas 8.0
#define P_thetan -32.0
#define P_sigman 8.0
#define P_thetah -39.0
#define P_sigmah -3.1
#define P_thetar -67.0
#define P_sigmar -2.0
#define P_thetab 0.4
#define P_sigmab -0.1
#define P_thetaH -39.0
#define P_sigmaH 8.0
#define P_taun0 1.0
#define P_taun1 100.0
#define P_thn -80.0
#define P_sigmant -26.0
#define P_tauh0 1.0
#define P_tauh1 500.0
#define P_thh -57.0
#define P_sigmaht -3.0
#define P_taur0 40.0
#define P_taur1 17.5
#define P_thr 68.0
#define P_sigmart -2.2
#define ISTN 25.0
//GPe parameters
#define P_phing 0.1
#define P_phihg 0.05
#define P_phirg 1.0
#define P_alphag 1.0
#define P_thetagg 20.0
#define P_betag 0.1
#define P_epsg 1.0e-4
#define P_kcag 15.0
#define P_gSG 0.6 //0.3 in Terman 2004
#define P_vSG 0.0
#define P_gGG 0.0005 //normal = 2.0, PD = 0.05 or 0.0005
#define P_vGG -80.0
#define P_glg 0.1
#define P_vlg -55.0
#define P_gkg 30.0
#define P_vkg -80.0
#define P_gnag 120.0
#define P_vnag 55.0
#define P_gtg 0.5
#define P_vcag 120.0
#define P_gcag 0.15
#define P_gahpg 30.0
#define P_k1g 30.0
#define P_thetamg -37.0
#define P_sigmamg 10.0
#define P_thetaag -57.0
#define P_sigmaag 2.0
#define P_thetasg -35.0
#define P_sigmasg 2.0
#define P_thetang -50.0
#define P_sigmang 14.0
#define P_thetahg -58.0
#define P_sigmahg -12.0
#define P_thetarg -70.0
#define P_sigmarg -2.0
#define P_thetaHg -57.0
#define P_sigmaHg 2.0
#define P_taun0g 0.05
#define P_taun1g 0.27
#define P_thng -40.0
#define P_sigmantg -12.0
#define P_tauh0g 0.05
#define P_tauh1g 0.27
#define P_thhg -40.0
#define P_sigmahtg -12.0
#define P_taurg 30.0
#define IGPE -13.3 //normal = -2.0, PD = -13.0
//GPi parameters
#define P_phini 0.1
#define P_phihi 0.05
#define P_phiri 1.0
#define P_alphai 2.0
#define P_thetagi 20.0
#define P_betai 0.08
#define P_epsi 1.0e-4
#define P_kcai 15.0
#define P_gSI 0.2 //0.3 in Terman 2004
#define P_vSI 0.0
#define P_gGI 1.7 //1.0 in Terman 2004
#define P_vGI -100.0
#define P_gli 0.1
#define P_vli -55.0
#define P_gki 30.0
#define P_vki -80.0
#define P_gnai 120.0
#define P_vnai 55.0
#define P_gti 0.5
#define P_vcai 120.0
#define P_gcai 0.15
#define P_gahpi 30.0
#define P_k1i 30.0
#define P_thetami -37.0
#define P_sigmami 10.0
#define P_thetaai -57.0
#define P_sigmaai 2.0
#define P_thetasi -35.0
#define P_sigmasi 2.0
#define P_thetani -50.0
#define P_sigmani 14.0
#define P_thetahi -58.0
#define P_sigmahi -12.0
#define P_thetari -70.0
#define P_sigmari -2.0
#define P_thetaHi -57.0
#define P_sigmaHi 2.0
#define P_taun0i 0.05
#define P_taun1i 0.27
#define P_thni -40.0
#define P_sigmanti -12.0
#define P_tauh0i 0.05
#define P_tauh1i 0.27
#define P_thhi -40.0
#define P_sigmahti -12.0
#define P_tauri 30.0
#define IGPI 3.0
//TC parameters (obtained from Terman 2004)
#define P_gIT 0.08 //0.06 in Terman 2004
#define P_vIT -85.0
#define P_glt 0.05
#define P_vlt -70.0
#define P_gkt 5.0
#define P_vkt -90.0
#define P_gnat 3.0
#define P_vnat 50.0
#define P_gtt 5.0
#define P_vcat 0.0
#define P_thetamt -37.0
#define P_sigmamt 7.0
#define P_thetaat -60.0
#define P_sigmaat 6.2
#define P_thetaht -41.0
#define P_sigmahtt -4.0
#define P_thetart -84.0
#define P_sigmartt -4.0
#define ISM 10.0 //5.0 in Terman 2004, 7.5-10 may be good
#define PSM 50.0 //25.0 in Terman 2004, 50 in his figures
#define DSM 5.0
#define SHIFT 2.5 //phase shift when counting spike number
#define THLD_TC -30.0 // threshold for TC voltage spikes
//integrator parameters
#define NEQ (NCELL*NEQN*3+NTCCELL*NTCEQN) //number of ODEs in total
#define LRW ( 22 + NEQ * ( NEQ + 9 ) )
#define LIW ( 20 + NEQ )
#define ITOL 2
#define RTOL 1.0e-10
#define ATOL 1.0e-12
#define ATOLV 1.0e-10 //absolute tolerance for v
#define ITASK 1
#define ISTATE 1
#define IOPT 0
#define JT 2
#define TOUT 0.1 //first output time
#define TSTEP 0.1 //output step size
#define NSTEP 65536
#define NDELAY 25000 //No. of steps before the control is applied
//current input parameters
#define ILOW 1.0 //lowest current pulse distance
#define ISTEP 5.0 //step size for increasing pulse distance
#define DEV_DBS 0.1 //white noise of I_{DBS}
#define DEV_NET 0.01 //heterogeneity in the network itself (0 to 1)
//parameters for counting spike time
#define NSPK 2000 //maximum no. of Gpi spikes
#define VTHLD 0.5 //Gpi synaptic output threshold for counting spikes
//parameters for calculating auto- and cross-correlations
#define NSTART 10000 //when to start calculating correlations (in steps)
#define NEND 55000 //when to end calculating correlations (in steps)
#define DURATION 8192 //duration for correlation calculation (in steps)
#define BIN_RESLN 100 //binning resolution of Gpi spikes (in steps)
#define NBIN DURATION/BIN_RESLN+1 //maximum no. of Gpi bins
#define COMBI NCELL-1 //no. of comparisons for crosscorrelation
//parameters for convolution
#define DURATION2 DURATION*2
#define MSTEP 21 //width of the response function (in bins)
#define WIDTH 1 //Gaussian width (in bins)
//parameters for cost function calculation
#define W1 200.0 //1.0 //weight for autocorrelation term
#define W2 100.0 //weight for crosscorrelation term
#define W3 0.0002 //weight for Gpi spike number in first bin term
#define W4 0 //0.1 //current cost weight
//function declaration
int integrator( int [], double [] ); //my integrator
extern "C" void lsoda_( void( int *, double *, double *, double * ), int *,
double *, double *, double *, int *, double *, double *, int *, int *, int *,
double *, int *, int *, int *, void( int *, double *, double *, int *, int *,
double *, int * ), int * ); //ODE integrator
void fex( int *, double *, double *, double * ); //ODE container
void jdum( int *, double *, double *, int *, int *, double *, int * ); //dummy
//STN currents
double Il( const double v );
double Ik( const double v, const double n );
double Ina( const double v, const double h );
double It( const double v, const double r );
double Ica( const double v );
double Iahp( const double v, const double ca );
//STN functions
double minf( const double v );
double ainf( const double v );
double sinf( const double v );
double ninf( const double v );
double hinf( const double v );
double rinf( const double v );
double binf( const double r );
double Hinf( const double v );
double taun( const double v );
double tauh( const double v );
double taur( const double v );
//GPe currents
double Ilg( const double v );
double Ikg( const double v, const double n );
double Inag( const double v, const double h );
double Itg( const double v, const double r );
double Icag( const double v );
double Iahpg( const double v, const double ca );
//GPe functions
double minfg( const double v );
double ainfg( const double v );
double sinfg( const double v );
double ninfg( const double v );
double hinfg( const double v );
double rinfg( const double v );
double Hinfg( const double v );
double taung( const double v );
double tauhg( const double v );
//GPi currents
double Ili( const double v );
double Iki( const double v, const double n );
double Inai( const double v, const double h );
double Iti( const double v, const double r );
double Icai( const double v );
double Iahpi( const double v, const double ca );
//GPi functions
double minfi( const double v );
double ainfi( const double v );
double sinfi( const double v );
double ninfi( const double v );
double hinfi( const double v );
double rinfi( const double v );
double Hinfi( const double v );
double tauni( const double v );
double tauhi( const double v );
//TC currents
double Ilt( const double v );
double Ikt( const double v, const double h );
double Inat( const double v, const double h );
double Itt( const double v, const double r );
double taurt( const double v );
double tauht( const double v );
//TC functions
double minft( const double v );
double ainft( const double v );
double hinft( const double v );
double rinft( const double v );
double Heaviside( const double var );
//other functions
//inhibition of GPe to STN
double IsynGS(const double, int, const double [ ]);
//neighboring inhibition of GPe to GPe
double IsynGG(const double, int, const double [ ] );
//synaptic current from GPi to TC
double IsynIT1( const double [], const double );
double IsynIT2( const double [], const double );
//DBS current to STN
double Idbs( const unsigned long int, const int );
//cost value calculation
//Ism current to TC
double Ism( const double );
//Ism current to TC with delay, used in peak counting
double Ism2( const double );
//heterogeneity / random number generator
double hetero( const int, const double );
//random square pulse DBS input generated from specified PDF
void randInput( const double [], double [] );
//random number generator for specified probability distribution function
double randGen2( const double [], const int, const double, const double );
//auto- and cross-correlation calculation
double costCal( const double [][NSTEP] );
//count the number of spikes
int spkCount ( const double [], double [] );
//calculate autocorrelation
void autoCal( const double [], int [] );
//convolution calculation
void convCal( const int [], float [] );
//calculate standard deviation
double stdCal ( const float[] );
void four1( float [], unsigned long int, int ); //fft
void realft( float [], unsigned long, int );
void twofft( float [], float [], float [], float [], unsigned long );
void convlv( float [], unsigned long, float [], unsigned long, int, float [] );