fix issue with first spike passing through adaptive synapse

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;

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 )
 {
 }
 

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 )
 {
 }
 

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.

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