From 9dc3f247e95f3928cfa563802448bdf29f323d39 Mon Sep 17 00:00:00 2001 From: sdiaz Date: Fri, 24 Nov 2017 14:18:41 +0100 Subject: [PATCH] Changed from calcium to firing rate --- nestkernel/archiving_node.cpp | 58 +++++----- nestkernel/archiving_node.h | 40 +++---- nestkernel/growth_curve.cpp | 28 ++--- nestkernel/growth_curve.h | 103 ++++++------------ nestkernel/nest_names.cpp | 6 +- nestkernel/nest_names.h | 6 +- nestkernel/node.h | 11 +- nestkernel/sp_manager.cpp | 8 +- nestkernel/synaptic_element.cpp | 8 +- nestkernel/synaptic_element.h | 37 +++---- pynest/examples/structural_plasticity.py | 54 ++++----- .../nest/tests/test_sp/test_get_sp_status.py | 4 +- .../nest/tests/test_sp/test_growth_curves.py | 100 ++++++++--------- .../test_sp/test_update_synaptic_elements.py | 12 +- 14 files changed, 216 insertions(+), 259 deletions(-) diff --git a/nestkernel/archiving_node.cpp b/nestkernel/archiving_node.cpp index 578a02ff06..aa31671b41 100644 --- a/nestkernel/archiving_node.cpp +++ b/nestkernel/archiving_node.cpp @@ -46,10 +46,10 @@ nest::Archiving_Node::Archiving_Node() , tau_minus_triplet_( 110.0 ) , tau_minus_triplet_inv_( 1. / tau_minus_triplet_ ) , last_spike_( -1.0 ) - , Ca_t_( 0.0 ) - , Ca_minus_( 0.0 ) - , tau_Ca_( 10000.0 ) - , beta_Ca_( 0.001 ) + , fr_t_( 0.0 ) + , fr_minus_( 0.0 ) + , tau_fr_( 10000.0 ) + , beta_fr_( 0.1 ) , synaptic_elements_map_() { } @@ -64,10 +64,10 @@ nest::Archiving_Node::Archiving_Node( const Archiving_Node& n ) , tau_minus_triplet_( n.tau_minus_triplet_ ) , tau_minus_triplet_inv_( n.tau_minus_inv_ ) , last_spike_( n.last_spike_ ) - , Ca_t_( n.Ca_t_ ) - , Ca_minus_( n.Ca_minus_ ) - , tau_Ca_( n.tau_Ca_ ) - , beta_Ca_( n.beta_Ca_ ) + , fr_t_( n.fr_t_ ) + , fr_minus_( n.fr_minus_ ) + , tau_fr_( n.tau_fr_ ) + , beta_fr_( n.beta_fr_ ) , synaptic_elements_map_( n.synaptic_elements_map_ ) { } @@ -178,7 +178,7 @@ nest::Archiving_Node::set_spiketime( Time const& t_sp, double offset ) { const double t_sp_ms = t_sp.get_ms() - offset; update_synaptic_elements( t_sp_ms ); - Ca_minus_ += beta_Ca_; + fr_minus_ += beta_fr_; if ( n_incoming_ ) { @@ -218,9 +218,9 @@ nest::Archiving_Node::get_status( DictionaryDatum& d ) const def< double >( d, names::t_spike, get_spiketime_ms() ); def< double >( d, names::tau_minus, tau_minus_ ); - def< double >( d, names::Ca, Ca_minus_ ); - def< double >( d, names::tau_Ca, tau_Ca_ ); - def< double >( d, names::beta_Ca, beta_Ca_ ); + def< double >( d, names::fr, fr_minus_ ); + def< double >( d, names::tau_fr, tau_fr_ ); + def< double >( d, names::beta_fr, beta_fr_ ); def< double >( d, names::tau_minus_triplet, tau_minus_triplet_ ); #ifdef DEBUG_ARCHIVER def< int >( d, names::archiver_length, history_.size() ); @@ -246,12 +246,12 @@ nest::Archiving_Node::set_status( const DictionaryDatum& d ) // We need to preserve values in case invalid values are set double new_tau_minus = tau_minus_; double new_tau_minus_triplet = tau_minus_triplet_; - double new_tau_Ca = tau_Ca_; - double new_beta_Ca = beta_Ca_; + double new_tau_fr = tau_fr_; + double new_beta_fr = beta_fr_; updateValue< double >( d, names::tau_minus, new_tau_minus ); updateValue< double >( d, names::tau_minus_triplet, new_tau_minus_triplet ); - updateValue< double >( d, names::tau_Ca, new_tau_Ca ); - updateValue< double >( d, names::beta_Ca, new_beta_Ca ); + updateValue< double >( d, names::tau_fr, new_tau_fr ); + updateValue< double >( d, names::beta_fr, new_beta_fr ); if ( new_tau_minus <= 0.0 || new_tau_minus_triplet <= 0.0 ) { @@ -263,19 +263,19 @@ nest::Archiving_Node::set_status( const DictionaryDatum& d ) tau_minus_inv_ = 1. / tau_minus_; tau_minus_triplet_inv_ = 1. / tau_minus_triplet_; - if ( new_tau_Ca <= 0.0 ) + if ( new_tau_fr <= 0.0 ) { throw BadProperty( "All time constants must be strictly positive." ); } - tau_Ca_ = new_tau_Ca; + tau_fr_ = new_tau_fr; - if ( new_beta_Ca <= 0.0 ) + if ( new_beta_fr <= 0.0 ) { throw BadProperty( - "For Ca to function as an integrator of the electrical activity, beta_ca " - "needs to be greater than 0." ); + "beta_fr needs to be greater than 0 to be able to calculate the" + "firing rate." ); } - beta_Ca_ = new_beta_Ca; + beta_fr_ = new_beta_fr; // check, if to clear spike history and K_minus bool clear = false; @@ -333,8 +333,8 @@ nest::Archiving_Node::clear_history() Kminus_ = 0.0; triplet_Kminus_ = 0.0; history_.clear(); - Ca_minus_ = 0.0; - Ca_t_ = 0.0; + fr_minus_ = 0.0; + fr_t_ = 0.0; } @@ -417,18 +417,18 @@ nest::Archiving_Node::get_synaptic_elements() const void nest::Archiving_Node::update_synaptic_elements( double t ) { - assert( t >= Ca_t_ ); + assert( t >= fr_t_ ); for ( std::map< Name, SynapticElement >::iterator it = synaptic_elements_map_.begin(); it != synaptic_elements_map_.end(); ++it ) { - it->second.update( t, Ca_t_, Ca_minus_, tau_Ca_ ); + it->second.update( t, fr_t_, fr_minus_, tau_fr_ ); } - // Update calcium concentration - Ca_minus_ = Ca_minus_ * std::exp( ( Ca_t_ - t ) / tau_Ca_ ); - Ca_t_ = t; + // Update the firing rate + fr_minus_ = fr_minus_ * std::exp( ( fr_t_ - t ) / tau_fr_ ); + fr_t_ = t; } void diff --git a/nestkernel/archiving_node.h b/nestkernel/archiving_node.h index 7ed0f4920e..9717a1591e 100644 --- a/nestkernel/archiving_node.h +++ b/nestkernel/archiving_node.h @@ -75,7 +75,7 @@ class Archiving_Node : public Node * \fn double get_Ca_minus() * return the current value of Ca_minus */ - double get_Ca_minus() const; + double get_fr_minus() const; /** * \fn double get_synaptic_elements(Name n) @@ -171,10 +171,10 @@ class Archiving_Node : public Node void set_status( const DictionaryDatum& d ); /** - * retrieve the current value of tau_Ca which defines the exponential decay - * constant of the intracellular calcium concentration + * retrieve the current value of tau_fr which defines an exponential decay + * constant for the calculation of the firing rate */ - double get_tau_Ca() const; + double get_tau_fr() const; protected: /** @@ -223,21 +223,21 @@ class Archiving_Node : public Node * Structural plasticity */ - // Time of the last update of the Calcium concentration in ms - double Ca_t_; + // Time of the last update of the firing rate in ms + double fr_t_; - // Value of the calcium concentration [Ca2+] at Ca_t_. Intracellular calcium - // concentration has a linear factor to mean electrical activity of 10^2, - // this means, for example, that a [Ca2+] of 0.2 is equivalent to a mean - // activity of 20Hz. - double Ca_minus_; + // Value of the firing rate calculated using a low pass filter algorithm. + // Each time the neuron fires, beta_FR_ is added to FR_minus_. + // Between spikes, in the function update_synaptic_elements, + // fr_minus is decreases exponentially with a time + // constant tau_FR_. + double fr_minus_; - // Time constant for exponential decay of the intracellular calcium - // concentration - double tau_Ca_; + // Time constant for exponential decay used to calculate the firing rate + double tau_fr_; - // Increase in calcium concentration [Ca2+] for each spike of the neuron - double beta_Ca_; + // Constant used to update the firing rate at each spike of the neuron + double beta_fr_; // Map of the synaptic elements std::map< Name, SynapticElement > synaptic_elements_map_; @@ -250,15 +250,15 @@ Archiving_Node::get_spiketime_ms() const } inline double -Archiving_Node::get_tau_Ca() const +Archiving_Node::get_tau_fr() const { - return tau_Ca_; + return tau_fr_; } inline double -Archiving_Node::get_Ca_minus() const +Archiving_Node::get_fr_minus() const { - return Ca_minus_; + return fr_minus_; } } // of namespace diff --git a/nestkernel/growth_curve.cpp b/nestkernel/growth_curve.cpp index 98eafa75fa..a9761ef24b 100644 --- a/nestkernel/growth_curve.cpp +++ b/nestkernel/growth_curve.cpp @@ -65,13 +65,13 @@ nest::GrowthCurveLinear::set( const DictionaryDatum& d ) double nest::GrowthCurveLinear::update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z_minus, - double tau_Ca, + double tau_fr, double growth_rate ) const { - const double Ca = Ca_minus * std::exp( ( t_minus - t ) / tau_Ca ); - const double z_value = growth_rate * tau_Ca * ( Ca - Ca_minus ) / eps_ + const double fr = fr_minus * std::exp( ( t_minus - t ) / tau_fr ); + const double z_value = growth_rate * tau_fr * ( fr - fr_minus ) / eps_ + growth_rate * ( t - t_minus ) + z_minus; return std::max( z_value, 0.0 ); @@ -106,9 +106,9 @@ nest::GrowthCurveGaussian::set( const DictionaryDatum& d ) double nest::GrowthCurveGaussian::update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z_minus, - double tau_Ca, + double tau_fr, double growth_rate ) const { // Numerical integration from t_minus to t @@ -118,13 +118,13 @@ nest::GrowthCurveGaussian::update( double t, const double xi = ( eta_ + eps_ ) / 2.0; double z_value = z_minus; - double Ca = Ca_minus; + double fr = fr_minus; for ( double lag = t_minus; lag < ( t - h / 2.0 ); lag += h ) { - Ca = Ca - ( ( Ca / tau_Ca ) * h ); + fr = fr - ( ( fr / tau_fr ) * h ); const double dz = - h * growth_rate * ( 2.0 * exp( -pow( ( Ca - xi ) / zeta, 2 ) ) - 1.0 ); + h * growth_rate * ( 2.0 * exp( -pow( ( fr - xi ) / zeta, 2 ) ) - 1.0 ); z_value = z_value + dz; } @@ -166,9 +166,9 @@ nest::GrowthCurveSigmoid::set( const DictionaryDatum& d ) double nest::GrowthCurveSigmoid::update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z_minus, - double tau_Ca, + double tau_fr, double growth_rate ) const { // Numerical integration from t_minus to t @@ -176,13 +176,13 @@ nest::GrowthCurveSigmoid::update( double t, const double h = Time::get_resolution().get_ms(); double z_value = z_minus; - double Ca = Ca_minus; + double fr = fr_minus; for ( double lag = t_minus; lag < ( t - h / 2.0 ); lag += h ) { - Ca = Ca - ( ( Ca / tau_Ca ) * h ); + fr = fr - ( ( fr / tau_fr ) * h ); const double dz = h * growth_rate - * ( ( 2.0 / ( 1.0 + exp( ( Ca - eps_ ) / psi_ ) ) ) - 1.0 ); + * ( ( 2.0 / ( 1.0 + exp( ( fr - eps_ ) / psi_ ) ) ) - 1.0 ); z_value = z_value + dz; } diff --git a/nestkernel/growth_curve.h b/nestkernel/growth_curve.h index 38235ad4ed..bb79df7bc9 100644 --- a/nestkernel/growth_curve.h +++ b/nestkernel/growth_curve.h @@ -43,7 +43,7 @@ namespace nest /** * \class GrowthCurve * Defines the way the number of synaptic elements changes through time - * according to the calcium concentration of the neuron. + * according to the firing rate of the neuron. */ class GrowthCurve { @@ -55,9 +55,9 @@ class GrowthCurve virtual void set( const DictionaryDatum& d ) = 0; virtual double update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z, - double tau_Ca, + double tau_fr, double growth_rate ) const = 0; virtual bool is( Name n ) @@ -87,26 +87,15 @@ class GrowthCurve when structural plasticity is enabled, allows the dynamic rewiring of the network during the simulation. This type of growth curve uses an exact integration method to update the - number of synaptic elements: dz/dt = nu (1 - (1/eps) * Ca(t)), where nu is - the growth rate [elements/ms] and eps is the desired average calcium - concentration. The growth rate nu is defined in the SynapticElement class. + number of synaptic elements: dz/dt = nu (1 - (1/eps) * fr(t)), where nu is + the growth rate [elements/ms] and eps is the desired firing rate. + The growth rate nu is defined in the SynapticElement class. Parameters: - eps double - The target calcium concentration that + eps double - The target firing rate in Hz that the neuron should look to achieve by creating or deleting synaptic elements. It should always be a - positive value. It is important to note that the - calcium concentration is linearly proportional to the - firing rate. This is because dCa/dt = - Ca(t)/tau_Ca - + beta_Ca if the neuron fires and dCa/dt = - - Ca(t)/tau_Ca otherwise, where tau_Ca is the calcium - concentration decay constant and beta_Ca is the - calcium intake constant (see SynapticElement class). - This means that eps also defines the desired firing - rate that the neuron should achieve. For example, an - eps = 0.05 [Ca2+] with tau_Ca = 10000.0 and beta_Ca = - 0.001 for a synaptic element means a desired firing - rate of 5Hz. + positive value. (see SynapticElement class). References: [1] Butz, Markus, Florentin Wörgötter, and Arjen van Ooyen. @@ -125,8 +114,8 @@ class GrowthCurve /** * \class GrowthCurveLinear * Uses an exact integration method to update the number of synaptic elements: - * dz/dt = nu (1 - (1/eps) * Ca(t)), where nu is the growth rate and - * eps is the desired average calcium concentration. + * dz/dt = nu (1 - (1/eps) * fr(t)), where nu is the growth rate and + * eps is the desired firing rate in Hz. */ class GrowthCurveLinear : public GrowthCurve { @@ -136,9 +125,9 @@ class GrowthCurveLinear : public GrowthCurve void set( const DictionaryDatum& d ); double update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z, - double tau_Ca, + double tau_fr, double growth_rate ) const; private: @@ -155,40 +144,29 @@ class GrowthCurveLinear : public GrowthCurve network during the simulation. This type of growth curve uses a forward Euler integration method to update the number of synaptic elements: - dz/dt = nu (2 * e^(- ((Ca(t) - xi)/z)^2 ) - 1) + dz/dt = nu (2 * e^(- ((fr(t) - xi)/z)^2 ) - 1) where xi = (eta + eps)/2, zeta = (eps - eta)/2 * sqrt(ln(2))), - eta is the minimum calcium concentration required for any synaptic element - to be created, eps is the target mean calcium concentration in the neuron + eta is the minimum firing rate required for any synaptic element + to be created, eps is the target firing rate in the neuron and nu is the growth rate in elements/ms. The growth rate nu is defined in the SynapticElement class. Parameters: - eta double - Minimum amount of calcium concentration that the + eta double - Minimum firing rate in Hz that the neuron needs to start creating synaptic elements. eta can have a negative value, making the growth curve move its maximum to the left. For example, if - eta=-0.5 and eps=0.5 [Ca2+], the maximum growth rate - (elements/ms) will be achieved at 0.0 [Ca2+]. If - eta=0.0 [Ca2+] and eps=0.5 [Ca2+] the maximum growth - rate will be achieved at 0.25 [Ca2+] while at 0.0 - [Ca+2] no new elements will be created. + eta=-5 and eps=5 Hz, the maximum growth rate + (elements/ms) will be achieved at 0.0 Hz. If + eta=0.0 and eps=5 Hz the maximum growth + rate will be achieved at 2.5 Hz while at 0.0 + Hz no new elements will be created. - eps double - The target calcium concentration that + eps double - The target firing rate in Hz that the neuron should look to achieve by creating or deleting synaptic elements. It should always be a - positive value. It is important to note that the - calcium concentration is linearly proportional to the - firing rate. This is because dCa/dt = - Ca(t)/tau_Ca - + beta_Ca if the neuron fires and dCa/dt = - - Ca(t)/tau_Ca otherwise, where tau_Ca is the calcium - concentration decay constant and beta_Ca is the - calcium intake constant (see SynapticElement class). - This means that eps can also be seen as the desired - firing rate that the neuron should achieve. For - example, an eps = 0.05 [Ca2+] with tau_Ca = 10000.0 - and beta_Ca = 0.001 for a synaptic element means a - desired firing rate of 5Hz. + positive value. nu double - Growth rate in elements/ms. The growth rate nu is defined in the SynapticElement class. Can be negative. @@ -214,8 +192,8 @@ class GrowthCurveLinear : public GrowthCurve * dz/dt = nu (2 * e^(- ((Ca(t) - xi)/zeta)^2 ) - 1) * where xi = (eta + eps)/2, * zeta = (eps - eta)/2 * sqrt(ln(2))), - * eta is the minimum calcium concentration required for any synaptic element - * to be created, eps is the target mean calcium concentration in the + * eta is the minimum firing rate in Hz required for any synaptic element + * to be created, eps is the target firing rate in Hz in the * neuron and nu is the growth rate. */ class GrowthCurveGaussian : public GrowthCurve @@ -226,9 +204,9 @@ class GrowthCurveGaussian : public GrowthCurve void set( const DictionaryDatum& d ); double update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z, - double tau_Ca, + double tau_fr, double growth_rate ) const; private: @@ -246,27 +224,16 @@ class GrowthCurveGaussian : public GrowthCurve network during the simulation. This type of growth curve uses a forward Euler integration method to update the number of synaptic elements: - dz/dt = nu ((2 / (1 + e^((Ca(t) - eps)/psi))) - 1) - eps is the target mean calcium concentration in the + dz/dt = nu ((2 / (1 + e^((fr(t) - eps)/psi))) - 1) + eps is the target firing rate in Hz in the neuron, psi controls the width of the sigmoid and nu is the growth rate in elements/ms. The growth rate nu is defined in the SynapticElement class. Parameters: - eps double - The target calcium concentration that + eps double - The target firing rate in Hz that the neuron should look to achieve by creating or deleting synaptic elements. It should always be a - positive value. It is important to note that the - calcium concentration is linearly proportional to the - firing rate. This is because dCa/dt = - Ca(t)/tau_Ca - + beta_Ca if the neuron fires and dCa/dt = - - Ca(t)/tau_Ca otherwise, where tau_Ca is the calcium - concentration decay constant and beta_Ca is the - calcium intake constant (see SynapticElement class). - This means that eps can also be seen as the desired - firing rate that the neuron should achieve. For - example, an eps = 0.05 [Ca2+] with tau_Ca = 10000.0 - and beta_Ca = 0.001 for a synaptic element means a - desired firing rate of 5Hz. + positive value. nu double - Growth rate in elements/ms. The growth rate nu is defined in the SynapticElement class. Can be negative. @@ -288,8 +255,8 @@ class GrowthCurveGaussian : public GrowthCurve * \class GrowthCurveSigmoid * Uses a forward Euler integration method to update the number of synaptic * elements: - * dz/dt = nu ((2 / (1 + e^((Ca(t) - eps)/psi))) - 1) - * eps is the target mean calcium concentration in the + * dz/dt = nu ((2 / (1 + e^((fr(t) - eps)/psi))) - 1) + * eps is the target firing rate in Hz in the * neuron, psi controls the width of the sigmoid * and nu is the growth rate. */ @@ -301,9 +268,9 @@ class GrowthCurveSigmoid : public GrowthCurve void set( const DictionaryDatum& d ); double update( double t, double t_minus, - double Ca_minus, + double fr_minus, double z, - double tau_Ca, + double tau_fr, double growth_rate ) const; private: diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 8b9ed49b30..8228742b4c 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -62,7 +62,7 @@ const Name autapses( "autapses" ); const Name b( "b" ); const Name beta( "beta" ); -const Name beta_Ca( "beta_Ca" ); +const Name beta_fr( "beta_fr" ); const Name binary( "binary" ); const Name c( "c" ); @@ -70,7 +70,6 @@ const Name c_1( "c_1" ); const Name c_2( "c_2" ); const Name c_3( "c_3" ); const Name C_m( "C_m" ); -const Name Ca( "Ca" ); const Name calibrate( "calibrate" ); const Name calibrate_node( "calibrate_node" ); const Name capacity( "capacity" ); @@ -164,6 +163,7 @@ const Name filename( "filename" ); const Name filenames( "filenames" ); const Name flush_after_simulate( "flush_after_simulate" ); const Name flush_records( "flush_records" ); +const Name fr( "fr" ); const Name frequency( "frequency" ); const Name frozen( "frozen" ); @@ -423,7 +423,7 @@ const Name tau_1( "tau_1" ); const Name tau_2( "tau_2" ); const Name tau_ahp( "tau_ahp" ); const Name tau_c( "tau_c" ); -const Name tau_Ca( "tau_Ca" ); +const Name tau_fr( "tau_fr" ); const Name tau_D_KNa( "tau_D_KNa" ); const Name tau_decay( "tau_decay" ); const Name tau_decay_AMPA( "tau_decay_AMPA" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index dbcea8548b..d8a7943ca3 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -76,7 +76,7 @@ extern const Name autapses; //!< Connectivity-related extern const Name b; //!< Specific to Brette & Gerstner 2005 (aeif_cond-*) extern const Name beta; //!< Specific to amat2_* extern const Name - beta_Ca; //!< Increment in calcium concentration with each spike + beta_fr; //!< Constant for firing rate calculations extern const Name binary; //!< Recorder parameter extern const Name c; //!< Specific to Izhikevich 2003 @@ -84,7 +84,6 @@ extern const Name c_1; //!< Specific to stochastic neuron pp_psc_delta extern const Name c_2; //!< Specific to stochastic neuron pp_psc_delta extern const Name c_3; //!< Specific to stochastic neuron pp_psc_delta extern const Name C_m; //!< Membrane capacitance -extern const Name Ca; //!< Calcium concentration extern const Name calibrate; //!< Command to calibrate the neuron (sli_neuron) extern const Name calibrate_node; //!< Command to calibrate the neuron (sli_neuron) @@ -191,6 +190,7 @@ extern const Name filenames; //!< Recorder parameter---keep, will disappear with NESTIO extern const Name flush_after_simulate; //!< Recorder parameter extern const Name flush_records; //!< Recorder parameter +extern const Name fr; //!< Firing rate extern const Name frequency; //!< Signal modulation frequency extern const Name frozen; //!< Node parameter @@ -483,7 +483,7 @@ extern const Name tau; //!< Used by stdp_connection_facetshw_hom extern const Name tau_1; //!< Specific to Kobayashi, Tsubo, Shinomoto 2009 extern const Name tau_2; //!< Specific to Kobayashi, Tsubo, Shinomoto 2009 extern const Name tau_ahp; //!< Specific to iaf_chxk_2008 neuron -extern const Name tau_Ca; //!< Rate of loss of calcium concentration +extern const Name tau_fr; //!< Decay constant firing rate calculations extern const Name tau_c; //!< Used by stdp_dopa_connection extern const Name tau_D_KNa; //!< specific to Hill & Tononi 2005 extern const Name tau_decay; //!< Synapse decay constant (beta fct decay) diff --git a/nestkernel/node.h b/nestkernel/node.h index 942c758e63..ab17118e16 100644 --- a/nestkernel/node.h +++ b/nestkernel/node.h @@ -609,20 +609,18 @@ class Node */ /** - * Return the Ca_minus value at time Ca_t which corresponds to the time of - * the last update in Calcium concentration which is performed each time - * a Node spikes. + * Return the firing rate at time fr_t * Return 0.0 if not overridden * @ingroup SP_functions */ virtual double - get_Ca_minus() const + get_fr_minus() const { return 0.0; } /** - * Get the number of synaptic element for the current Node at Ca_t which + * Get the number of synaptic element for the current Node at fr_t which * corresponds to the time of the last spike. * Return 0.0 if not overridden * @ingroup SP_functions @@ -665,8 +663,7 @@ class Node /** * Triggers the update of all SynapticElements - * stored in the synaptic_element_map_. It also updates the calcium - * concentration. + * stored in the synaptic_element_map_. It also updates the firing rate. * @param t double time when the update is being performed * @ingroup SP_functions */ diff --git a/nestkernel/sp_manager.cpp b/nestkernel/sp_manager.cpp index d6041095d3..9228ff3e52 100644 --- a/nestkernel/sp_manager.cpp +++ b/nestkernel/sp_manager.cpp @@ -425,11 +425,11 @@ SPManager::update_structural_plasticity( SPBuilder* sp_builder ) sp_builder->get_pre_synaptic_element_name(), sp_builder->get_post_synaptic_element_name() ); // update the number of synaptic elements - get_synaptic_elements( sp_builder->get_pre_synaptic_element_name(), + /*get_synaptic_elements( sp_builder->get_pre_synaptic_element_name(), pre_vacant_id, pre_vacant_n, pre_deleted_id, - pre_deleted_n ); + pre_deleted_n );*/ } // Get post synaptic elements data from local nodes get_synaptic_elements( sp_builder->get_post_synaptic_element_name(), @@ -450,7 +450,7 @@ SPManager::update_structural_plasticity( SPBuilder* sp_builder ) sp_builder->get_synapse_model(), sp_builder->get_pre_synaptic_element_name(), sp_builder->get_post_synaptic_element_name() ); - get_synaptic_elements( sp_builder->get_pre_synaptic_element_name(), + /*get_synaptic_elements( sp_builder->get_pre_synaptic_element_name(), pre_vacant_id, pre_vacant_n, pre_deleted_id, @@ -459,7 +459,7 @@ SPManager::update_structural_plasticity( SPBuilder* sp_builder ) post_vacant_id, post_vacant_n, post_deleted_id, - post_deleted_n ); + post_deleted_n );*/ } // Communicate vacant elements diff --git a/nestkernel/synaptic_element.cpp b/nestkernel/synaptic_element.cpp index e6c049857a..15bdb5e813 100644 --- a/nestkernel/synaptic_element.cpp +++ b/nestkernel/synaptic_element.cpp @@ -149,15 +149,15 @@ nest::SynapticElement::set( const DictionaryDatum& d ) void nest::SynapticElement::update( double t, double t_minus, - double Ca_minus, - double tau_Ca ) + double fr_minus, + double tau_fr ) { if ( z_t_ != t_minus ) { throw KernelException( - "Last update of the calcium concentration does not match the last update " + "Last update of the firing rate does not match the last update " "of the synaptic element" ); } - z_ = growth_curve_->update( t, t_minus, Ca_minus, z_, tau_Ca, growth_rate_ ); + z_ = growth_curve_->update( t, t_minus, fr_minus, z_, tau_fr, growth_rate_ ); z_t_ = t; } diff --git a/nestkernel/synaptic_element.h b/nestkernel/synaptic_element.h index ee68fda931..a393a587ae 100644 --- a/nestkernel/synaptic_element.h +++ b/nestkernel/synaptic_element.h @@ -36,37 +36,31 @@ and deletion of synapses. Description: - This class represents synaptic element of a node (like Axonl boutons or + This class represents synaptic element of a node (like Axonal boutons or dendritic spines) used for structural plasticity. The synaptic elements represent connection points between two neurons. They grow according to a homeostatic growth rule. The dynamics of the - number of synaptic elements is driven by the average electrical activity of - the neuron (indirectly measured through the Calcium concentration of the - node). The probability of two neurons creating a new synapse between them, + number of synaptic elements is driven by the average firing rate of the + neuron. The probability of two neurons creating a new synapse between them, depends on the number of available synaptic elements of each neuron. Parameters: z double - Current number of synaptic elements. Stored as a double variable but the actual usable number of - synaptic elements is an integer truncated from - this + synaptic elements is an int truncated from this double value. An standard value for the growth of - a - synaptic element is around 0.0001 elements/ms. + synaptic elements is around 0.0001 elements/ms. continuous boolean - Defines if the number of synaptic elements should be treated as a continuous double number or as an integer value. Default is false. growth_rate double - The maximum amount by which the synaptic elements - will change between time steps. In elements/ms. tau_vacant double - Rate at which vacant synaptic elements will decay. Typical is 0.1 which represents a loss of 10% of the vacant synaptic elements each - time - the structural_plasticity_update_interval is - reached by the simulation time. - growth_curve GrowthCurve* - Rule which defines the dynamics of this - synaptic element. + time structural_plasticity_update_interval is + called by the simulation. + growth_curve GrowthCurve* - Defines the changes of this synaptic element. References: [1] Butz, Markus, Florentin Wörgötter, and Arjen van Ooyen. @@ -101,9 +95,8 @@ namespace nest * The synaptic elements represent connection points between two neurons that * grow according to a homeostatic growth rule. Basically, the dynamics of the * number of synaptic elements is driven by the average electrical activity of - * the neuron (indirectly measured through the Calcium concentration of the - * node). The probability of two neurons creating a new synapse between them, - * depends on the number of available synaptic elements of each neuron. + * the neuron. The probability of two neurons creating a new synapse between + * them, depends on the number of available synaptic elements of each neuron. */ class SynapticElement { @@ -154,14 +147,14 @@ class SynapticElement /* - * Updates the number of available synaptic elements according to the mean - * calcium concentration of the neuron at time t. + * Updates the number of available synaptic elements according to the + * firing rate of the neuron at time t. * @param t Current time (in ms) * @param t_minus Time of last update - * @param Ca_minus Calcium concentration at time t_minus - * @param tau_Ca change in the calcium concentration on each spike + * @param fr_minus firing rate at time t_minus + * @param tau_fr decay constant for firing rate calculations */ - void update( double t, double t_minus, double Ca_minus, double tau_Ca ); + void update( double t, double t_minus, double fr_minus, double tau_fr ); /** * \fn double get_z_value(Archiving_Node const *a, double t) const diff --git a/pynest/examples/structural_plasticity.py b/pynest/examples/structural_plasticity.py index 8cbea654e1..762ad6d48d 100644 --- a/pynest/examples/structural_plasticity.py +++ b/pynest/examples/structural_plasticity.py @@ -27,14 +27,14 @@ 20% inhibitory. The simulation starts without any connectivity. A set of homeostatic rules are defined, according to which structural plasticity will create and delete synapses dynamically during the simulation until a desired -level of electrical activity is reached. The model of structural plasticity +firing rate is reached. The model of structural plasticity used here corresponds to the formulation presented in Butz, M., & van Ooyen, A. (2013). A simple rule for dendritic spine and axonal bouton formation can account for cortical reorganization after focal retinal lesions. PLoS Comput. Biol. 9 (10), e1003259. At the end of the simulation, a plot of the evolution of the connectivity -in the network and the average calcium concentration in the neurons is created. +in the network and the firing rate in the neurons is created. ''' import nest @@ -79,8 +79,8 @@ def __init__(self): 'growth_curve': "gaussian", 'growth_rate': 0.0001, # (elements/ms) 'continuous': False, - 'eta': 0.0, # Ca2+ - 'eps': 0.05, # Ca2+ + 'eta': 0.0, # Hz + 'eps': 5.0, # Hz } # Inhibitory synaptic elements of excitatory neurons @@ -88,8 +88,8 @@ def __init__(self): 'growth_curve': "gaussian", 'growth_rate': 0.0001, # (elements/ms) 'continuous': False, - 'eta': 0.0, # Ca2+ - 'eps': self.growth_curve_e_e['eps'], # Ca2+ + 'eta': 0.0, # Hz + 'eps': self.growth_curve_e_e['eps'], # Hz } # Excitatory synaptic elements of inhibitory neurons @@ -97,8 +97,8 @@ def __init__(self): 'growth_curve': "gaussian", 'growth_rate': 0.0004, # (elements/ms) 'continuous': False, - 'eta': 0.0, # Ca2+ - 'eps': 0.2, # Ca2+ + 'eta': 0.0, # Hz + 'eps': 20.0, # Hz } # Inhibitory synaptic elements of inhibitory neurons @@ -106,8 +106,8 @@ def __init__(self): 'growth_curve': "gaussian", 'growth_rate': 0.0001, # (elements/ms) 'continuous': False, - 'eta': 0.0, # Ca2+ - 'eps': self.growth_curve_i_e['eps'] # Ca2+ + 'eta': 0.0, # Hz + 'eps': self.growth_curve_i_e['eps'] # Hz } ''' @@ -127,8 +127,8 @@ def __init__(self): self.nodes_e = None self.nodes_i = None - self.mean_ca_e = [] - self.mean_ca_i = [] + self.mean_fr_e = [] + self.mean_fr_i = [] self.total_connections_e = [] self.total_connections_i = [] @@ -236,18 +236,18 @@ def connect_external_input(self): {'weight': self.psc_ext, 'delay': 1.0}) ''' - In order to save the amount of average calcium concentration in each - population through time we create the function record_ca. Here we use the + In order to save the average firing rate in each + population through time we create the function record_fr. Here we use the GetStatus function to retrieve the value of Ca for every neuron in the network and then store the average. ''' - def record_ca(self): - ca_e = nest.GetStatus(self.nodes_e, 'Ca'), # Calcium concentration - self.mean_ca_e.append(numpy.mean(ca_e)) + def record_fr(self): + fr_e = nest.GetStatus(self.nodes_e, 'fr'), # Calcium concentration + self.mean_fr_e.append(numpy.mean(fr_e)) - ca_i = nest.GetStatus(self.nodes_i, 'Ca'), # Calcium concentration - self.mean_ca_i.append(numpy.mean(ca_i)) + fr_i = nest.GetStatus(self.nodes_i, 'fr'), # Calcium concentration + self.mean_fr_i.append(numpy.mean(fr_i)) ''' In order to save the state of the connectivity in the network through time @@ -275,15 +275,15 @@ def plot_data(self): fig, ax1 = pl.subplots() ax1.axhline(self.growth_curve_e_e['eps'], linewidth=4.0, color='#9999FF') - ax1.plot(self.mean_ca_e, 'b', - label='Ca Concentration Excitatory Neurons', linewidth=2.0) + ax1.plot(self.mean_fr_e, 'b', + label='Firing rate excitatory neurons', linewidth=2.0) ax1.axhline(self.growth_curve_i_e['eps'], linewidth=4.0, color='#FF9999') - ax1.plot(self.mean_ca_i, 'r', - label='Ca Concentration Inhibitory Neurons', linewidth=2.0) - ax1.set_ylim([0, 0.275]) + ax1.plot(self.mean_fr_i, 'r', + label='Firing rate inhibitory neurons', linewidth=2.0) + ax1.set_ylim([0, 27.5]) ax1.set_xlabel("Time in [s]") - ax1.set_ylabel("Ca concentration") + ax1.set_ylabel("Firing rate [Hz]") ax2 = ax1.twinx() ax2.plot(self.total_connections_e, 'm', label='Excitatory connections', linewidth=2.0, linestyle='--') @@ -300,7 +300,7 @@ def plot_data(self): function we first enable structural plasticity in the network and then we simulate in steps. On each step we record the calcium concentration and the connectivity. At the end of the simulation, the plot of connections and - calcium concentration through time is generated. + firing rate through time is generated. ''' def simulate(self): @@ -312,7 +312,7 @@ def simulate(self): sim_steps = numpy.arange(0, self.t_sim, self.record_interval) for i, step in enumerate(sim_steps): nest.Simulate(self.record_interval) - self.record_ca() + self.record_fr() self.record_connectivity() if i % 20 == 0: print("Progress: " + str(i / 2) + "%") diff --git a/pynest/nest/tests/test_sp/test_get_sp_status.py b/pynest/nest/tests/test_sp/test_get_sp_status.py index bb51f4bd26..c581bac1c7 100644 --- a/pynest/nest/tests/test_sp/test_get_sp_status.py +++ b/pynest/nest/tests/test_sp/test_get_sp_status.py @@ -50,8 +50,8 @@ class TestGetStructuralPlasticityStatus(unittest.TestCase): 'growth_curve': "gaussian", 'growth_rate': 0.0001, # (elements/ms) 'continuous': False, - 'eta': 0.0, # Ca2+ - 'eps': 0.05 + 'eta': 0.0, # Hz + 'eps': 5.0 } ''' diff --git a/pynest/nest/tests/test_sp/test_growth_curves.py b/pynest/nest/tests/test_sp/test_growth_curves.py index e0a628b243..d0d0862b92 100644 --- a/pynest/nest/tests/test_sp/test_growth_curves.py +++ b/pynest/nest/tests/test_sp/test_growth_curves.py @@ -34,19 +34,19 @@ class SynapticElementIntegrator(object): """ Generic class which describes how to compute the number of - Synaptic Element based on Ca value + Synaptic Element based on the firing rate Each derived class should overwrite the get_se(self, t) method """ - def __init__(self, tau_ca=10000.0, beta_ca=0.001): + def __init__(self, tau_fr=10000.0, beta_fr=0.1): """ Constructor - :param tau_ca (float): time constant of Ca decay - :param beta_ca (float): each spike increase Ca value by this value + :param tau_fr (float): time constant of Ca decay + :param beta_fr (float): each spike increase Ca value by this value """ - self.tau_ca = tau_ca - self.beta_ca = beta_ca + self.tau_fr = tau_fr + self.beta_fr = beta_fr self.t_minus = 0 self.ca_minus = 0 @@ -59,7 +59,7 @@ def reset(self): def handle_spike(self, t): """ - Add beta_ca to the value of Ca at t = spike time + Add beta_fr to the value of fr at t = spike time Also update the number of synaptic element :param t (float): spike time @@ -67,20 +67,20 @@ def handle_spike(self, t): assert t >= self.t_minus # Update the number of synaptic element self.se_minus = self.get_se(t) - # update Ca value - self.ca_minus = self.get_ca(t) + self.beta_ca + # update fr value + self.fr_minus = self.get_fr(t) + self.beta_fr self.t_minus = t - def get_ca(self, t): + def get_fr(self, t): """ :param t (float): current time - :return: Ca value + :return: firing rate (fr) value """ assert t >= self.t_minus - ca = self.ca_minus * math.exp((self.t_minus - t) / self.tau_ca) - if ca > 0: - return ca + fr = self.fr_minus * math.exp((self.t_minus - t) / self.tau_fr) + if fr > 0: + return fr else: return 0 @@ -99,7 +99,7 @@ class LinearExactSEI(SynapticElementIntegrator): """ Compute the number of synaptic element corresponding to a linear growth curve - dse/dCa = nu * (1 - Ca/eps) + dse/dfr = nu * (1 - fr/eps) Use the exact solution """ @@ -123,8 +123,8 @@ def get_se(self, t): """ assert t >= self.t_minus se = 1 / self.eps * ( - self.growth_rate * self.tau_ca * ( - self.get_ca(t) - self.ca_minus + self.growth_rate * self.tau_fr * ( + self.get_fr(t) - self.fr_minus ) + self.growth_rate * self.eps * (t - self.t_minus) ) + self.se_minus if se > 0: @@ -137,7 +137,7 @@ class LinearNumericSEI(SynapticElementIntegrator): """ Compute the number of synaptic element corresponding to a linear growth curve - dse/dCa = nu * (1 - Ca/eps) + dse/dfr = nu * (1 - fr/eps) Use numerical integration (see scipy.integrate.quad) """ @@ -167,14 +167,14 @@ def get_se(self, t): return 0 def growth_curve(self, t): - return self.growth_rate * (1.0 - (self.get_ca(t) / self.eps)) + return self.growth_rate * (1.0 - (self.get_fr(t) / self.eps)) class GaussianNumericSEI(SynapticElementIntegrator): """ Compute the number of synaptic element corresponding to a linear growth curve - dse/dCa = nu * (2 * exp( ((Ca - xi)/zeta)^2 ) - 1) + dse/dfr = nu * (2 * exp( ((fr - xi)/zeta)^2 ) - 1) with: xi = (eta + eps) / 2.0 zeta = (eta - eps) / (2.0 * sqrt(ln(2.0))) @@ -212,7 +212,7 @@ def get_se(self, t): def growth_curve(self, t): return self.growth_rate * ( 2 * math.exp( - - math.pow((self.get_ca(t) - self.xi) / self.zeta, 2) + - math.pow((self.get_fr(t) - self.xi) / self.zeta, 2) ) - 1 ) @@ -222,7 +222,7 @@ class SigmoidNumericSEI(SynapticElementIntegrator): """ Compute the number of synaptic element corresponding to a sigmoid growth curve - dse/dCa = nu * ((2.0 / exp( (Ca - eps)/psi)) - 1.0) + dse/dfr = nu * ((2.0 / exp( (fr - eps)/psi)) - 1.0) Use numerical integration (see scipy.integrate.quad) """ @@ -257,7 +257,7 @@ def get_se(self, t): def growth_curve(self, t): return self.growth_rate * ( (2.0 / (1.0 + math.exp( - (self.get_ca(t) - self.eps) / self.psi + (self.get_fr(t) - self.eps) / self.psi ))) - 1.0 ) @@ -282,8 +282,8 @@ def setUp(self): self.se_integrator = [] self.sim_steps = None - self.ca_nest = None - self.ca_python = None + self.fr_nest = None + self.fr_python = None self.se_nest = None self.se_python = None @@ -298,9 +298,9 @@ def setUp(self): def simulate(self): self.sim_steps = numpy.arange(0, self.sim_time, self.sim_step) - self.ca_nest = numpy.zeros( + self.fr_nest = numpy.zeros( (len(self.local_nodes), len(self.sim_steps))) - self.ca_python = numpy.zeros( + self.fr_python = numpy.zeros( (len(self.se_integrator), len(self.sim_steps))) self.se_nest = numpy.zeros( (len(self.local_nodes), len(self.sim_steps))) @@ -310,8 +310,8 @@ def simulate(self): start = time.clock() for t_i, t in enumerate(self.sim_steps): for n_i, n in enumerate(self.local_nodes): - self.ca_nest[n_i][t_i], synaptic_elements = nest.GetStatus( - [n], ('Ca', 'synaptic_elements'))[0] + self.fr_nest[n_i][t_i], synaptic_elements = nest.GetStatus( + [n], ('fr', 'synaptic_elements'))[0] self.se_nest[n_i][t_i] = synaptic_elements['se']['z'] nest.Simulate(self.sim_step) @@ -329,12 +329,12 @@ def simulate(self): for sei in self.se_integrator] spike_i += 1 for sei_i, sei in enumerate(self.se_integrator): - self.ca_python[sei_i, t_i] = sei.get_ca(t) + self.fr_python[sei_i, t_i] = sei.get_fr(t) self.se_python[sei_i, t_i] = sei.get_se(t) for sei_i, sei in enumerate(self.se_integrator): testing.assert_almost_equal( - self.ca_nest[n_i], self.ca_python[sei_i], decimal=5) + self.fr_nest[n_i], self.fr_python[sei_i], decimal=5) testing.assert_almost_equal( self.se_nest[n_i], self.se_python[sei_i], decimal=5) @@ -343,9 +343,9 @@ def plot(self): for i, sei in enumerate(self.se_integrator): pylab.figure() pylab.subplot(1, 2, 1) - pylab.title('Ca') - pylab.plot(self.sim_steps, self.ca_nest[0, :]) - pylab.plot(self.sim_steps, self.ca_python[i]) + pylab.title('Firing Rate in Hz') + pylab.plot(self.sim_steps, self.fr_nest[0, :]) + pylab.plot(self.sim_steps, self.fr_python[i]) pylab.legend(('nest', sei.__class__.__name__)) pylab.subplot(1, 2, 2) pylab.title('Synaptic Element') @@ -357,15 +357,15 @@ def plot(self): pylab.savefig('sp_raster_plot.png') def test_linear_growth_curve(self): - beta_ca = 0.0001 - tau_ca = 10000.0 + beta_fr = 0.1 + tau_fr = 10000.0 growth_rate = 0.0001 eps = 0.10 nest.SetStatus( self.local_nodes, { - 'beta_Ca': beta_ca, - 'tau_Ca': tau_ca, + 'beta_fr': beta_fr, + 'tau_fr': tau_fr, 'synaptic_elements': { 'se': { 'growth_curve': 'linear', @@ -377,9 +377,9 @@ def test_linear_growth_curve(self): } ) self.se_integrator.append(LinearExactSEI( - tau_ca=tau_ca, beta_ca=beta_ca, eps=eps, growth_rate=growth_rate)) + tau_fr=tau_fr, beta_fr=beta_fr, eps=eps, growth_rate=growth_rate)) self.se_integrator.append(LinearNumericSEI( - tau_ca=tau_ca, beta_ca=beta_ca, eps=eps, growth_rate=growth_rate)) + tau_fr=tau_fr, beta_fr=beta_fr, eps=eps, growth_rate=growth_rate)) self.simulate() @@ -399,16 +399,16 @@ def test_linear_growth_curve(self): decimal=8) def test_gaussian_growth_curve(self): - beta_ca = 0.0001 - tau_ca = 10000.0 + beta_fr = 0.1 + tau_fr = 10000.0 growth_rate = 0.0001 eta = 0.05 eps = 0.10 nest.SetStatus( self.local_nodes, { - 'beta_Ca': beta_ca, - 'tau_Ca': tau_ca, + 'beta_fr': beta_fr, + 'tau_fr': tau_fr, 'synaptic_elements': { 'se': { 'growth_curve': 'gaussian', @@ -419,7 +419,7 @@ def test_gaussian_growth_curve(self): } ) self.se_integrator.append( - GaussianNumericSEI(tau_ca=tau_ca, beta_ca=beta_ca, + GaussianNumericSEI(tau_fr=tau_fr, beta_fr=beta_fr, eta=eta, eps=eps, growth_rate=growth_rate)) self.simulate() @@ -439,16 +439,16 @@ def test_gaussian_growth_curve(self): decimal=5) def test_sigmoid_growth_curve(self): - beta_ca = 0.0001 - tau_ca = 10000.0 + beta_fr = 0.1 + tau_fr = 10000.0 growth_rate = 0.0001 eps = 0.10 psi = 0.10 nest.SetStatus( self.local_nodes, { - 'beta_Ca': beta_ca, - 'tau_Ca': tau_ca, + 'beta_fr': beta_fr, + 'tau_fr': tau_fr, 'synaptic_elements': { 'se': { 'growth_curve': 'sigmoid', @@ -459,7 +459,7 @@ def test_sigmoid_growth_curve(self): } ) self.se_integrator.append( - SigmoidNumericSEI(tau_ca=tau_ca, beta_ca=beta_ca, + SigmoidNumericSEI(tau_fr=tau_fr, beta_fr=beta_fr, eps=eps, psi=psi, growth_rate=growth_rate)) self.simulate() diff --git a/pynest/nest/tests/test_sp/test_update_synaptic_elements.py b/pynest/nest/tests/test_sp/test_update_synaptic_elements.py index 37f8391b7f..9367e0ac4b 100644 --- a/pynest/nest/tests/test_sp/test_update_synaptic_elements.py +++ b/pynest/nest/tests/test_sp/test_update_synaptic_elements.py @@ -32,14 +32,14 @@ def test_update_synaptic_elements(self): growth_curve_axonal = { 'growth_curve': "gaussian", 'growth_rate': 0.0001, # Beta (elements/ms) - 'eta': 0.4, - 'eps': 0.7, + 'eta': 4.0, #Hz + 'eps': 7.0, #Hz } growth_curve_dendritic_E = { 'growth_curve': "gaussian", 'growth_rate': 0.0001, # Beta (elements/ms) - 'eta': 0.1, - 'eps': 0.7, + 'eta': 1.0, #Hz + 'eps': 7.0, #Hz } structural_p_elements = { 'Den_ex': growth_curve_dendritic_E, @@ -47,8 +47,8 @@ def test_update_synaptic_elements(self): } new_growth_curve_axonal = { - 'eta': 0.0, - 'eps': 0.0, + 'eta': 0.0, #Hz + 'eps': 0.0, #Hz } elements_to_update = { 'Axon': new_growth_curve_axonal