/
terub_stn_neuron.nestml
170 lines (136 loc) · 6.92 KB
/
terub_stn_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
"""
terub_stn - Terman Rubin neuron model
#####################################
Description
+++++++++++
terub_stn 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_stn_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 STN Neuron
#time constants for slow gating variables
inline tau_n_0 ms = 1.0 ms
inline tau_n_1 ms = 100.0 ms
inline theta_n_tau mV = -80.0 mV
inline sigma_n_tau mV = -26.0 mV
inline tau_h_0 ms = 1.0 ms
inline tau_h_1 ms = 500.0 ms
inline theta_h_tau mV = -57.0 mV
inline sigma_h_tau mV = -3.0 mV
inline tau_r_0 ms = 7.1 ms # Guo 7.1 Terman02 40.0
inline tau_r_1 ms = 17.5 ms
inline theta_r_tau mV = 68.0 mV
inline sigma_r_tau mV = -2.2 mV
#steady state values for gating variables
inline theta_a mV = -63.0 mV
inline sigma_a mV = 7.8 mV
inline theta_h mV = -39.0 mV
inline sigma_h mV = -3.1 mV
inline theta_m mV = -30.0 mV
inline sigma_m mV = 15.0 mV
inline theta_n mV = -32.0 mV
inline sigma_n mV = 8.0 mV
inline theta_r mV = -67.0 mV
inline sigma_r mV = -2.0 mV
inline theta_s mV = -39.0 mV
inline sigma_s mV = 8.0 mV
inline theta_b real = 0.25 # Guo 0.25 Terman02 0.4
inline sigma_b real = 0.07 # Guo 0.07 Terman02 -0.1
#time evolvement of gating variables
inline phi_h real = 0.75
inline phi_n real = 0.75
inline phi_r real = 0.5 # Guo 0.5 Terman02 0.2
# Calcium concentration and afterhyperpolarization current
inline epsilon 1/ms = 0.00005 / ms # 1/ms Guo 0.00005 Terman02 0.0000375
inline k_Ca real = 22.5
inline k1 real = 15.0
inline I_exc_mod pA = -convolve(g_exc, exc_spikes) * nS * V_m
inline I_inh_mod pA = convolve(g_inh, inh_spikes) * nS * (V_m - E_gs)
inline tau_n ms = tau_n_0 + tau_n_1 / (1. + exp(-(V_m-theta_n_tau)/sigma_n_tau))
inline tau_h ms = tau_h_0 + tau_h_1 / (1. + exp(-(V_m-theta_h_tau)/sigma_h_tau))
inline tau_r ms = tau_r_0 + tau_r_1 / (1. + exp(-(V_m-theta_r_tau)/sigma_r_tau))
inline a_inf real = 1. / (1. +exp(-(V_m-theta_a)/sigma_a))
inline h_inf real = 1. / (1. + exp(-(V_m-theta_h)/sigma_h));
inline m_inf real = 1. / (1. + exp(-(V_m-theta_m)/sigma_m))
inline n_inf real = 1. / (1. + exp(-(V_m-theta_n)/sigma_n))
inline r_inf real = 1. / (1. + exp(-(V_m-theta_r)/sigma_r))
inline s_inf real = 1. / (1. + exp(-(V_m-theta_s)/sigma_s))
inline b_inf real = 1. / (1. + exp((gate_r-theta_b)/sigma_b)) - 1. / (1. + exp(-theta_b/sigma_b))
inline I_Na pA = g_Na * m_inf * m_inf * m_inf * gate_h * (V_m - E_Na)
inline I_K pA = g_K * gate_n * gate_n * gate_n * gate_n * (V_m - E_K )
inline I_L pA = g_L * (V_m - E_L )
inline I_T pA = g_T *a_inf*a_inf*a_inf*b_inf*b_inf* (V_m - E_Ca)
inline I_Ca pA = g_Ca * s_inf * s_inf * (V_m - E_Ca)
inline I_ahp pA = g_ahp * (Ca_con / (Ca_con + 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) + I_e + I_stim + I_exc_mod + I_inh_mod) / C_m
#channel dynamics
gate_h' = phi_h *((h_inf - gate_h) / tau_h) # h-variable
gate_n' = phi_n *((n_inf - gate_n) / tau_n) # n-variable
gate_r' = phi_r *((r_inf - gate_r) / tau_r) # r-variable
# Calcium concentration
Ca_con' = epsilon*( (-I_Ca - I_T ) / pA - 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 = -60 mV # Resting membrane potential
g_L nS = 2.25 nS # Leak conductance
C_m pF = 1 pF # Capacity of the membrane
E_Na mV = 55 mV # Sodium reversal potential
g_Na nS = 37.5 nS # Sodium peak conductance
E_K mV = -80 mV # Potassium reversal potential
g_K nS = 45 nS # Potassium peak conductance
E_Ca mV = 140 mV # Calcium reversal potential
g_Ca nS = 0.5 nS # Calcium peak conductance
g_T nS = 0.5 nS # T-type Calcium channel peak conductance
g_ahp nS = 9 nS # Afterpolarization current peak conductance
tau_syn_exc ms = 1 ms # Rise time of the excitatory synaptic alpha function
tau_syn_inh ms = 0.08 ms # Rise time of the inhibitory synaptic alpha function
E_gs mV = -85 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:
if is_refractory:
refr_t -= resolution()
V_m_old = V_m
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