fix convolutions integration and enhance integrate_odes() to integrat…

…e specific ODEs

models/neurons/aeif_cond_alpha_neuron.nestml

@@ -117,7 +117,8 @@ model aeif_cond_alpha_neuron:
         w += b
         emit_spike()
 
-    onCondition(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false
 

models/neurons/aeif_cond_exp_neuron.nestml

@@ -112,6 +112,7 @@ model aeif_cond_exp_neuron:
         w += b
         emit_spike()
 
-    onCondition(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false

models/neurons/hh_cond_exp_destexhe_neuron.nestml

@@ -33,11 +33,13 @@ hh_cond_exp_traub
 model hh_cond_exp_destexhe_neuron:
 
     state:
-        r integer = 0 # counts number of tick during the refractory period
         g_noise_exc uS = g_noise_exc0
         g_noise_inh uS = g_noise_inh0
 
         V_m mV = E_L #  Membrane potential
+        V_m_old mV = E_L #  Membrane potential at the 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 )
@@ -114,11 +116,12 @@ model hh_cond_exp_destexhe_neuron:
         alpha_p_init 1/ms = 0.0001/(ms * mV) * (V_m + 30. mV) / (1. - exp(-(V_m + 30. mV) / 9. mV))
         beta_p_init 1/ms = -0.0001/(ms * mV) * (V_m + 30. mV) / (1. - exp( (V_m + 30. mV) / 9. mV ))
 
+        refr_T ms = 2 ms          # Duration of refractory period
+
         # constant external input current
         I_e pA = 0 pA
 
     internals:
-        RefractoryCounts integer = 20
         D_exc uS**2/ms = 2 * sigma_noise_exc**2 / tau_syn_exc
         D_inh uS**2/ms = 2 * sigma_noise_inh**2 / tau_syn_inh
         A_exc uS = ((D_exc * tau_syn_exc / 2) * (1 - exp(-2 * resolution() / tau_syn_exc )))**.5
@@ -133,16 +136,24 @@ model hh_cond_exp_destexhe_neuron:
         spike
 
     update:
-        U_old mV = V_m
-
-        integrate_odes()
-
-        g_noise_exc = g_noise_exc0 + (g_noise_exc - g_noise_exc0) * exp(-resolution() / tau_syn_exc) + A_exc * random_normal(0, 1)
-        g_noise_inh = g_noise_inh0 + (g_noise_inh - g_noise_inh0) * exp(-resolution() / tau_syn_inh) + A_inh * random_normal(0, 1)
-
-        # sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
-        if r > 0:
-            r -= 1
-        elif V_m > V_T + 30mV and U_old > V_m:
-            r = RefractoryCounts
-            emit_spike()
+        if is_refractory:
+            # neuron is absolute refractory, do not evolve ODEs
+            refr_t -= resolution()
+        else:
+            # neuron not refractory, so evolve all ODEs (including V_m)
+
+            g_noise_exc = g_noise_exc0 + (g_noise_exc - g_noise_exc0) * exp(-resolution() / tau_syn_exc) + A_exc * random_normal(0, 1)
+            g_noise_inh = g_noise_inh0 + (g_noise_inh - g_noise_inh0) * exp(-resolution() / tau_syn_inh) + A_inh * random_normal(0, 1)
+
+            V_m_old mV = V_m
+            integrate_odes()
+
+    onCondition(not is_refractory and V_m > V_T + 30mV and V_m_old > V_m):
+        refr_t = refr_T    # start of the refractory period
+        is_refractory = true
+        emit_spike()
+
+    onCondition(refr_t <= resolution() / 2):
+        # end of refractory period
+        refr_t = 0 ms
+        is_refractory = false

models/neurons/hh_cond_exp_traub_neuron.nestml

@@ -50,9 +50,10 @@ hh_psc_alpha
 model hh_cond_exp_traub_neuron:
 
     state:
-        r integer = 0    # counts number of tick during the refractory period
-
         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 )
@@ -97,7 +98,7 @@ model hh_cond_exp_traub_neuron:
                                   # parameters, V_T = -63 mV results in a threshold around -50 mV.
         tau_syn_exc ms = 5 ms     # Synaptic time constant of excitatory synapse
         tau_syn_inh ms = 10 ms    # Synaptic time constant of inhibitory synapse
-        t_ref ms = 2 ms           # Refractory period
+        refr_T ms = 2 ms          # Duration of refractory period
         E_exc mV = 0 mV           # Excitatory synaptic reversal potential
         E_inh mV = -80 mV         # Inhibitory synaptic reversal potential
 
@@ -111,9 +112,6 @@ model hh_cond_exp_traub_neuron:
         # constant external input current
         I_e pA = 0 pA
 
-    internals:
-        RefractoryCounts integer = steps(t_ref)
-
     input:
         inh_spikes <- inhibitory spike
         exc_spikes <- excitatory spike
@@ -123,13 +121,20 @@ model hh_cond_exp_traub_neuron:
         spike
 
     update:
-        U_old mV = V_m
-
-        integrate_odes()
-
-        # sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
-        if r > 0:
-            r -= 1
-        elif V_m > V_T + 30 mV and U_old > V_m:
-            r = RefractoryCounts
-            emit_spike()
+        if is_refractory:
+            # neuron is absolute refractory, do not evolve ODEs
+            refr_t -= resolution()
+        else:
+            # neuron not refractory, so evolve all ODEs (including V_m)
+            V_m_old mV = V_m
+            integrate_odes()
+
+    onCondition(not is_refractory 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(refr_t <= resolution() / 2):
+        # end of refractory period
+        refr_t = 0 ms
+        is_refractory = false

models/neurons/hh_moto_5ht_neuron.nestml

@@ -32,8 +32,10 @@ hh_psc_alpha
 """
 model hh_moto_5ht_neuron:
     state:
-        r integer = 0    # number of steps in the current refractory phase
         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)
@@ -70,7 +72,7 @@ model hh_moto_5ht_neuron:
         Ca_in'= (0.01 / s) * (-alpha * (I_Ca_N + I_Ca_L) - 4. * Ca_in)
 
     parameters:
-        t_ref ms = 2.0 ms                # Refractory period
+        refr_T ms = 2 ms                 # Duration of refractory period
 
         g_Na nS = 5000.0 nS              # Sodium peak conductance
         g_L nS = 200.0 nS                # Leak conductance
@@ -105,9 +107,6 @@ model hh_moto_5ht_neuron:
 
         alpha mmol/pA = 1E-5 mmol/pA
 
-    internals:
-        RefractoryCounts integer = steps(t_ref) # refractory time in steps
-
     input:
         inh_spikes <- inhibitory spike
         exc_spikes <- excitatory spike
@@ -117,14 +116,26 @@ model hh_moto_5ht_neuron:
         spike
 
     update:
-        U_old mV = V_m
+        if is_refractory:
+            # neuron is absolute refractory, do not evolve ODEs
+            refr_t -= resolution()
+        else:
+            # neuron not refractory, evolve all ODEs
+            integrate_odes()
+
+        V_m_old = V_m
         integrate_odes()
-        # sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
-        if r > 0: # is refractory?
-            r -= 1
-        elif V_m > 0 mV and U_old > V_m: # threshold && maximum
-            r = RefractoryCounts
-            emit_spike()
+
+    onCondition(not is_refractory 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(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.))

models/neurons/hh_psc_alpha_neuron.nestml

@@ -46,7 +46,7 @@ hh_cond_exp_traub
 model hh_psc_alpha_neuron:
     state:
         V_m mV = V_m_init    # Membrane potential
-        U_old mV = V_m_init    # Membrane potential at previous timestep for threshold check
+        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
 
@@ -116,18 +116,19 @@ model hh_psc_alpha_neuron:
         spike
 
     update:
-        U_old = V_m
+        V_m_old = V_m
         if is_refractory:
             refr_t -= resolution()
         integrate_odes()
 
-    onCondition(not is_refractory and V_m > 0 mV and U_old > V_m):
+    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(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false
 

models/neurons/hill_tononi_neuron.nestml

@@ -36,7 +36,8 @@ References
 """
 model hill_tononi_neuron:
     state:
-        r_potassium integer = 0
+        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
@@ -179,7 +180,6 @@ model hill_tononi_neuron:
 
         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 )
-        PotassiumRefractoryCounts integer = steps(t_spike)
 
     input:
         AMPA <- spike
@@ -192,26 +192,30 @@ model hill_tononi_neuron:
         spike
 
     update:
-        integrate_odes()
+        if is_refractory:
+            # neuron is absolute refractory, do not evolve ODEs
+            potassium_refr_t -= resolution()
+        else:
+            # neuron not refractory, so evolve all ODEs
+            integrate_odes()
 
+    onCondition(potassium_refr_t <= resolution() / 2):
         # Deactivate potassium current after spike time have expired
-        if (r_potassium > 0) and (r_potassium-1 == 0):
-            g_spike = false # Deactivate potassium current.
-        r_potassium -= 1
+        g_spike = false # Deactivate potassium current.
 
-    onCondition((not g_spike) and V_m >= Theta):
+    onCondition(not g_spike and V_m >= Theta):
         # Set V and Theta to the sodium reversal potential.
         V_m = E_Na
         Theta = E_Na
 
         # Activate fast potassium current. Drives the
         # membrane potential towards the potassium reversal
         # potential (activate only if duration is non-zero).
-        if PotassiumRefractoryCounts > 0:
+        if is_refractory:
             g_spike = true
         else:
             g_spike = false
-        r_potassium = PotassiumRefractoryCounts
+        potassium_refr_t = t_spike
 
         emit_spike()
 

models/neurons/iaf_chxk_2008_neuron.nestml

@@ -37,6 +37,7 @@ model iaf_chxk_2008_neuron:
 
     state:
         V_m mV = E_L   # membrane potential
+        V_m_prev mV = E_L   # membrane potential
         g_ahp nS = 0 nS      # AHP conductance
         g_ahp' nS/ms = 0 nS/ms   # AHP conductance
 
@@ -87,26 +88,25 @@ model iaf_chxk_2008_neuron:
         spike
 
     update:
-        vm_prev mV = V_m
+        V_m_prev = V_m
         integrate_odes()
-        if vm_prev < V_th and V_m >= V_th:
-            # Neuron is not absolute refractory
 
-            # Find precise spike time using linear interpolation
-            sigma real = ( V_m - V_th ) * resolution() / ( V_m - vm_prev ) / ms
+    onCondition(V_m_prev < V_th and V_m >= V_th):
+        # Find precise spike time using linear interpolation
+        sigma real = ( V_m - V_th ) * resolution() / ( V_m - V_m_prev ) / ms
 
-            alpha real = exp( -sigma / tau_ahp )
+        alpha real = exp( -sigma / tau_ahp )
 
-            delta_dg_ahp real = PSConInit_AHP * alpha
-            delta_g_ahp real = delta_dg_ahp * sigma  #PSConInit_AHP * sigma * alpha
+        delta_dg_ahp real = PSConInit_AHP * alpha
+        delta_g_ahp real = delta_dg_ahp * sigma
 
-            if ahp_bug == true:
-                # Bug in original code ignores AHP conductance from previous spikes
-                g_ahp = delta_g_ahp * nS
-                g_ahp' = delta_dg_ahp * nS/ms
-            else:
-                # Correct implementation adds initial values for new AHP to AHP history
-                g_ahp  += delta_g_ahp * nS
-                g_ahp' += delta_dg_ahp * nS/ms
+        if ahp_bug == true:
+            # Bug in original code ignores AHP conductance from previous spikes
+            g_ahp = delta_g_ahp * nS
+            g_ahp' = delta_dg_ahp * nS/ms
+        else:
+            # Correct implementation adds initial values for new AHP to AHP history
+            g_ahp  += delta_g_ahp * nS
+            g_ahp' += delta_dg_ahp * nS/ms
 
-            emit_spike()
+        emit_spike()

models/neurons/iaf_cond_alpha_neuron.nestml

@@ -90,7 +90,8 @@ model iaf_cond_alpha_neuron:
         V_m = V_reset
         emit_spike()
 
-    onCondition(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false
 

models/neurons/iaf_cond_beta_neuron.nestml

@@ -122,7 +122,8 @@ model iaf_cond_beta_neuron:
         V_m = V_reset
         emit_spike()
 
-    onCondition(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false
 

models/neurons/iaf_cond_exp_neuron.nestml

@@ -81,7 +81,8 @@ model iaf_cond_exp_neuron:
         V_m = V_reset
         emit_spike()
 
-    onCondition(refr_t <= 0 ms):
+    onCondition(refr_t <= resolution() / 2):
         # end of refractory period
+        refr_t = 0 ms
         is_refractory = false