Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix issue with first spike passing through adaptive synapse #3384

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Fix issue with first spike passing through adaptive synapse #3384

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b907e0d
fix issue with first spike passing through adaptive synapse
Dec 16, 2024
df296bb
Merge remote-tracking branch 'origin/master' into first_syn_spike
Feb 10, 2025
081b237
add documentation about overriding initial value in Tsodyks2 synapse
Feb 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 21 additions & 16 deletions models/quantal_stp_synapse.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion models/quantal_stp_synapse_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
{
}

Expand Down
26 changes: 18 additions & 8 deletions models/tsodyks2_synapse.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,10 @@ The parameter ``A_se`` from the publications is represented by the
synaptic weight. The variable x in the synapse properties is the
factor that scales the synaptic weight.

Please note that the initial value of ``u`` should be equal to the value of
``U``. Thus, when setting a new value for ``U`` before the start of the
simulation, make sure to set ``u`` to the same value.

.. warning::

This synaptic plasticity rule does not take
Expand Down Expand Up @@ -227,9 +231,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_ );
Comment on lines +240 to +241
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
double x_decay = std::exp( -h / tau_rec_ );
double u_decay = ( tau_fac_ < 1.0e-10 ) ? 0.0 : std::exp( -h / tau_fac_ );
const double x_decay = std::exp( -h / tau_rec_ );
const double u_decay = ( tau_fac_ == 0 ) ? 0.0 : std::exp( -h / tau_fac_ ); // tau_fac == 0 disables facilitation

I find this check here a bit confusing. Since tau_fac_ == 0 is allowed to turn off facilitation (actually the default), I think it is most sensible to actually test for zero.


// 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 );
Expand All @@ -239,10 +253,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;
Expand All @@ -257,7 +267,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 )
{
}

Expand Down
4 changes: 2 additions & 2 deletions pynest/examples/evaluate_tsodyks2_synapse.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find it strange that both dep and fac parameters are given above, but only one of the sets is used, without any further explanation. It is also not clear to me why the default is now changed from fac to dep. Wouldn't it be most illustrative to show the effect of both facilitation and depression by creating a network with in total four plastic synapses?


###############################################################################
# Create three neurons.
Expand Down
23 changes: 12 additions & 11 deletions testsuite/pytests/test_tsodyks2_synapse.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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.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"])
Comment on lines +137 to +140
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on the comment in the C++ code above, I would suggest here

Suggested change
if self.synapse_parameters["tau_fac"] < 1e-10:
u_decay = 0.0
else:
u_decay = np.exp(-h / self.synapse_parameters["tau_fac"])
if self.synapse_parameters["tau_fac"] == 0:
u_decay = 0.0
else:
u_decay = math.exp( - h / self.synapse_parameters["tau_fac"] )

Do we actually test for both cases, i.e. tau_fac == 0 and tau_fac != 0?
Also, since the exponent looks like a scalar, math.exp() is faster than the NumPy version.


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
Expand Down
Loading