/
terub_gpe_neuron.nestml
163 lines (133 loc) · 6.66 KB
/
terub_gpe_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
"""
terub_gpe - Terman Rubin neuron model
#####################################
Description
+++++++++++
terub_gpe is an implementation of a spiking neuron using the Terman Rubin model
based on the Hodgkin-Huxley formalism.
(1) **Post-syaptic currents:** Incoming spike events induce a post-synaptic change of current modelled
by an alpha function. The alpha function is normalised such that an event of
weight 1.0 results in a peak current of 1 pA.
(2) **Spike Detection:** 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.
References
++++++++++
.. [1] Terman, D. and Rubin, J.E. and Yew, A. C. and Wilson, C.J.
Activity Patterns in a Model for the Subthalamopallidal Network
of the Basal Ganglia
The Journal of Neuroscience, 22(7), 2963-2976 (2002)
.. [2] Rubin, J.E. and Terman, D.
High Frequency Stimulation of the Subthalamic Nucleus Eliminates
Pathological Thalamic Rhythmicity in a Computational Model
Journal of Computational Neuroscience, 16, 211-235 (2004)
"""
model terub_gpe_neuron:
state:
V_m mV = E_L # 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
gate_h real = 0.0 # gating variable h
gate_n real = 0.0 # gating variable n
gate_r real = 0.0 # gating variable r
Ca_con real = 0.0 # calcium concentration
equations:
# Parameters for Terman Rubin GPe Neuron
inline g_tau_n_0 ms = 0.05 ms
inline g_tau_n_1 ms = 0.27 ms
inline g_theta_n_tau mV = -40.0 mV
inline g_sigma_n_tau mV = -12.0 mV
inline g_tau_h_0 ms = 0.05 ms
inline g_tau_h_1 ms = 0.27 ms
inline g_theta_h_tau mV = -40.0 mV
inline g_sigma_h_tau mV = -12.0 mV
inline g_tau_r ms = 30.0 ms
# steady state values for gating variables
inline g_theta_a mV = -57.0 mV
inline g_sigma_a mV = 2.0 mV
inline g_theta_h mV = -58.0 mV
inline g_sigma_h mV = -12.0 mV
inline g_theta_m mV = -37.0 mV
inline g_sigma_m mV = 10.0 mV
inline g_theta_n mV = -50.0 mV
inline g_sigma_n mV = 14.0 mV
inline g_theta_r mV = -70.0 mV
inline g_sigma_r mV = -2.0 mV
inline g_theta_s mV = -35.0 mV
inline g_sigma_s mV = 2.0 mV
# time evolvement of gating variables
inline g_phi_h real = 0.05
inline g_phi_n real = 0.1 #Report: 0.1, Terman Rubin 2002: 0.05
inline g_phi_r real = 1.0
# Calcium concentration and afterhyperpolarization current
inline g_epsilon 1/ms = 0.0001 /ms
inline g_k_Ca real = 15.0 #Report:15, Terman Rubin 2002: 20.0
inline g_k1 real = 30.0
inline I_exc_mod real = -convolve(g_exc, exc_spikes) * nS * V_m
inline I_inh_mod real = convolve(g_inh, inh_spikes) * nS * (V_m-E_gg)
inline tau_n real = g_tau_n_0 + g_tau_n_1 / (1. + exp(-(V_m-g_theta_n_tau)/g_sigma_n_tau))
inline tau_h real = g_tau_h_0 + g_tau_h_1 / (1. + exp(-(V_m-g_theta_h_tau)/g_sigma_h_tau))
inline tau_r real = g_tau_r
inline a_inf real = 1. / (1. + exp(-(V_m-g_theta_a)/g_sigma_a))
inline h_inf real = 1. / (1. + exp(-(V_m-g_theta_h)/g_sigma_h))
inline m_inf real = 1. / (1. + exp(-(V_m-g_theta_m)/g_sigma_m))
inline n_inf real = 1. / (1. + exp(-(V_m-g_theta_n)/g_sigma_n))
inline r_inf real = 1. / (1. + exp(-(V_m-g_theta_r)/g_sigma_r))
inline s_inf real = 1. / (1. + exp(-(V_m-g_theta_s)/g_sigma_s))
inline I_Na real = g_Na * m_inf * m_inf * m_inf * gate_h * (V_m - E_Na)
inline I_K real = g_K * gate_n * gate_n * gate_n * gate_n * (V_m - E_K )
inline I_L real = g_L * (V_m - E_L )
inline I_T real = g_T * a_inf* a_inf * a_inf * gate_r * (V_m - E_Ca)
inline I_Ca real = g_Ca * s_inf * s_inf * (V_m - E_Ca)
inline I_ahp real = g_ahp * (Ca_con / (Ca_con + g_k1)) * (V_m - E_K )
# V dot -- synaptic input are currents, inhib current is negative
V_m' = ( -(I_Na + I_K + I_L + I_T + I_Ca + I_ahp) * pA + I_e + I_stim + I_exc_mod * pA + I_inh_mod * pA) / C_m
# channel dynamics
gate_h' = g_phi_h * ((h_inf - gate_h) / tau_h) / ms # h-variable
gate_n' = g_phi_n * ((n_inf - gate_n) / tau_n) / ms # n-variable
gate_r' = g_phi_r * ((r_inf - gate_r) / tau_r) / ms # r-variable
# Calcium concentration
Ca_con' = g_epsilon*(-I_Ca - I_T - g_k_Ca * Ca_con)
# synapses: alpha functions
## alpha function for the g_inh
kernel g_inh = (e/tau_syn_inh) * t * exp(-t/tau_syn_inh)
## alpha function for the g_exc
kernel g_exc = (e/tau_syn_exc) * t * exp(-t/tau_syn_exc)
parameters:
E_L mV = -55 mV # Resting membrane potential
g_L nS = 0.1 nS # Leak conductance
C_m pF = 1 pF # Capacitance of the membrane
E_Na mV = 55 mV # Sodium reversal potential
g_Na nS = 120 nS # Sodium peak conductance
E_K mV = -80.0 mV # Potassium reversal potential
g_K nS = 30.0 nS # Potassium peak conductance
E_Ca mV = 120 mV # Calcium reversal potential
g_Ca nS = 0.15 nS # Calcium peak conductance
g_T nS = 0.5 nS # T-type Calcium channel peak conductance
g_ahp nS = 30 nS # Afterpolarization current peak conductance
tau_syn_exc ms = 1 ms # Rise time of the excitatory synaptic alpha function
tau_syn_inh ms = 12.5 ms # Rise time of the inhibitory synaptic alpha function
E_gg mV = -100 mV # Reversal potential for inhibitory input (from GPe)
refr_T ms = 2 ms # Duration of refractory period
# constant external input current
I_e pA = 0 pA
input:
exc_spikes <- excitatory spike
inh_spikes <- inhibitory spike
I_stim pA <- continuous
output:
spike
update:
V_m_old = V_m
if is_refractory:
refr_t -= resolution()
integrate_odes()
onCondition(not is_refractory and V_m > 0 mV 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