add timestep() predefined function akin to resolution()

nest · Sep 21, 2024 · ed55a2c · ed55a2c
1 parent b6e73ce
commit ed55a2c
Show file tree

Hide file tree

Showing 22 changed files with 282 additions and 255 deletions.
diff --git a/models/neurons/hh_cond_exp_traub_neuron.nestml b/models/neurons/hh_cond_exp_traub_neuron.nestml
@@ -53,7 +53,6 @@ model hh_cond_exp_traub_neuron:
         V_m mV = E_L    #  Membrane potential
         V_m_old mV = E_L    #  Membrane potential at previous timestep
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
         Act_m real = alpha_m_init / ( alpha_m_init + beta_m_init )
         Act_h real = alpha_h_init / ( alpha_h_init + beta_h_init )
@@ -72,6 +71,7 @@ model hh_cond_exp_traub_neuron:
         inline I_syn_inh pA = convolve(g_inh, inh_spikes) * nS * ( V_m - E_inh )
 
         V_m' = ( -I_Na - I_K - I_L - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
         # channel dynamics
         inline V_rel mV = V_m - V_T
@@ -123,18 +123,13 @@ model hh_cond_exp_traub_neuron:
     update:
         # Hodgkin-Huxley type model: ODEs are always integrated, regardless of refractory state
         V_m_old = V_m
