From 263595df27b6dfb1d6e5ff4cfada4790c5ce7f07 Mon Sep 17 00:00:00 2001 From: "C.A.P. Linssen" Date: Mon, 13 Nov 2023 09:45:59 -0800 Subject: [PATCH] fix convolutions integration and enhance integrate_odes() to integrate specific ODEs --- models/neurons/aeif_cond_alpha_neuron.nestml | 3 +- models/neurons/aeif_cond_exp_neuron.nestml | 3 +- .../hh_cond_exp_destexhe_neuron.nestml | 41 +++++++++------ .../neurons/hh_cond_exp_traub_neuron.nestml | 37 ++++++++------ models/neurons/hh_moto_5ht_neuron.nestml | 35 ++++++++----- models/neurons/hh_psc_alpha_neuron.nestml | 9 ++-- models/neurons/hill_tononi_neuron.nestml | 22 ++++---- models/neurons/iaf_chxk_2008_neuron.nestml | 34 ++++++------- models/neurons/iaf_cond_alpha_neuron.nestml | 3 +- models/neurons/iaf_cond_beta_neuron.nestml | 3 +- models/neurons/iaf_cond_exp_neuron.nestml | 3 +- .../neurons/iaf_cond_exp_sfa_rr_neuron.nestml | 39 +++++++++------ models/neurons/iaf_psc_alpha_neuron.nestml | 3 +- models/neurons/iaf_psc_delta_neuron.nestml | 3 +- models/neurons/iaf_psc_exp_dend_neuron.nestml | 29 +++++++---- models/neurons/iaf_psc_exp_htum_neuron.nestml | 4 ++ models/neurons/iaf_psc_exp_neuron.nestml | 5 +- .../izhikevich_psc_alpha_neuron.nestml | 36 +++++++------ models/neurons/mat2_psc_exp_neuron.nestml | 50 +++++++++++-------- models/neurons/terub_gpe_neuron.nestml | 10 ++-- models/neurons/terub_stn_neuron.nestml | 11 ++-- .../neurons/traub_cond_multisyn_neuron.nestml | 11 ++-- models/neurons/traub_psc_alpha_neuron.nestml | 9 ++-- models/neurons/wb_cond_exp_neuron.nestml | 9 ++-- models/neurons/wb_cond_multisyn_neuron.nestml | 11 ++-- 25 files changed, 251 insertions(+), 172 deletions(-) diff --git a/models/neurons/aeif_cond_alpha_neuron.nestml b/models/neurons/aeif_cond_alpha_neuron.nestml index cf7b53e58..bf1e2f330 100644 --- a/models/neurons/aeif_cond_alpha_neuron.nestml +++ b/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 diff --git a/models/neurons/aeif_cond_exp_neuron.nestml b/models/neurons/aeif_cond_exp_neuron.nestml index 43440d36d..67b976212 100644 --- a/models/neurons/aeif_cond_exp_neuron.nestml +++ b/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 diff --git a/models/neurons/hh_cond_exp_destexhe_neuron.nestml b/models/neurons/hh_cond_exp_destexhe_neuron.nestml index b8ecfed54..09512aced 100644 --- a/models/neurons/hh_cond_exp_destexhe_neuron.nestml +++ b/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 diff --git a/models/neurons/hh_cond_exp_traub_neuron.nestml b/models/neurons/hh_cond_exp_traub_neuron.nestml index 4d91f906c..a366a83d7 100644 --- a/models/neurons/hh_cond_exp_traub_neuron.nestml +++ b/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 diff --git a/models/neurons/hh_moto_5ht_neuron.nestml b/models/neurons/hh_moto_5ht_neuron.nestml index a02381aa4..9770f53b9 100644 --- a/models/neurons/hh_moto_5ht_neuron.nestml +++ b/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.)) diff --git a/models/neurons/hh_psc_alpha_neuron.nestml b/models/neurons/hh_psc_alpha_neuron.nestml index 35fe62732..5940d429e 100644 --- a/models/neurons/hh_psc_alpha_neuron.nestml +++ b/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 diff --git a/models/neurons/hill_tononi_neuron.nestml b/models/neurons/hill_tononi_neuron.nestml index 0ba4eb720..7c435b805 100644 --- a/models/neurons/hill_tononi_neuron.nestml +++ b/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,14 +192,18 @@ 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 @@ -207,11 +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 PotassiumRefractoryCounts > 0: + if is_refractory: g_spike = true else: g_spike = false - r_potassium = PotassiumRefractoryCounts + potassium_refr_t = t_spike emit_spike() diff --git a/models/neurons/iaf_chxk_2008_neuron.nestml b/models/neurons/iaf_chxk_2008_neuron.nestml index 559d1dffc..ea8bc6eb4 100644 --- a/models/neurons/iaf_chxk_2008_neuron.nestml +++ b/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() diff --git a/models/neurons/iaf_cond_alpha_neuron.nestml b/models/neurons/iaf_cond_alpha_neuron.nestml index 12ef44895..10416c39d 100644 --- a/models/neurons/iaf_cond_alpha_neuron.nestml +++ b/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 diff --git a/models/neurons/iaf_cond_beta_neuron.nestml b/models/neurons/iaf_cond_beta_neuron.nestml index a4bee9a02..a5906c9cc 100644 --- a/models/neurons/iaf_cond_beta_neuron.nestml +++ b/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 diff --git a/models/neurons/iaf_cond_exp_neuron.nestml b/models/neurons/iaf_cond_exp_neuron.nestml index cbd53d846..e1dbfc895 100644 --- a/models/neurons/iaf_cond_exp_neuron.nestml +++ b/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 diff --git a/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml b/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml index 0d700f34e..7463f2c15 100644 --- a/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml +++ b/models/neurons/iaf_cond_exp_sfa_rr_neuron.nestml @@ -36,9 +36,9 @@ aeif_cond_alpha, aeif_cond_exp, iaf_chxk_2008 model iaf_cond_exp_sfa_rr_neuron: state: - r integer = 0 # counts number of tick during the refractory period - 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 @@ -60,7 +60,7 @@ model iaf_cond_exp_sfa_rr_neuron: parameters: V_th mV = -57.0 mV # Threshold potential V_reset mV = -70.0 mV # Reset potential - t_ref ms = 0.5 ms # Refractory period + refr_T ms = .5 ms # Duration of refractory period g_L nS = 28.95 nS # Leak conductance C_m pF = 289.5 pF # Membrane capacitance E_exc mV = 0 mV # Excitatory reversal potential @@ -78,9 +78,6 @@ model iaf_cond_exp_sfa_rr_neuron: # constant external input current I_e pA = 0 pA - internals: - RefractoryCounts integer = steps(t_ref) # refractory time in steps - input: inh_spikes <- inhibitory spike exc_spikes <- excitatory spike @@ -90,13 +87,23 @@ model iaf_cond_exp_sfa_rr_neuron: spike update: - integrate_odes() - if r != 0: # neuron is absolute refractory - r = r - 1 - V_m = V_reset # clamp potential - elif V_m >= V_th: # neuron is not absolute refractory - r = RefractoryCounts - V_m = V_reset # clamp potential - g_sfa += q_sfa - g_rr += q_rr - emit_spike() + if is_refractory: + # neuron is absolute refractory, do not evolve ODEs + refr_t -= resolution() + else: + # neuron not refractory, so evolve all ODEs + integrate_odes() + + onCondition(not is_refractory 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(refr_t <= resolution() / 2): + # end of refractory period + refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/iaf_psc_alpha_neuron.nestml b/models/neurons/iaf_psc_alpha_neuron.nestml index d6029ae02..d6adfd45d 100644 --- a/models/neurons/iaf_psc_alpha_neuron.nestml +++ b/models/neurons/iaf_psc_alpha_neuron.nestml @@ -105,6 +105,7 @@ model iaf_psc_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 diff --git a/models/neurons/iaf_psc_delta_neuron.nestml b/models/neurons/iaf_psc_delta_neuron.nestml index 5d49ea7cb..eea64c377 100644 --- a/models/neurons/iaf_psc_delta_neuron.nestml +++ b/models/neurons/iaf_psc_delta_neuron.nestml @@ -94,6 +94,7 @@ model iaf_psc_delta_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 diff --git a/models/neurons/iaf_psc_exp_dend_neuron.nestml b/models/neurons/iaf_psc_exp_dend_neuron.nestml index 9d9a43157..141c59f7e 100644 --- a/models/neurons/iaf_psc_exp_dend_neuron.nestml +++ b/models/neurons/iaf_psc_exp_dend_neuron.nestml @@ -38,9 +38,10 @@ iaf_cond_exp model iaf_psc_exp_dend: state: - r integer = 0 # Counts number of ticks during the refractory period V_m mV = E_L # Membrane potential I_dend pA = 0 pA # Third factor, to be read out by synapse during weight update + refr_t ms = 0 ms # Refractory period timer + is_refractory boolean = false equations: kernel I_kernel_inh = exp(-t/tau_syn_inh) @@ -53,7 +54,7 @@ model iaf_psc_exp_dend: tau_m ms = 10 ms # Membrane time constant tau_syn_inh ms = 2 ms # Time constant of inhibitory synaptic current tau_syn_exc ms = 2 ms # Time constant of excitatory synaptic current - t_ref ms = 2 ms # Duration of refractory period + refr_T ms = 2 ms # Duration of refractory period E_L mV = -70 mV # Resting potential V_reset mV = -70 mV # Reset potential of the membrane V_th mV = -55 mV # Spike threshold potential @@ -61,9 +62,6 @@ model iaf_psc_exp_dend: # constant external input current I_e pA = 0 pA - internals: - RefractoryCounts integer = steps(t_ref) # refractory time in steps - input: exc_spikes <- excitatory spike inh_spikes <- inhibitory spike @@ -73,14 +71,23 @@ model iaf_psc_exp_dend: spike update: - I_dend *= .95 + I_dend *= .95 # I_dend is always evolved, even when refractory - if r == 0: # neuron not refractory, so evolve V - integrate_odes() + if is_refractory: + # neuron is absolute refractory, do not evolve V_m + refr_t -= resolution() else: - r = r - 1 # neuron is absolute refractory + # neuron not refractory, so evolve all ODEs (including V_m) + integrate_odes() - onCondition(V_m >= V_th): # threshold crossing - r = RefractoryCounts + onCondition(not is_refractory 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 emit_spike() + + onCondition(refr_t <= resolution() / 2): + # end of refractory period + refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/iaf_psc_exp_htum_neuron.nestml b/models/neurons/iaf_psc_exp_htum_neuron.nestml index cfe62759f..c79da9558 100644 --- a/models/neurons/iaf_psc_exp_htum_neuron.nestml +++ b/models/neurons/iaf_psc_exp_htum_neuron.nestml @@ -19,6 +19,10 @@ membrane potential reaches threshold. The total refractory time must be larger or equal to the absolute refractory time. If equal, the refractoriness of the model if equivalent to the other models of NEST. +.. note:: + This neuron model can only be used in combination with a fixed + simulation resolution (timestep size). + .. note:: If tau_m is very close to tau_syn_exc or tau_syn_inh, numerical problems may arise due to singularities in the propagator matrics. If this is diff --git a/models/neurons/iaf_psc_exp_neuron.nestml b/models/neurons/iaf_psc_exp_neuron.nestml index 719261775..d8b1db5fa 100644 --- a/models/neurons/iaf_psc_exp_neuron.nestml +++ b/models/neurons/iaf_psc_exp_neuron.nestml @@ -102,13 +102,14 @@ model iaf_psc_exp_neuron: onReceive(inh_spikes): I_syn_inh += inh_spikes * pA * s - onCondition(V_m >= V_th): + onCondition(not is_refractory and 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(refr_t <= 0 ms): + onCondition(refr_t <= resolution() / 2): # end of refractory period + refr_t = 0 ms is_refractory = false diff --git a/models/neurons/izhikevich_psc_alpha_neuron.nestml b/models/neurons/izhikevich_psc_alpha_neuron.nestml index e15b683f8..2e171fd74 100644 --- a/models/neurons/izhikevich_psc_alpha_neuron.nestml +++ b/models/neurons/izhikevich_psc_alpha_neuron.nestml @@ -38,9 +38,10 @@ References """ model izhikevich_psc_alpha_neuron: state: - r integer = 0 # Number of steps in the current refractory phase V_m mV = -65 mV # Membrane potential U_m pA = 0 pA # Membrane potential recovery variable + refr_t ms = 0 ms # Refractory period timer + is_refractory boolean = false equations: # synapses: alpha functions @@ -65,14 +66,11 @@ model izhikevich_psc_alpha_neuron: V_peak mV = 0 mV # Spike detection threshold (reset condition) tau_syn_exc ms = 0.2 ms # Synaptic time constant of excitatory synapse tau_syn_inh ms = 2 ms # Synaptic time constant of inhibitory synapse - t_ref ms = 2 ms # Refractory period + refr_T ms = 2 ms # Duration of refractory period # constant external input current I_e pA = 0 pA - internals: - RefractoryCounts integer = steps(t_ref) # refractory time in steps - input: inh_spikes <- inhibitory spike exc_spikes <- excitatory spike @@ -82,13 +80,21 @@ model izhikevich_psc_alpha_neuron: spike update: - integrate_odes() - - # refractoriness and threshold crossing - if r > 0: # is refractory? - r -= 1 - elif V_m >= V_peak: - V_m = c - U_m += d - emit_spike() - r = RefractoryCounts + if is_refractory: + # neuron is absolute refractory, do not evolve V_m and U_m + refr_t -= resolution() + else: + integrate_odes() + + onCondition(not is_refractory and V_m >= V_peak): + # threshold crossing + refr_t = refr_T # start of the refractory period + is_refractory = true + V_m = c + U_m += d + emit_spike() + + onCondition(refr_t <= resolution() / 2): + # end of refractory period + refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/mat2_psc_exp_neuron.nestml b/models/neurons/mat2_psc_exp_neuron.nestml index 790409111..19700eaf5 100644 --- a/models/neurons/mat2_psc_exp_neuron.nestml +++ b/models/neurons/mat2_psc_exp_neuron.nestml @@ -47,8 +47,9 @@ model mat2_psc_exp_neuron: V_th_alpha_1 mV = 0 mV # Two-timescale adaptive threshold V_th_alpha_2 mV = 0 mV # Two-timescale adaptive threshold - r integer = 0 # counts number of tick during the refractory period V_m mV = E_L # Absolute membrane potential. + refr_t ms = 0 ms # Refractory period timer + is_refractory boolean = false equations: kernel I_kernel_inh = exp(-t/tau_syn_inh) @@ -60,10 +61,10 @@ model mat2_psc_exp_neuron: parameters: tau_m ms = 5 ms # Membrane time constant C_m pF = 100 pF # Capacitance of the membrane - t_ref ms = 2 ms # Duration of absolute refractory period (no spiking) + refr_T ms = 2 ms # Duration of refractory period E_L mV = -70 mV # Resting potential - tau_syn_exc ms = 1 ms # Time constant of postsynaptic excitatory currents - tau_syn_inh ms = 3 ms # Time constant of postsynaptic inhibitory currents + tau_syn_exc ms = 1 ms # Time constant of postsynaptic excitatory currents + tau_syn_inh ms = 3 ms # Time constant of postsynaptic inhibitory currents tau_1 ms = 10 ms # Short time constant of adaptive threshold tau_2 ms = 200 ms # Long time constant of adaptive threshold alpha_1 mV = 37 mV # Amplitude of short time threshold adaption [3] @@ -78,8 +79,6 @@ model mat2_psc_exp_neuron: P11th real = exp( -h / tau_1 ) P22th real = exp( -h / tau_2 ) - RefractoryCounts integer = steps(t_ref) # refractory time in steps - input: exc_spikes <- excitatory spike inh_spikes <- inhibitory spike @@ -89,21 +88,32 @@ model mat2_psc_exp_neuron: spike update: - # evolve membrane potential - integrate_odes() + if is_refractory: + # neuron is absolute refractory, do not evolve ODEs + refr_t -= resolution() + else: + # neuron not refractory - # evolve adaptive threshold - V_th_alpha_1 = V_th_alpha_1 * P11th - V_th_alpha_2 = V_th_alpha_2 * P22th + # evolve adaptive threshold + V_th_alpha_1 = V_th_alpha_1 * P11th + V_th_alpha_2 = V_th_alpha_2 * P22th - if r == 0: # not refractory - if V_m >= E_L + omega + V_th_alpha_1 + V_th_alpha_2: # threshold crossing - r = RefractoryCounts + # evolve V_m + integrate_odes() - # procedure for adaptive potential - V_th_alpha_1 += alpha_1 # short time - V_th_alpha_2 += alpha_2 # long time - emit_spike() - else: - r = r - 1 + onCondition(not is_refractory and V_m >= E_L + omega + V_th_alpha_1 + V_th_alpha_2): + # threshold crossing + refr_t = refr_T # start of the refractory period + is_refractory = true + + # procedure for adaptive potential + V_th_alpha_1 += alpha_1 # short time + V_th_alpha_2 += alpha_2 # long time + + emit_spike() + + onCondition(refr_t <= resolution() / 2): + # end of refractory period + refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/terub_gpe_neuron.nestml b/models/neurons/terub_gpe_neuron.nestml index fe6286a2c..a39172028 100644 --- a/models/neurons/terub_gpe_neuron.nestml +++ b/models/neurons/terub_gpe_neuron.nestml @@ -32,7 +32,7 @@ References model terub_gpe_neuron: state: V_m mV = E_L # Membrane potential - U_old mV = E_L # Membrane potential at previous timestep for threshold check + 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 @@ -146,18 +146,18 @@ model terub_gpe_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 - is_refractory = false refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/terub_stn_neuron.nestml b/models/neurons/terub_stn_neuron.nestml index 0c3504cfa..8833395e6 100644 --- a/models/neurons/terub_stn_neuron.nestml +++ b/models/neurons/terub_stn_neuron.nestml @@ -28,7 +28,7 @@ References model terub_stn_neuron: state: V_m mV = E_L # Membrane potential - U_old mV = E_L # Membrane potential at previous timestep for threshold check + 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 @@ -152,18 +152,19 @@ model terub_stn_neuron: spike update: - U_old = V_m if is_refractory: refr_t -= resolution() + + V_m_old = V_m 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 - is_refractory = false refr_t = 0 ms + is_refractory = false diff --git a/models/neurons/traub_cond_multisyn_neuron.nestml b/models/neurons/traub_cond_multisyn_neuron.nestml index cac8b79ea..8e4a0898d 100644 --- a/models/neurons/traub_cond_multisyn_neuron.nestml +++ b/models/neurons/traub_cond_multisyn_neuron.nestml @@ -28,7 +28,7 @@ hh_cond_exp_traub model traub_cond_multisyn_neuron: state: V_m mV = -70. mV # Membrane potential - U_old mV = E_L # Membrane potential at previous timestep for threshold check + 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 @@ -150,21 +150,22 @@ model traub_cond_multisyn_neuron: spike update: - U_old = V_m if is_refractory: refr_t -= resolution() + + V_m_old = V_m integrate_odes() - onCondition(not is_refractory and V_m > V_Tr and U_old > V_m): + 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(refr_t <= 0 ms): + onCondition(refr_t <= resolution() / 2): # end of refractory period - is_refractory = false refr_t = 0 ms + is_refractory = false function compute_synapse_constant(Tau_1 ms, Tau_2 ms, g_peak real) real: # Factor used to account for the missing 1/((1/Tau_2)-(1/Tau_1)) term diff --git a/models/neurons/traub_psc_alpha_neuron.nestml b/models/neurons/traub_psc_alpha_neuron.nestml index 7fdc85390..1643cb099 100644 --- a/models/neurons/traub_psc_alpha_neuron.nestml +++ b/models/neurons/traub_psc_alpha_neuron.nestml @@ -23,7 +23,7 @@ hh_cond_exp_traub model traub_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 @@ -94,18 +94,19 @@ model traub_psc_alpha_neuron: spike update: - U_old = V_m if is_refractory: refr_t -= resolution() + V_m_old = V_m 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 diff --git a/models/neurons/wb_cond_exp_neuron.nestml b/models/neurons/wb_cond_exp_neuron.nestml index 2264ec4f8..682f28959 100644 --- a/models/neurons/wb_cond_exp_neuron.nestml +++ b/models/neurons/wb_cond_exp_neuron.nestml @@ -29,7 +29,7 @@ hh_cond_exp_traub, wb_cond_multisyn model wb_cond_exp_neuron: state: V_m mV = E_L # Membrane potential - U_old mV = E_L # Membrane potential at previous timestep for threshold check + 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 @@ -86,19 +86,20 @@ model wb_cond_exp_neuron: spike update: - U_old = V_m if is_refractory: refr_t -= resolution() + V_m_old = V_m 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 function _subexpr(V_m mV) real: diff --git a/models/neurons/wb_cond_multisyn_neuron.nestml b/models/neurons/wb_cond_multisyn_neuron.nestml index c69f7524d..22ecef08d 100644 --- a/models/neurons/wb_cond_multisyn_neuron.nestml +++ b/models/neurons/wb_cond_multisyn_neuron.nestml @@ -29,7 +29,7 @@ wb_cond_multisyn model wb_cond_multisyn_neuron: state: V_m mV = -65. mV # Membrane potential - U_old mV = E_L # Membrane potential at previous timestep for threshold check + 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 @@ -138,21 +138,22 @@ model wb_cond_multisyn_neuron: spike update: - U_old = V_m if is_refractory: refr_t -= resolution() + + V_m_old = V_m integrate_odes() - onCondition(not is_refractory and V_m > V_Tr and U_old > V_m): + 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(refr_t <= 0 ms): + onCondition(refr_t <= resolution() / 2): # end of refractory period - is_refractory = false 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