From edc9f650413562406b0c5928832980250bc454aa Mon Sep 17 00:00:00 2001 From: "C.A.P. Linssen" Date: Mon, 16 Dec 2024 11:38:29 +0100 Subject: [PATCH] fix issue with first spike passing through adaptive synapse --- models/quantal_stp_synapse.h | 37 +++++++++++--------- models/quantal_stp_synapse_impl.h | 2 +- models/tsodyks2_synapse.h | 22 +++++++----- pynest/examples/evaluate_tsodyks2_synapse.py | 4 +-- testsuite/pytests/test_tsodyks2_synapse.py | 23 ++++++------ 5 files changed, 50 insertions(+), 38 deletions(-) diff --git a/models/quantal_stp_synapse.h b/models/quantal_stp_synapse.h index 3c866a3674..41b2fbf51d 100644 --- a/models/quantal_stp_synapse.h +++ b/models/quantal_stp_synapse.h @@ -209,11 +209,28 @@ inline bool quantal_stp_synapse< targetidentifierT >::send( Event& e, size_t t, const CommonSynapseProperties& ) { const double t_spike = e.get_stamp().get_ms(); - const double h = t_spike - t_lastspike_; - // Compute the decay factors, based on the time since the last spike. - const double p_decay = std::exp( -h / tau_rec_ ); - const double u_decay = ( tau_fac_ < 1.0e-10 ) ? 0.0 : std::exp( -h / tau_fac_ ); + if ( t_lastspike_ >= 0.0 ) + { + // only update a and u if this is not the first spike to pass through the synapse + + // Compute the decay factors, based on the time since the last spike. + const double h = t_spike - t_lastspike_; + const double p_decay = std::exp( -h / tau_rec_ ); + const double u_decay = ( tau_fac_ < 1.0e-10 ) ? 0.0 : std::exp( -h / tau_fac_ ); + + // Compute release probability + u_ = U_ + u_ * ( 1. - U_ ) * u_decay; // Eq. 4 from [2]_ + + // Compute number of sites that recovered during the interval. + for ( int depleted = n_ - a_; depleted > 0; --depleted ) + { + if ( get_vp_specific_rng( t )->drand() < ( 1.0 - p_decay ) ) + { + ++a_; + } + } + } // Compute number of released sites int n_release = 0; @@ -237,18 +254,6 @@ quantal_stp_synapse< targetidentifierT >::send( Event& e, size_t t, const Common a_ -= n_release; } - // Compute release probability - u_ = U_ + u_ * ( 1. - U_ ) * u_decay; // Eq. 4 from [2]_ - - // Compute number of sites that recovered during the interval. - for ( int depleted = n_ - a_; depleted > 0; --depleted ) - { - if ( get_vp_specific_rng( t )->drand() < ( 1.0 - p_decay ) ) - { - ++a_; - } - } - t_lastspike_ = t_spike; return send_spike; diff --git a/models/quantal_stp_synapse_impl.h b/models/quantal_stp_synapse_impl.h index 6e494313c8..1df262df46 100644 --- a/models/quantal_stp_synapse_impl.h +++ b/models/quantal_stp_synapse_impl.h @@ -46,7 +46,7 @@ quantal_stp_synapse< targetidentifierT >::quantal_stp_synapse() , tau_fac_( 0.0 ) , n_( 1 ) , a_( n_ ) - , t_lastspike_( 0.0 ) + , t_lastspike_( -1.0 ) { } diff --git a/models/tsodyks2_synapse.h b/models/tsodyks2_synapse.h index a57f300244..e831497148 100644 --- a/models/tsodyks2_synapse.h +++ b/models/tsodyks2_synapse.h @@ -227,9 +227,19 @@ tsodyks2_synapse< targetidentifierT >::send( Event& e, size_t t, const CommonSyn { Node* target = get_target( t ); const double t_spike = e.get_stamp().get_ms(); - const double h = t_spike - t_lastspike_; - double x_decay = std::exp( -h / tau_rec_ ); - double u_decay = ( tau_fac_ < 1.0e-10 ) ? 0.0 : std::exp( -h / tau_fac_ ); + + if ( t_lastspike_ >= 0.0 ) + { + // only update x and u if this is not the first spike to pass through the synapse + + const double h = t_spike - t_lastspike_; + double x_decay = std::exp( -h / tau_rec_ ); + double u_decay = ( tau_fac_ < 1.0e-10 ) ? 0.0 : std::exp( -h / tau_fac_ ); + + // now we compute spike number n+1 + x_ = 1. + ( x_ - x_ * u_ - 1. ) * x_decay; // Eq. 5 from reference [3]_ + u_ = U_ + u_ * ( 1. - U_ ) * u_decay; // Eq. 4 from [3]_ + } // We use the current values for the spike number n. e.set_receiver( *target ); @@ -239,10 +249,6 @@ tsodyks2_synapse< targetidentifierT >::send( Event& e, size_t t, const CommonSyn e.set_rport( get_rport() ); e(); - // now we compute spike number n+1 - x_ = 1. + ( x_ - x_ * u_ - 1. ) * x_decay; // Eq. 5 from reference [3]_ - u_ = U_ + u_ * ( 1. - U_ ) * u_decay; // Eq. 4 from [3]_ - t_lastspike_ = t_spike; return true; @@ -257,7 +263,7 @@ tsodyks2_synapse< targetidentifierT >::tsodyks2_synapse() , x_( 1.0 ) , tau_rec_( 800.0 ) , tau_fac_( 0.0 ) - , t_lastspike_( 0.0 ) + , t_lastspike_( -1.0 ) { } diff --git a/pynest/examples/evaluate_tsodyks2_synapse.py b/pynest/examples/evaluate_tsodyks2_synapse.py index d458fc1772..7bfdc06c1c 100644 --- a/pynest/examples/evaluate_tsodyks2_synapse.py +++ b/pynest/examples/evaluate_tsodyks2_synapse.py @@ -89,8 +89,8 @@ ############################################################################### # Now we assign the parameter set to the synapse models. -tsodyks_params = dict(fac_params, synapse_model="tsodyks_synapse") # for tsodyks_synapse -tsodyks2_params = dict(fac_params, synapse_model="tsodyks2_synapse") # for tsodyks2_synapse +tsodyks_params = dict(dep_params, synapse_model="tsodyks_synapse") # for tsodyks_synapse +tsodyks2_params = dict(dep_params, synapse_model="tsodyks2_synapse") # for tsodyks2_synapse ############################################################################### # Create three neurons. diff --git a/testsuite/pytests/test_tsodyks2_synapse.py b/testsuite/pytests/test_tsodyks2_synapse.py index fdc8803586..b024ee8e91 100644 --- a/testsuite/pytests/test_tsodyks2_synapse.py +++ b/testsuite/pytests/test_tsodyks2_synapse.py @@ -121,7 +121,7 @@ def reproduce_weight_drift(self, _pre_spikes, absolute_weight=1.0): n_steps = 1 + int(np.ceil(self.simulation_duration / self.resolution)) w_log = [] - t_lastspike = 0.0 + t_lastspike = -1.0 R_ = 1.0 # fraction of synaptic resources available for transmission in the range [0..1] u_ = self.synapse_parameters["U"] for time_in_simulation_steps in range(n_steps): @@ -130,20 +130,21 @@ def reproduce_weight_drift(self, _pre_spikes, absolute_weight=1.0): # Adjusting the current time to make it exact. t_spike = _pre_spikes[pre_spikes_forced_to_grid.index(time_in_simulation_steps)] - # Evaluating the depression rule. - h = t_spike - t_lastspike - R_decay = np.exp(-h / self.synapse_parameters["tau_rec"]) - if self.synapse_parameters["tau_fac"] < 1e-10: - u_decay = 0.0 - else: - u_decay = np.exp(-h / self.synapse_parameters["tau_fac"]) + if t_lastspike >= 0.: + # Evaluating the depression rule. + h = t_spike - t_lastspike + R_decay = np.exp(-h / self.synapse_parameters["tau_rec"]) + if self.synapse_parameters["tau_fac"] < 1e-10: + u_decay = 0.0 + else: + u_decay = np.exp(-h / self.synapse_parameters["tau_fac"]) + + R_ = 1.0 + (R_ - R_ * u_ - 1.0) * R_decay + u_ = self.synapse_parameters["U"] + u_ * (1.0 - self.synapse_parameters["U"]) * u_decay w = R_ * u_ * absolute_weight w_log.append(w) - R_ = 1.0 + (R_ - R_ * u_ - 1.0) * R_decay - u_ = self.synapse_parameters["U"] + u_ * (1.0 - self.synapse_parameters["U"]) * u_decay - t_lastspike = t_spike return w_log