-        integrate_odes()
-
-        if is_refractory:
-            # neuron is absolute refractory, decrease refractoriness timer
-            refr_t -= resolution()
-
-    onCondition(not is_refractory and V_m > V_T + 30 mV and V_m_old > V_m):
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(V_m, Inact_n, Act_h, Act_m, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(V_m, Inact_n, Act_h, Act_m)
+
+    onCondition(refr_t <= 0 ms and V_m > V_T + 30 mV and V_m_old > V_m):
         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
diff --git a/models/neurons/hh_moto_5ht_neuron.nestml b/models/neurons/hh_moto_5ht_neuron.nestml
@@ -35,7 +35,6 @@ model hh_moto_5ht_neuron:
         V_m mV = V_m_init    # Membrane potential
         V_m_old mV = V_m_init    # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
         Ca_in mmol = Ca_in_init    # Inside Calcium concentration
         Act_m real =    alpha_m(V_m_init) / ( alpha_m(V_m_init) + beta_m(V_m_init) )
         Act_h real = h_inf(V_m_init)
@@ -61,6 +60,7 @@ model hh_moto_5ht_neuron:
         inline I_K_Ca pA = g_K_Ca_5ht * g_K_Ca * ((Ca_in * Ca_in) / (Ca_in * Ca_in + 0.014 * 0.014)) * (V_m - E_K)
 
         V_m' =( -( I_Na + I_K + I_L + I_Ca_N + I_Ca_L + I_K_Ca ) + I_stim + I_e + I_syn_inh + I_syn_exc ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
         Inact_n' = (n_inf(V_m) - Inact_n) / n_tau(V_m)
         Act_m' = alpha_m(V_m) * (1. - Act_m) - beta_m(V_m) * Act_m
@@ -69,7 +69,7 @@ model hh_moto_5ht_neuron:
         Act_mc' = (mc_inf(V_m) - Act_mc) / mc_tau
         Act_hc' = (hc_inf(V_m) - Act_hc) / hc_tau
 
-        Ca_in'= (0.01 / s) * (-alpha * (I_Ca_N + I_Ca_L) - 4. * Ca_in)
+        Ca_in' = (0.01 / s) * (-alpha * (I_Ca_N + I_Ca_L) - 4. * Ca_in)
 
     parameters:
         refr_T ms = 2 ms                 # Duration of refractory period
@@ -118,23 +118,18 @@ model hh_moto_5ht_neuron:
     update:
         # Hodgkin-Huxley type model: ODEs are always integrated, regardless of refractory state
         V_m_old = V_m
-        integrate_odes()
-
-        if is_refractory:
-            # neuron is absolute refractory, decrease refractoriness timer
-            refr_t -= resolution()
-
-    onCondition(not is_refractory and V_m > 0 mV and V_m_old > V_m):
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(V_m, Inact_n, Act_m, Act_h, Act_p, Act_mc, Act_hc, Ca_in, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(V_m, Inact_n, Act_m, Act_h, Act_p, Act_mc, Act_hc, Ca_in)
+
+    onCondition(refr_t <= 0 ms and V_m > 0 mV and V_m_old > V_m):
         # threshold && 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 h_inf(V_m mV) real:
         return 1. / (1. + exp((V_m + 65.) / 7.))
 

diff --git a/models/neurons/hh_psc_alpha_neuron.nestml b/models/neurons/hh_psc_alpha_neuron.nestml
@@ -48,7 +48,6 @@ model hh_psc_alpha_neuron:
         V_m mV = V_m_init    # Membrane potential
         V_m_old mV = V_m_init    # Membrane potential at previous timestep for threshold check
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
         Act_m real =  alpha_m_init / ( alpha_m_init + beta_m_init )   # Activation variable m for Na
         Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init )    # Inactivation variable h for Na
@@ -81,6 +80,7 @@ model hh_psc_alpha_neuron:
         Inact_h' = ( alpha_h * ( 1 - Inact_h ) - beta_h * Inact_h ) / ms # h-variable
 
         V_m' = ( -( I_Na + I_K + I_L ) + I_e + I_stim + I_syn_exc - I_syn_inh ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         V_m_init mV = -65 mV     # Initial membrane potential
@@ -118,20 +118,14 @@ model hh_psc_alpha_neuron:
     update:
         # Hodgkin-Huxley type model: ODEs are always integrated, regardless of refractory state
         V_m_old = V_m
-        integrate_odes()
-
-        if is_refractory:
-            # neuron is absolute refractory, decrease refractoriness timer
-            refr_t -= resolution()
-
-    onCondition(not is_refractory and V_m > 0 mV and V_m_old > V_m):
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(V_m, Inact_h, Act_n, Act_m, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(V_m, Inact_h, Act_n, Act_m)
+
+    onCondition(refr_t <= 0 ms 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
-
diff --git a/models/neurons/hill_tononi_neuron.nestml b/models/neurons/hill_tononi_neuron.nestml
@@ -37,7 +37,6 @@ References
 model hill_tononi_neuron:
     state:
         potassium_refr_t ms = 0 ms
-        is_refractory boolean = false
         g_spike boolean = false
 
         V_m mV = ( g_NaL * E_Na + g_KL * E_K ) / ( g_NaL + g_KL ) # membrane potential
@@ -86,6 +85,7 @@ model hill_tononi_neuron:
         # The spike current is only activate immediately after a spike.
         inline I_spike mV = (g_spike) ? -( V_m - E_K ) / Tau_spike * ms : 0 mV
         V_m'  = ( ( I_Na + I_K + I_syn + I_NaP + I_KNa + I_T + I_h + I_e + I_stim ) / Tau_m + I_spike * pA/(ms * mV) ) * s/nF
+        potassium_refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``potassium_refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
         #############
         # Intrinsic currents
@@ -192,14 +192,14 @@ model hill_tononi_neuron:
         spike
 
     update:
-        if is_refractory:
-            # neuron is absolute refractory, do not evolve ODEs
-            potassium_refr_t -= resolution()
+        if potassium_refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(potassium_refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m, Theta, IKNa_D, IT_m, IT_h, Ih_m)
 
-    onCondition(potassium_refr_t <= resolution() / 2):
+    onCondition(potassium_refr_t <= 0):
         # Deactivate potassium current after spike time have expired
         g_spike = false # Deactivate potassium current.
 
@@ -211,10 +211,11 @@ model hill_tononi_neuron:
         # Activate fast potassium current. Drives the
         # membrane potential towards the potassium reversal
         # potential (activate only if duration is non-zero).
-        if is_refractory:
+        if potassium_refr_t >= 0 ms:
             g_spike = true
         else:
             g_spike = false
+
         potassium_refr_t = t_spike
 
         emit_spike()

diff --git a/models/neurons/iaf_cond_alpha_neuron.nestml b/models/neurons/iaf_cond_alpha_neuron.nestml
@@ -39,7 +39,6 @@ model iaf_cond_alpha_neuron:
     state:
         V_m mV = E_L    # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
     equations:
         kernel g_inh = (e/tau_syn_inh) * t * exp(-t/tau_syn_inh)
@@ -50,6 +49,7 @@ model iaf_cond_alpha_neuron:
         inline I_leak pA = g_L * ( V_m - E_L )
 
         V_m' = ( -I_leak - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         C_m pF = 250 pF          # Membrane capacitance
@@ -76,22 +76,15 @@ model iaf_cond_alpha_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
+            integrate_odes(refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m)
 
     onCondition(V_m >= V_th):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-
diff --git a/models/neurons/iaf_cond_beta_neuron.nestml b/models/neurons/iaf_cond_beta_neuron.nestml
@@ -47,7 +47,6 @@ model iaf_cond_beta_neuron:
     state:
         V_m mV = E_L    # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
         # inputs from the inhibitory conductance
         g_in real = 0
@@ -69,6 +68,7 @@ model iaf_cond_beta_neuron:
         inline I_leak pA = g_L * (V_m - E_L)  # pA = nS * mV
 
         V_m' =  (-I_leak - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         C_m pF = 250 pF              # Capacitance of the membrane
@@ -108,22 +108,15 @@ model iaf_cond_beta_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
+            integrate_odes(refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m)
 
-    onCondition(V_m >= V_th):
+    onCondition(V_m >= V_th and refr_t <= 0 ms):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-
diff --git a/models/neurons/iaf_cond_exp_neuron.nestml b/models/neurons/iaf_cond_exp_neuron.nestml
@@ -30,7 +30,6 @@ model iaf_cond_exp_neuron:
     state:
         V_m mV = E_L    # Membrane potential
         refr_t ms = 0 ms    # Refractory period timer
-        is_refractory boolean = false
 
     equations:
         kernel g_inh = exp(-t/tau_syn_inh) # inputs from the inh conductance
@@ -41,6 +40,7 @@ model iaf_cond_exp_neuron:
         inline I_leak pA = g_L * ( V_m - E_L )
 
         V_m' = ( -I_leak - I_syn_exc - I_syn_inh + I_e + I_stim ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         C_m pF = 250 pF         # Membrane capacitance
@@ -67,22 +67,15 @@ model iaf_cond_exp_neuron:
         spike
 
     update:
-        if is_refractory:
+        if refr_t > 0 ms:
             # neuron is absolute refractory, do not evolve V_m
-            refr_t -= resolution()
+            integrate_odes(refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs (including V_m)
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m)
 
     onCondition(V_m >= V_th):
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false
-
diff --git a/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml b/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml
@@ -38,7 +38,6 @@ model iaf_cond_exp_sfa_rr_neuron:
     state:
         V_m mV = E_L # membrane potential
         refr_t ms = 0 ms          # Refractory period timer
-        is_refractory boolean = false
         g_sfa nS = 0 nS     # inputs from the sfa conductance
         g_rr nS = 0 nS      # inputs from the rr conductance
 
@@ -56,6 +55,7 @@ model iaf_cond_exp_sfa_rr_neuron:
         inline I_rr pA = g_rr * ( V_m - E_rr )
 
         V_m' = ( -I_L + I_e + I_stim - I_syn_exc - I_syn_inh - I_sfa - I_rr ) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
 
     parameters:
         V_th mV = -57.0 mV       # Threshold potential
@@ -87,23 +87,17 @@ model iaf_cond_exp_sfa_rr_neuron:
         spike
 
     update:
-        if is_refractory:
-            # neuron is absolute refractory, do not evolve ODEs
-            refr_t -= resolution()
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(refr_t)
         else:
-            # neuron not refractory, so evolve all ODEs
-            integrate_odes()
+            # neuron not refractory
+            integrate_odes(V_m, g_rr, g_sfa)
 
-    onCondition(not is_refractory and V_m >= V_th): # threshold crossing
+    onCondition(refr_t <= 0 ms and V_m >= V_th): # threshold crossing
         # threshold crossing
         refr_t = refr_T    # start of the refractory period
-        is_refractory = true
         V_m = V_reset
         g_sfa += q_sfa
         g_rr += q_rr
         emit_spike()
-
-    onCondition(is_refractory and refr_t <= resolution() / 2):
-        # end of refractory period
-        refr_t = 0 ms
-        is_refractory = false