-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathSD.ode
450 lines (383 loc) · 17 KB
/
SD.ode
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
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
########################################################################################
# This code is archived with “The Role of Glutamate in Neuronal Ion Homeostasis: A Case Study of Spreading# Depolarization” by Hubel et al. PLoS Computational Biology (2017).
# The software used is AUTO.
# Feel free to use/modify this code for your own research. Please, cite our paper if you use this code.
########################################################################################
# RATE EQUATIONS FOR SD WITH GLUTAMATE: #
# 'v' membrane potential in [mV] #
# 'n/h_HH' Hodgkin-Huxley gating variables in [1] #
# 'r_NMDA/AMPA' NMDA/AMPA gating variables in [1] #
# 'nki/ncli' amount of K/Cl in ICS in [fmol] #
# 'dnk_b/g' amount of K that is exchanged with external bath/glia in [fmol] #
# 'ngl_c/e' amount of glutamate in all clefts/ECS in [fmol] #
########################################################################################
v' = 1000. * v_DOT
n_HH' = 1000. * n_HH_DOT
h_HH' = 1000. * h_HH_DOT
r_AMPA' = 1000. * r_AMPA_DOT * f_receptor
r_NMDA' = 1000. * r_NMDA_DOT * f_receptor
nki' = 1000. * nki_DOT
ncli' = 1000. * ncli_DOT
dnk_b' = 1000. * dnk_b_DOT
dnk_g' = 1000. * dnk_g_DOT
ngl_c' = 1000. * ngl_c_DOT
ngl_e' = 1000. * ngl_e_DOT
init v=-71.1
init n_HH=0.07143807
init h_HH=0.9774849
init r_AMPA=0.
init r_NMDA=0.
init nki=1050.
init ncli=67.5
init dnk_b=0.
init dnk_g=0.
init ngl_c=0.
init ngl_e=0.
k_rep = 0.001
ngl_rel' = 1000. * (J_rel - k_rep / ngl_tot * (ngl_c_i + ngl_e_i + ngl_c_g + ngl_e_g))
ngl_c_i' = 1000. * (V_c_i - k_rep / ngl_tot * ngl_c_i)
ngl_c_g' = 1000. * (V_c_g - k_rep / ngl_tot * ngl_c_g)
ngl_c_e' = 1000. * Jglu_diff
ngl_e_i' = 1000. * (V_e_i - k_rep / ngl_tot * ngl_e_i)
ngl_e_g' = 1000. * (V_e_g - k_rep / ngl_tot * ngl_e_g)
init ngl_rel=0
init ngl_c_i=0
init ngl_c_g=0
init ngl_c_e=0
init ngl_e_i=0
init ngl_e_g=0
#############################################################################
# GLUTAMATE UPTAKE/RELEASE PARAMETERS: #
# 'exp_rel' exponent for 'v'-dependence of glutamate release in [1] #
# 'Vmax' maximal velocity of glutamate uptake by neuron in [mM/msec] #
#############################################################################
par f_upt=1
par f_co=1
par Vmax=0.03
exp_rel = 2
################################
# switch for receptor dynamics #
################################
f_receptor = 1
#######################################
# SD from extrac. potassium elevation #
#######################################
f_OGD = 1
par KB=15.
###############
# SD from OGD #
###############
# par KB=4.
# par delta_t=15
# par t_OGD=40
#
# t1 = 20
# t2 = t1 + delta_t
# t3 = t2 + t_OGD
# t4 = t3 + delta_t
# f01_OGD = 1
# f12_OGD = heav(t - t1) * (-(t - t1) / (t2 - t1))
# f23_OGD = heav(t - t2) * (t - t2) / (t2 - t1)
# f34_OGD = heav(t - t3) * (t - t3) / (t4 - t3)
# f4x_OGD = heav(t - t4) * (-(t - t4) / (t4 - t3))
# f_OGD = f01_OGD + f12_OGD + f23_OGD + f34_OGD + f4x_OGD
###############################################
# GATING FUNCTIONS: #
# Hodgkin-Huxley (HH) and NMDA/AMPA receptors #
# (gating time constants in [msec]) #
###############################################
an_HH = 0.01 * (v + 34.0) / (1.0 - exp(-0.1 * (v + 34.0)))
bn_HH = 0.125 * exp(-(v + 44.0)/80.0)
am_HH = 0.1 * (v + 30.0) / (1.0 - exp(-0.1 * (v + 30.0)))
bm_HH = 4.0 * exp(-(v + 55.0)/18.0)
ah_HH = 0.07 * exp(-(v + 44.0)/20.0)
bh_HH = 1.0 / (1.0 + exp(-0.1 * (v + 14.0)))
m_HH = am_HH / (am_HH + bm_HH)
ar_AMPA = 1.1
br_AMPA = 0.19
ar_NMDA = 0.072
br_NMDA = 0.0066
#############################################################################################
# MORPHOLOGY: #
# 'vl_i0/e0/g0' initial/resting volume of ICS, ECS and glia in [um^3] #
# 'A_m' membrane surface area in [um^2] (same for neuron and glia) #
# 'h','r' height and radius of synaptic cleft in [um] #
# 'vl_en' volume inside the glial envelope in [um^3] #
# 'A_sig' cross section area for glutamate diffusion from cleft into ECS in [um^2] #
#############################################################################################
vl_i0 = 7500.
vl_e0 = 2500.
vl_g0 = 7500.
vl_tot = vl_i0 + vl_e0 + vl_g0
A_m = 18000.
r = 100e-3
h = 20e-3
vl_en = 6 * pi * r**2 * h
A_sig = 4 * pi * r**2 * 0.05
##################################################################################
# ION CONENTRATIONS: #
# 'ioni/e0' initial/resting conc. in ICS/ECS in [mM=mmol/l] #
# 'xi' conc. of impermeant uncharged matter in ICS in [mM] #
# 'nioni/e0' initial/resting ion amount in ICS/ECS in [fmol] #
# 'nani/e' amount of impermeant anions in ICS/ECS in [fmol] #
# 'nxi/e' amount of impermeant uncharged particles in ICS/ECS in [fmol] #
# 'ni0/e0/g0' initial amount of all particles in ICS/ECS/glia in [fmol] #
# 'dnna_b/g' amount of Na that is exchanged with external bath/glia in [fmol] #
# 'dncl_b/g' amount of Cl that is exchanged with external bath/glia in [fmol] #
# 'NNAI/E' current amount of Na in ICS/ECS in [fmol] #
# 'NKE/CLE' current amount of K and Cl in ECS in [fmol] #
# 'ni/e/g' current amount of all particles in ICS/ECS/glia in [fmol] #
# 'n_tot' total amount of particles in the whole system in [fmol] #
# 'vl_i/e/g' current volume of ICS/ECS/glia in [um^3] #
# 'IONI/E' current ion conc. in ICS/ECS in [mM=mmol/l] #
# (conversion factor 1e3 from [fmol/um^3] to [mM] for 'KE') #
##################################################################################
nai0 = 15.
nae0 = 144.
ki0 = 140.
ke0 = 4.
cli0 = 9.
cle0 = 130.
xi = 5.
nki0 = ki0 * vl_i0 * 1e-3
nke0 = ke0 * vl_e0 * 1e-3
nnai0 = nai0 * vl_i0 * 1e-3
nnae0 = nae0 * vl_e0 * 1e-3
ncli0 = cli0 * vl_i0 * 1e-3
ncle0 = cle0 * vl_e0 * 1e-3
nani = nki0 + nnai0 - ncli0
nane = nke0 + nnae0 - ncle0
nxi = xi * vl_i0 * 1e-3
nxe = vl_e0/vl_i0 * (nxi + nani + nki0 + nnai0 + ncli0) - (nane + nke0 + nnae0 + ncle0)
ni0 = nnai0 + nki0 + ncli0 + nani + nxi
ne0 = nnae0 + nke0 + ncle0 + nane + nxe
ng0 = ni0 / vl_i0 * vl_g0
dncl_g = 0.8 * dnk_g
dnna_g =-0.2 * dnk_g
NNAI = nnai0 + nki0 - nki - ncli0 + ncli
NNAE = nnae0 + nnai0 - NNAI - dnna_g
NKE = nke0 + nki0 - nki - dnk_g - dnk_b
NCLE = ncle0 + ncli0 - ncli - dncl_g
ni = NNAI + nki + ncli + nani + nxi
ne = NNAE + NKE + NCLE + nane + nxe
ng = ng0 + dnna_g + dnk_g + dncl_g
n_tot = ni + ne + ng
vl_i = ni / n_tot * vl_tot
vl_e = ne / n_tot * vl_tot
vl_g = ng / n_tot * vl_tot
NAI = NNAI / vl_i * 1e3
KI = NKI / vl_i * 1e3
CLI = NCLI / vl_i * 1e3
NAE = NNAE / vl_e * 1e3
KE = NKE / vl_e * 1e3
CLE = NCLE / vl_e * 1e3
#######################################################
# NERNST POTENTIALS: #
# 'EK/NA/CL' Nernst potetnials for K/Na/Cl in [mV] #
#######################################################
EK = 26.64 * log(KE / KI)
ENA = 26.64 * log(NAE / NAI)
ECL =-26.64 * log(CLE / CLI)
#######################################################################################
# CURRENTS ON NEURAL MEMBRANE: #
# 'gk/na_g' maximal conductances of gated HH-like K/Na channels in [mS/cm^2] #
# 'gk/na/cl_l' conductances of K/Na/Cl leak channels in [mS/cm^2] #
# 'nsyn_AP/SD' number of synapses activated in action potential (AP)/SD in [1] #
# 'g_NMDA/AMPA_AP' effective NMDA/AMPA receptor conductance in single AP in [mS] #
# 'g_NMDA/AMPA' effective NMDA/AMPA conductance density for SD in [mS/cm^2] #
# 'max_p' maximal pump rate for K/Na-exchange in [uA/cm^2] #
# 'mg' Mg concentration in ECS in [mM] #
# 'f_NMDA' factor for voltage-dependent Mg block in NMDA channels in [1] #
# 'IK/NA_g' gated K/Na current in [uA/cm^2] #
# 'IK/NA/CL_l' leak K/Na/Cl current in [uA/cm^2] #
# 'IK/NA_NMDA' K/Na current through NMDA receptor gates in [uA/cm^2] #
# 'IK/NA_AMPA' K/Na current through AMPA receptor gates in [uA/cm^2] #
# 'IP' K/Na-exchange pump current in [uA/cm^2] #
# 'INA/K' sum of all Na/K currents in [uA/cm^2] #
#######################################################################################
gk_g = 40.
gna_g = 100.
gna_l = 0.0135
gk_l = 0.05
gcl_l = 0.05
nsyn_AP = 20
nsyn_SD = 5000
par g_NMDA_AP=1e-7
par g_AMPA_AP=3.5e-7
g_AMPA = g_AMPA_AP / (A_m * 1e-8) * nsyn_SD / nsyn_AP
g_NMDA = g_NMDA_AP / (A_m * 1e-8) * nsyn_SD / nsyn_AP
mg = 1.2
f_NMDA = 1.0 / (1.0 + 0.33 * mg * exp(-(0.07*v + 0.7)))
max_p = 6.46
INA_g = gna_g * m_HH**3 * h_HH * (v-ENA)
IK_g = gk_g * n_HH**4 * (v-EK)
INA_l = gna_l * (v-ENA)
IK_l = gk_l * (v-EK)
ICL_l = gcl_l * (v-ECL)
IK_NMDA = g_NMDA * r_NMDA * (v-EK) * f_NMDA
INA_NMDA = g_NMDA * r_NMDA * (v-ENA) * f_NMDA
IK_AMPA = g_AMPA * r_AMPA * (v-EK)
INA_AMPA = g_AMPA * r_AMPA * (v-ENA)
IP = max_p / (1.0 + exp((15 - nai)/3.)) / (1. + exp(5.5 - ke)) * f_OGD
INA = INA_l + INA_g + INA_AMPA + INA_NMDA + 3. * IP
IK = IK_l + IK_g + IK_AMPA + IK_NMDA - 2. * IP
#############################################################################################
# GLUTAMATE RELEASE IN ONE SYNAPSE #
# 'ngl_tot/i' total/intracellular amount of glutamate in [fmol] #
# 'v_cr' critical mebrane potential for glutamate release in [mV] #
# 'v_hi' highest value of membrane potential seen in the model in [mV] #
# 'rel_max' maximal (i.e. for 'v=v_hi') release rate in [fmol/msec] #
# 'f_rel' factor for 'V'-dependence of glutamate release in [1] #
# 'J_rel' glutamate release flux into all synaptic clefts involved in SD [fmol/msec] #
#############################################################################################
rel_max = 1.5e-5
par ngl_tot=10
ngl_avl = ngl_tot - ngl_rel
v_cr =-50
v_hi = 50
f_rel = heav(v-v_cr) * ((v-v_cr)/(v_hi-v_cr))**exp_rel
J_rel = nsyn_SD * rel_max * f_rel * ngl_avl / ngl_tot
####################################################################################
# GLUTAMATE UPTAKE: #
# 'gl_c/e' glutamate concentration in clefts/ECS in [mM] #
# 'k_m' affinity pf uptake system in [mM] #
# 'Vmax_c/e_i' maximal uptake velocity from cleft/ECS into ICS in [mM/msec] #
# 'Vmax_c/e_g' maximal uptake velocity from cleft/ECS into glia in [mM/msec] #
# 'F' Faraday's constant in [A*sec/mol] #
# 'conv' current-to-flux conversion s.th. 'conv/vl_i*INA' is in [mM/msec] #
# 'nsyn_AP/SD' number of synapses activated in action potential (AP)/SD in [1] #
# 'V_c/e_i' uptake velocity from cleft/ECS into ICS in [fmol/msec] #
# 'V_c/e_g' uptake velocity from cleft/ECS into glia in [fmol/msec] #
# 'JNa/Cl/K_co' glutamate cotransport flux of Na/Cl/K into ICS in [mM/msec] #
# 'INa/Cl/K_co' cotransport Na/Cl/K currents on neural membrane in [uA/cm^2] #
####################################################################################
gl_c = ngl_c / (vl_en * nsyn_SD) * 1e3
gl_e = ngl_e / vl_e * 1e3
k_m = 3e-2
Vmax_c_i = f_upt * Vmax
Vmax_e_i = f_upt * Vmax * 0.12 * vl_e0 / vl_e
Vmax_c_g = f_upt * Vmax * 4
Vmax_e_g = f_upt * Vmax * 4 * 0.24 * vl_e0 / vl_e
F = 96485
conv = A_m / F * 10
V_c_i = Vmax_c_i * gl_c / (gl_c + k_m) * vl_en * nsyn_SD * 1e-3 * f_OGD
V_e_i = Vmax_e_i * gl_e / (gl_e + k_m) * vl_e * 1e-3 * f_OGD
V_c_g = Vmax_c_g * gl_c / (gl_c + k_m) * vl_en * nsyn_SD * 1e-3 * f_OGD
V_e_g = Vmax_e_g * gl_e / (gl_e + k_m) * vl_e * 1e-3 * f_OGD
JNa_co = (V_c_i + V_e_i) * 3 * f_co
JCl_co = (V_c_i + V_e_i) * f_co
JK_co =-(V_c_i + V_e_i) * f_co
INa_co = JNa_co / conv * 1e3
ICl_co =-JCl_co / conv * 1e3
IK_co = JK_co / conv * 1e3
################################################################################
# DIFFUSION PROCESSES: #
# 'lambda' strength of bath coupling in [1/msec] #
# 'kbath' level of K concentration of bath in [mM] #
# 'JK_diff' diffusive K flux between ECS and bath in [mM/msec] #
# 'D_glu' glutamate diffusion constatn in [um^2/msec] #
# 'dx' distance between cleft and staionary bath concentraion in [um] #
# 'Jglu_diff' diffusive glutamate flux from cleft to ECS in [fmol/msec] #
################################################################################
par lambda=0.0001
JK_diff =-lambda * (KB - KE) * f_OGD
D_glu = 0.3
dx = 20
Jglu_diff =-D_glu * nsyn_SD / dx * A_sig * (ngl_e/vl_e - ngl_c/(vl_en*nsyn_SD))
###################################################################
# GLIAL BUFFERING: #
# 'k1' buffering constant in [mM/msec] #
# 'k2' backward buffering rate in [mM/msec] #
# 'dnk_max' maximal amount of K that can go into glia in [fmol] #
# 'JK_glia' K uptake flux in [fmol/msec] #
###################################################################
k1 = 1.44e-2
k2 = 5.1e-3
dnk_max = 350
JK_glia =-(k2 - k1 / (1.0 + exp((5.5-KE)/2.5)) * (dnk_max - dnk_g)/dnk_max) * vl_e0 * 1e-3 * f_OGD
##########################################################
# RATE FUNCTIONS USED BY SOLVER: #
# 'C_m' capacitance of neural membrane in [uF/cm^2] #
# 'phi' conventional time scale parameter in [1/msec] #
##########################################################
C_m = 1
phi = 3
v_DOT =-1 / C_m * (INA + IK + ICL_l + INA_co + IK_co + ICL_co)
n_HH_DOT = phi * (an_HH * (1-n_HH) - bn_HH * n_HH)
h_HH_DOT = phi * (ah_HH * (1-h_HH) - bh_HH * h_HH)
r_AMPA_dot = gl_c * ar_AMPA * (1.0 - r_AMPA) - br_AMPA * r_AMPA
r_NMDA_dot = gl_c * ar_NMDA * (1.0 - r_NMDA) - br_NMDA * r_NMDA
nki_DOT =-conv * 1e-3 * (IK + IK_co)
ncli_DOT = conv * 1e-3 * (ICL_l + ICl_co)
dnk_b_DOT = JK_diff * vl_e * 1e-3
dnk_g_DOT = JK_glia
ngl_c_DOT = J_rel - Jglu_diff - V_c_i - V_c_g
ngl_e_DOT = Jglu_diff - V_e_i - V_e_g
####################################
# AUXILIARY VARIABLES FOR PLOTTING #
####################################
aux _vl_i = vl_i
aux _vl_e = vl_e
aux _vl_g = vl_g
aux _gl_c = gl_c
aux _gl_e = gl_e
aux _EK = EK
aux _ENA = ENA
aux _ECL = ECL
aux _KI = ki
aux _KE = KE
aux _NAI = NAI
aux _NAE = NAE
aux _CLI = cli
aux _CLE = CLE
##############################
# NUMERICS AND PLOT SETTINGS #
##############################
@ maxstor=10000000, bounds=10000000
@ bell=0
# @ meth=Runge-Kutta
# @ dt=5e-5
#@ meth=stiff
#@ dt=1e-4
@ meth=stiff
@ dt=1e-3
#@ total=800
#@ xhi=800
@ total=400
@ xhi=400
#@ ylo=-300, yhi=300
#@ nplot=6
#@ yp1=v
#@ yp2=_ECL
#@ yp3=_CLI
#@ yp4=_CLE
#@ yp5=dnk_b
#@ yp6=dnk_g
@ ylo=-150, yhi=160
@ nplot=9
@ yp1=v
@ yp2=_EK
@ yp3=_ENA
@ yp4=_ECL
@ yp5=_vl_i
@ yp6=_vl_e
@ yp7=_vl_g
@ yp8=_gl_c
@ yp9=_gl_e
#@ ylo=-0.5, yhi=200
#@ nplot=8
#@ yp1=_gl_e
#@ yp2=_gl_c
#@ yp3=_KI
#@ yp4=_KE
#@ yp5=_NAI
#@ yp6=_NAE
#@ yp7=_CLI
#@ yp8=_CLE
#@ ylo=0, yhi=10000
#@ nplot=3
#@ yp1=_vl_i
#@ yp2=_vl_e
#@ yp3=_vl_g
done