/
wb_cond_multisyn_neuron.nestml
190 lines (147 loc) · 7.56 KB
/
wb_cond_multisyn_neuron.nestml
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
"""
wb_cond_multisyn - Wang-Buzsaki model with multiple synapses
############################################################
Description
+++++++++++
wb_cond_multisyn is an implementation of a modified Hodkin-Huxley model.
Spike detection is done by a combined threshold-and-local-maximum search: if
there is a local maximum above a certain threshold of the membrane potential,
it is considered a spike.
AMPA, NMDA, GABA_A, and GABA_B conductance-based synapses with
beta-function (difference of two exponentials) time course.
References
++++++++++
.. [1] Wang, X.J. and Buzsaki, G., (1996) Gamma oscillation by synaptic
inhibition in a hippocampal interneuronal network model. Journal of
Neuroscience, 16(20), pp.6402-6413.
See also
++++++++
wb_cond_multisyn
"""
model wb_cond_multisyn_neuron:
state:
V_m mV = -65. mV # Membrane potential
V_m_old mV = E_L # Membrane potential at previous timestep for threshold check
refr_t ms = 0 ms # Refractory period timer
is_refractory boolean = false
Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init ) # Inactivation variable h for Na
Act_n real = alpha_n_init / (alpha_n_init + beta_n_init) # Activation variable n for K
g_AMPA real = 0
g_NMDA real = 0
g_GABAA real = 0
g_GABAB real = 0
g_AMPA$ real = AMPAInitialValue
g_NMDA$ real = NMDAInitialValue
g_GABAA$ real = GABA_AInitialValue
g_GABAB$ real = GABA_BInitialValue
equations:
recordable inline I_syn_ampa pA = -convolve(g_AMPA, AMPA) * nS * ( V_m - AMPA_E_rev )
recordable inline I_syn_nmda pA = -convolve(g_NMDA, NMDA) * nS * ( V_m - NMDA_E_rev ) / ( 1 + exp( ( NMDA_Vact - V_m ) / NMDA_Sact ) )
recordable inline I_syn_gaba_a pA = -convolve(g_GABAA, GABA_A) * nS * ( V_m - GABA_A_E_rev )
recordable inline I_syn_gaba_b pA = -convolve(g_GABAB, GABA_B) * nS * ( V_m - GABA_B_E_rev )
recordable inline I_syn pA = I_syn_ampa + I_syn_nmda + I_syn_gaba_a + I_syn_gaba_b
inline I_Na pA = g_Na * Act_m_inf(V_m)**3 * Inact_h * ( V_m - E_Na )
inline I_K pA = g_K * Act_n**4 * ( V_m - E_K )
inline I_L pA = g_L * ( V_m - E_L )
Inact_h' = ( alpha_h(V_m) * ( 1 - Inact_h ) - beta_h(V_m) * Inact_h ) / ms # h-variable
Act_n' = ( alpha_n(V_m) * ( 1 - Act_n ) - beta_n(V_m) * Act_n ) / ms # n-variable
V_m' = ( -( I_Na + I_K + I_L ) + I_e + I_stim + I_syn ) / C_m
#############
# Synapses
#############
kernel g_AMPA' = g_AMPA$ - g_AMPA / AMPA_Tau_2,
g_AMPA$' = -g_AMPA$ / AMPA_Tau_1
kernel g_NMDA' = g_NMDA$ - g_NMDA / NMDA_Tau_2,
g_NMDA$' = -g_NMDA$ / NMDA_Tau_1
kernel g_GABAA' = g_GABAA$ - g_GABAA / GABA_A_Tau_2,
g_GABAA$' = -g_GABAA$ / GABA_A_Tau_1
kernel g_GABAB' = g_GABAB$ - g_GABAB / GABA_B_Tau_2,
g_GABAB$' = -g_GABAB$ / GABA_B_Tau_1
parameters:
g_Na nS = 3500.0 nS # Sodium peak conductance
g_K nS = 900.0 nS # Potassium peak conductance
g_L nS = 10 nS # Leak conductance
C_m pF = 100.0 pF # Membrane Capacitance
E_Na mV = 55.0 mV # Sodium reversal potential
E_K mV = -90.0 mV # Potassium reversal potentia
E_L mV = -65.0 mV # Leak reversal Potential (aka resting potential)
V_Tr mV = -55.0 mV # Spike Threshold
refr_T ms = 2 ms # Duration of refractory period
# Parameters for synapse of type AMPA, GABA_A, GABA_B and NMDA
AMPA_g_peak nS = 0.1 nS # peak conductance
AMPA_E_rev mV = 0.0 mV # reversal potential
AMPA_Tau_1 ms = 0.5 ms # rise time
AMPA_Tau_2 ms = 2.4 ms # decay time, Tau_1 < Tau_2
NMDA_g_peak nS = 0.075 nS # peak conductance
NMDA_Tau_1 ms = 4.0 ms # rise time
NMDA_Tau_2 ms = 40.0 ms # decay time, Tau_1 < Tau_2
NMDA_E_rev mV = 0.0 mV # reversal potential
NMDA_Vact mV = -58.0 mV # inactive for V << Vact, inflection of sigmoid
NMDA_Sact mV = 2.5 mV # scale of inactivation
GABA_A_g_peak nS = 0.33 nS # peak conductance
GABA_A_Tau_1 ms = 1.0 ms # rise time
GABA_A_Tau_2 ms = 7.0 ms # decay time, Tau_1 < Tau_2
GABA_A_E_rev mV = -70.0 mV # reversal potential
GABA_B_g_peak nS = 0.0132 nS # peak conductance
GABA_B_Tau_1 ms = 60.0 ms # rise time
GABA_B_Tau_2 ms = 200.0 ms # decay time, Tau_1 < Tau_2
GABA_B_E_rev mV = -90.0 mV # reversal potential for intrinsic current
# constant external input current
I_e pA = 0 pA
internals:
AMPAInitialValue real = compute_synapse_constant( AMPA_Tau_1, AMPA_Tau_2, AMPA_g_peak )
NMDAInitialValue real = compute_synapse_constant( NMDA_Tau_1, NMDA_Tau_2, NMDA_g_peak )
GABA_AInitialValue real = compute_synapse_constant( GABA_A_Tau_1, GABA_A_Tau_2, GABA_A_g_peak )
GABA_BInitialValue real = compute_synapse_constant( GABA_B_Tau_1, GABA_B_Tau_2, GABA_B_g_peak )
alpha_n_init real = -0.05 * (V_m / mV + 34.0) / (exp(-0.1 * (V_m / mV + 34.0)) - 1.0)
beta_n_init real = 0.625 * exp(-(V_m / mV + 44.0) / 80.0)
alpha_m_init real = 0.1 * (V_m / mV + 35.0) / (1.0 - exp(-0.1 * (V_m / mV + 35.)))
beta_m_init real = 4.0 * exp(-(V_m / mV + 60.0) / 18.0)
alpha_h_init real = 0.35 * exp(-(V_m / mV + 58.0) / 20.0)
beta_h_init real = 5.0 / (exp(-0.1 * (V_m / mV + 28.0)) + 1.0)
input:
AMPA <- spike
NMDA <- spike
GABA_A <- spike
GABA_B <- spike
I_stim pA <- continuous
output:
spike
update:
if is_refractory:
refr_t -= resolution()
V_m_old = V_m
integrate_odes()
onCondition(not is_refractory and V_m > V_Tr and V_m_old > V_m):
# threshold crossing and maximum
refr_t = refr_T # start of the refractory period
is_refractory = true
emit_spike()
onCondition(is_refractory and refr_t <= resolution() / 2):
# end of refractory period
refr_t = 0 ms
is_refractory = false
function compute_synapse_constant(Tau_1 ms, Tau_2 ms, g_peak nS) real:
# Factor used to account for the missing 1/((1/Tau_2)-(1/Tau_1)) term
# in the ht_neuron_dynamics integration of the synapse terms.
# See: Exact digital simulation of time-invariant linear systems
# with applications to neuronal modeling, Rotter and Diesmann,
# section 3.1.2.
exact_integration_adjustment real = ( ( 1 / Tau_2 ) - ( 1 / Tau_1 ) ) * ms
t_peak ms = ( Tau_2 * Tau_1 ) * ln( Tau_2 / Tau_1 ) / ( Tau_2 - Tau_1 )
normalisation_factor real = 1 / ( exp( -t_peak / Tau_1 ) - exp( -t_peak / Tau_2 ) )
return (g_peak / nS) * normalisation_factor * exact_integration_adjustment
function Act_m_inf(V_m mV) real:
return alpha_m(V_m) / ( alpha_m(V_m) + beta_m(V_m) )
function alpha_m(V_m mV) real:
return 0.1 * (V_m / mV + 35.0) / (1.0 - exp(-0.1 * (V_m / mV + 35.)))
function beta_m(V_m mV) real:
return 4.0 * exp(-(V_m / mV + 60.0) / 18.0)
function alpha_n(V_m mV) real:
return -0.05 * (V_m / mV + 34.0) / (exp(-0.1 * (V_m / mV + 34.0)) - 1.0)
function beta_n(V_m mV) real:
return 0.625 * exp(-(V_m / mV + 44.0) / 80.0)
function alpha_h(V_m mV) real:
return 0.35 * exp(-(V_m / mV + 58.0) / 20.0)
function beta_h(V_m mV) real:
return 5.0 / (exp(-0.1 * (V_m / mV + 28.0)) + 1.0)