From eee6c1c35eb828eb01cef42c826e0bc7d6a610fb Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 18 Jul 2023 15:51:06 -0700 Subject: [PATCH 01/45] Adding GLIF with double alpha psc. --- models/glif_psc_double_alpha.cpp | 711 +++++++++++++++++++++++++++++++ models/glif_psc_double_alpha.h | 494 +++++++++++++++++++++ modelsets/full | 1 + nestkernel/nest_names.cpp | 2 + nestkernel/nest_names.h | 2 + 5 files changed, 1210 insertions(+) create mode 100644 models/glif_psc_double_alpha.cpp create mode 100644 models/glif_psc_double_alpha.h diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp new file mode 100644 index 0000000000..b48b75d7b4 --- /dev/null +++ b/models/glif_psc_double_alpha.cpp @@ -0,0 +1,711 @@ +/* + * glif_psc_double_alpha.cpp + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#include "glif_psc_double_alpha.h" + +// C++ includes: +#include +#include + +// Includes from libnestutil: +#include "dict_util.h" +#include "exceptions.h" +#include "iaf_propagator.h" +#include "kernel_manager.h" +#include "universal_data_logger_impl.h" + +// Includes from sli: +#include "dict.h" +#include "dictutils.h" + +using namespace nest; + +nest::RecordablesMap< nest::glif_psc_double_alpha > nest::glif_psc_double_alpha::recordablesMap_; + +namespace nest +{ +// Override the create() method with one call to RecordablesMap::insert_() +// for each quantity to be recorded. +template <> +void +RecordablesMap< nest::glif_psc_double_alpha >::create() +{ + insert_( names::V_m, &nest::glif_psc_double_alpha::get_V_m_ ); + insert_( names::ASCurrents_sum, &nest::glif_psc_double_alpha::get_ASCurrents_sum_ ); + insert_( names::I, &nest::glif_psc_double_alpha::get_I_ ); + insert_( names::I_syn, &nest::glif_psc_double_alpha::get_I_syn_ ); + insert_( names::threshold, &nest::glif_psc_double_alpha::get_threshold_ ); + insert_( names::threshold_spike, &nest::glif_psc_double_alpha::get_threshold_spike_ ); + insert_( names::threshold_voltage, &nest::glif_psc_double_alpha::get_threshold_voltage_ ); +} +} + +/* ---------------------------------------------------------------- + * Default constructors defining default parameters and state + * ---------------------------------------------------------------- */ + +nest::glif_psc_double_alpha::Parameters_::Parameters_() + : G_( 9.43 ) // in nS + , E_L_( -78.85 ) // in mV + , th_inf_( -51.68 - E_L_ ) // in mv, rel to E_L_, - 51.68 - E_L_, i.e., 27.17 + , C_m_( 58.72 ) // in pF + , t_ref_( 3.75 ) // in ms + , V_reset_( 0.0 ) // in mV, rel to E_L_, -78.85 - E_L_ + , th_spike_add_( 0.37 ) // in mV + , th_spike_decay_( 0.009 ) // in 1/ms + , voltage_reset_fraction_( 0.20 ) + , voltage_reset_add_( 18.51 ) // in mV + , th_voltage_index_( 0.005 ) // in 1/ms + , th_voltage_decay_( 0.09 ) // in 1/ms + , asc_init_( std::vector< double >( 2, 0.0 ) ) // in pA + , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms + , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA + , asc_r_( std::vector< double >( 2, 1.0 ) ) + , tau_syn_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms + , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless + , has_connections_( false ) + , has_theta_spike_( false ) + , has_asc_( false ) + , has_theta_voltage_( false ) +{ +} + +nest::glif_psc_double_alpha::State_::State_( const Parameters_& p ) + : U_( 0.0 ) // in mV + , threshold_( p.th_inf_ ) // in mV + , threshold_spike_( 0.0 ) // in mV + , threshold_voltage_( 0.0 ) // in mV + , I_( 0.0 ) // in pA + , I_syn_( 0.0 ) // in pA + , I_syn_fast_( 0.0 ) // in pA + , I_syn_slow_( 0.0 ) // in pA + , ASCurrents_( p.asc_init_ ) // in pA + , ASCurrents_sum_( 0.0 ) // in pA + , refractory_steps_( 0 ) +{ + for ( std::size_t a = 0; a < p.asc_init_.size(); ++a ) + { + ASCurrents_sum_ += ASCurrents_[ a ]; + } + y1_.clear(); + y2_.clear(); + y1_slow_.clear(); + y2_slow_.clear(); +} + +/* ---------------------------------------------------------------- + * Parameter and state extractions and manipulation functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::Parameters_::get( DictionaryDatum& d ) const +{ + def< double >( d, names::V_th, th_inf_ + E_L_ ); + def< double >( d, names::g, G_ ); + def< double >( d, names::E_L, E_L_ ); + def< double >( d, names::C_m, C_m_ ); + def< double >( d, names::t_ref, t_ref_ ); + def< double >( d, names::V_reset, V_reset_ + E_L_ ); + + def< double >( d, names::th_spike_add, th_spike_add_ ); + def< double >( d, names::th_spike_decay, th_spike_decay_ ); + def< double >( d, names::voltage_reset_fraction, voltage_reset_fraction_ ); + def< double >( d, names::voltage_reset_add, voltage_reset_add_ ); + + def< double >( d, names::th_voltage_index, th_voltage_index_ ); + def< double >( d, names::th_voltage_decay, th_voltage_decay_ ); + + def< std::vector< double > >( d, names::asc_init, asc_init_ ); + def< std::vector< double > >( d, names::asc_decay, asc_decay_ ); + def< std::vector< double > >( d, names::asc_amps, asc_amps_ ); + def< std::vector< double > >( d, names::asc_r, asc_r_ ); + + ArrayDatum tau_syn_ad( tau_syn_ ); + def< ArrayDatum >( d, names::tau_syn, tau_syn_ad ); + + ArrayDatum tau_syn_slow_ad( tau_syn_slow_ ); + def< ArrayDatum >( d, names::tau_syn_slow, tau_syn_slow_ad ); + + ArrayDatum amp_slow_ad( amp_slow_ ); + def< ArrayDatum >( d, names::amp_slow, amp_slow_ad ); + + def< bool >( d, names::has_connections, has_connections_ ); + def< bool >( d, names::spike_dependent_threshold, has_theta_spike_ ); + def< bool >( d, names::after_spike_currents, has_asc_ ); + def< bool >( d, names::adapting_threshold, has_theta_voltage_ ); +} + +double +nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* node ) +{ + // if E_L_ is changed, we need to adjust all variables defined relative to + // E_L_ + const double ELold = E_L_; + updateValueParam< double >( d, names::E_L, E_L_, node ); + const double delta_EL = E_L_ - ELold; + + if ( updateValueParam< double >( d, names::V_reset, V_reset_, node ) ) + { + V_reset_ -= E_L_; + } + else + { + V_reset_ -= delta_EL; + } + + if ( updateValueParam< double >( d, names::V_th, th_inf_, node ) ) + { + th_inf_ -= E_L_; + } + else + { + th_inf_ -= delta_EL; + } + + updateValueParam< double >( d, names::g, G_, node ); + updateValueParam< double >( d, names::C_m, C_m_, node ); + updateValueParam< double >( d, names::t_ref, t_ref_, node ); + + updateValueParam< double >( d, names::th_spike_add, th_spike_add_, node ); + updateValueParam< double >( d, names::th_spike_decay, th_spike_decay_, node ); + updateValueParam< double >( d, names::voltage_reset_fraction, voltage_reset_fraction_, node ); + updateValueParam< double >( d, names::voltage_reset_add, voltage_reset_add_, node ); + + updateValueParam< double >( d, names::th_voltage_index, th_voltage_index_, node ); + updateValueParam< double >( d, names::th_voltage_decay, th_voltage_decay_, node ); + + updateValue< std::vector< double > >( d, names::asc_init, asc_init_ ); + updateValue< std::vector< double > >( d, names::asc_decay, asc_decay_ ); + updateValue< std::vector< double > >( d, names::asc_amps, asc_amps_ ); + updateValue< std::vector< double > >( d, names::asc_r, asc_r_ ); + + // set model mechanisms + updateValueParam< bool >( d, names::spike_dependent_threshold, has_theta_spike_, node ); + updateValueParam< bool >( d, names::after_spike_currents, has_asc_, node ); + updateValueParam< bool >( d, names::adapting_threshold, has_theta_voltage_, node ); + + // check model mechanisms parameter + if ( not( ( not has_theta_spike_ and not has_asc_ and not has_theta_voltage_ ) or // glif1 + ( has_theta_spike_ and not has_asc_ and not has_theta_voltage_ ) or // glif2 + ( not has_theta_spike_ and has_asc_ and not has_theta_voltage_ ) or // glif3 + ( has_theta_spike_ and has_asc_ and not has_theta_voltage_ ) or // glif4 + ( has_theta_spike_ and has_asc_ and has_theta_voltage_ ) // glif5 + ) ) + { + throw BadProperty( + "Incorrect model mechanism combination setting." + "See documentation for setting of model mechanism parameters:" + "spike_dependent_threshold, after_spike_currents, adapting_threshold." ); + } + + // check after ASC parameters' sizes and values + if ( has_asc_ ) + { + // check size + const size_t asc_size = asc_decay_.size(); + if ( not( + ( asc_init_.size() == asc_size ) and ( asc_amps_.size() == asc_size ) and ( asc_r_.size() == asc_size ) ) ) + { + throw BadProperty( + "All after spike current parameters (i.e., asc_init, k, asc_amps, r) " + "must have the same size." ); + } + + // check values + for ( std::size_t a = 0; a < asc_decay_.size(); ++a ) + { + if ( asc_decay_[ a ] <= 0.0 ) + { + throw BadProperty( "After-spike current time constant must be strictly positive." ); + } + + if ( asc_r_[ a ] < 0.0 or asc_r_[ a ] > 1.0 ) + { + throw BadProperty( "After spike current fraction following spike coefficients r must be within [0.0, 1.0]." ); + } + } + } + + if ( V_reset_ >= th_inf_ ) + { + throw BadProperty( "Reset potential must be smaller than threshold." ); + } + + if ( C_m_ <= 0.0 ) + { + throw BadProperty( "Capacitance must be strictly positive." ); + } + + if ( G_ <= 0.0 ) + { + throw BadProperty( "Membrane conductance must be strictly positive." ); + } + + if ( t_ref_ <= 0.0 ) + { + throw BadProperty( "Refractory time constant must be strictly positive." ); + } + + if ( has_theta_voltage_ ) + { + if ( th_voltage_decay_ <= 0.0 ) + { + throw BadProperty( "Voltage-induced threshold time constant must be strictly positive." ); + } + } + + // check spike component parameters + if ( has_theta_spike_ ) + { + if ( th_spike_decay_ <= 0.0 ) + { + throw BadProperty( "Spike induced threshold time constant must be strictly positive." ); + } + + if ( voltage_reset_fraction_ < 0.0 or voltage_reset_fraction_ > 1.0 ) + { + throw BadProperty( "Voltage fraction coefficient following spike must be within [0.0, 1.0]." ); + } + } + + const size_t old_n_receptors = this->n_receptors_(); + if ( updateValue< std::vector< double > >( d, names::tau_syn, tau_syn_ ) ) + { + if ( this->n_receptors_() != old_n_receptors and has_connections_ ) + { + throw BadProperty( + "The neuron has connections, therefore the number of ports cannot be " + "reduced." ); + } + for ( size_t i = 0; i < tau_syn_.size(); ++i ) + { + if ( tau_syn_[ i ] <= 0 ) + { + throw BadProperty( "All synaptic time constants must be strictly positive." ); + } + } + } + if ( updateValue< std::vector< double > >( d, names::tau_syn_slow, tau_syn_slow_ ) ) + { + for ( size_t i = 0; i < tau_syn_slow_.size(); ++i ) + { + if ( tau_syn_slow_[ i ] <= 0 ) + { + throw BadProperty( "All slow synaptic time constants must be strictly positive." ); + } + } + } + // for amp_slow + if ( updateValue< std::vector< double > >( d, names::amp_slow, amp_slow_ ) ) + { + for ( size_t i = 0; i < amp_slow_.size(); ++i ) + { + if ( amp_slow_[ i ] <= 0 ) + { + throw BadProperty( "All slow synaptic amplitudes must be strictly positive." ); + } + } + } + + return delta_EL; +} + +void +nest::glif_psc_double_alpha::State_::get( DictionaryDatum& d, const Parameters_& p ) const +{ + def< double >( d, names::V_m, U_ + p.E_L_ ); + def< std::vector< double > >( d, names::ASCurrents, ASCurrents_ ); + def< double >( d, names::threshold_spike, threshold_spike_ ); + def< double >( d, names::threshold_voltage, threshold_voltage_ ); +} + +void +nest::glif_psc_double_alpha::State_::set( const DictionaryDatum& d, const Parameters_& p, double delta_EL, Node* node ) +{ + if ( updateValueParam< double >( d, names::V_m, U_, node ) ) + { + U_ -= p.E_L_; + } + else + { + U_ -= delta_EL; + } + + bool asc_flag = updateValue< std::vector< double > >( d, names::ASCurrents, ASCurrents_ ); + if ( asc_flag and not p.has_asc_ ) + { + throw BadProperty( "After spike currents are not supported or settable in the current model mechanisms." ); + } + + const size_t asc_size = p.asc_decay_.size(); + if ( asc_flag ) + { + if ( ASCurrents_.size() != asc_size ) + { + throw BadProperty( "After spike current values must have have the same size (" + std::to_string( asc_size ) + + ") of its parameters (i.e., asc_init, k, asc_amps, r)." ); + } + } + + if ( updateValueParam< double >( d, names::threshold_spike, threshold_spike_, node ) and not p.has_theta_spike_ ) + { + throw BadProperty( "Threshold spike component is not supported or settable in the current model mechanisms." ); + } + + if ( updateValueParam< double >( d, names::threshold_voltage, threshold_voltage_, node ) + and not p.has_theta_voltage_ ) + { + throw BadProperty( "Threshold voltage component is not supported or settable in the current model mechanisms." ); + } +} + +nest::glif_psc_double_alpha::Buffers_::Buffers_( glif_psc_double_alpha& n ) + : logger_( n ) +{ +} + +nest::glif_psc_double_alpha::Buffers_::Buffers_( const Buffers_&, glif_psc_double_alpha& n ) + : logger_( n ) +{ +} + +/* ---------------------------------------------------------------- + * Default and copy constructor for node + * ---------------------------------------------------------------- */ + +nest::glif_psc_double_alpha::glif_psc_double_alpha() + : ArchivingNode() + , P_() + , S_( P_ ) + , B_( *this ) +{ + recordablesMap_.create(); +} + +nest::glif_psc_double_alpha::glif_psc_double_alpha( const glif_psc_double_alpha& n ) + : ArchivingNode( n ) + , P_( n.P_ ) + , S_( n.S_ ) + , B_( n.B_, *this ) +{ +} + +/* ---------------------------------------------------------------- + * Node initialization functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::init_buffers_() +{ + B_.spikes_.clear(); // includes resize + B_.currents_.clear(); // include resize + B_.logger_.reset(); // includes resize +} + +void +nest::glif_psc_double_alpha::pre_run_hook() +{ + B_.logger_.init(); + + const double h = Time::get_resolution().get_ms(); // in ms + + // pre-computing of decay parameters + if ( P_.has_theta_spike_ ) + { + V_.theta_spike_decay_rate_ = std::exp( -P_.th_spike_decay_ * h ); + V_.theta_spike_refractory_decay_rate_ = std::exp( -P_.th_spike_decay_ * P_.t_ref_ ); + } + + if ( P_.has_asc_ ) + { + V_.asc_decay_rates_.resize( P_.asc_decay_.size() ); + V_.asc_stable_coeff_.resize( P_.asc_decay_.size() ); + V_.asc_refractory_decay_rates_.resize( P_.asc_decay_.size() ); + for ( std::size_t a = 0; a < P_.asc_decay_.size(); ++a ) + { + V_.asc_decay_rates_[ a ] = std::exp( -P_.asc_decay_[ a ] * h ); + V_.asc_stable_coeff_[ a ] = ( ( 1.0 / P_.asc_decay_[ a ] ) / h ) * ( 1.0 - V_.asc_decay_rates_[ a ] ); + V_.asc_refractory_decay_rates_[ a ] = P_.asc_r_[ a ] * std::exp( -P_.asc_decay_[ a ] * P_.t_ref_ ); + } + } + + if ( P_.has_theta_voltage_ ) + { + V_.potential_decay_rate_ = std::exp( -P_.G_ * h / P_.C_m_ ); + V_.theta_voltage_decay_rate_inverse_ = 1 / ( std::exp( P_.th_voltage_decay_ * h ) ); + V_.phi = P_.th_voltage_index_ / ( P_.th_voltage_decay_ - P_.G_ / P_.C_m_ ); + V_.abpara_ratio_voltage_ = P_.th_voltage_index_ / P_.th_voltage_decay_; + } + + // postsynaptic currents + V_.P11_.resize( P_.n_receptors_() ); + V_.P21_.resize( P_.n_receptors_() ); + V_.P22_.resize( P_.n_receptors_() ); + V_.P31_.resize( P_.n_receptors_() ); + V_.P32_.resize( P_.n_receptors_() ); + + // Additional variables for the slow component + V_.P11_slow_.resize( P_.n_receptors_() ); + V_.P21_slow_.resize( P_.n_receptors_() ); + V_.P22_slow_.resize( P_.n_receptors_() ); + V_.P31_slow_.resize( P_.n_receptors_() ); + V_.P32_slow_.resize( P_.n_receptors_() ); + + S_.y1_.resize( P_.n_receptors_() ); + S_.y2_.resize( P_.n_receptors_() ); + S_.y1_slow_.resize( P_.n_receptors_() ); // Slow component + S_.y2_slow_.resize( P_.n_receptors_() ); // Slow component + V_.PSCInitialValues_.resize( P_.n_receptors_() ); + V_.PSCInitialValues_slow_.resize( P_.n_receptors_() ); // Slow component + + B_.spikes_.resize( P_.n_receptors_() ); + + double Tau_ = P_.C_m_ / P_.G_; // in ms + V_.P33_ = std::exp( -h / Tau_ ); + V_.P30_ = 1 / P_.C_m_ * ( 1 - V_.P33_ ) * Tau_; + + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + // these P are independent + V_.P11_[ i ] = V_.P22_[ i ] = std::exp( -h / P_.tau_syn_[ i ] ); + + V_.P21_[ i ] = h * V_.P11_[ i ]; + + // these are determined according to a numeric stability criterion + // input time parameter shall be in ms, capacity in pF + std::tie( V_.P31_[ i ], V_.P32_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + + // Slow component + V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); + V_.P21_slow_[ i ] = h * V_.P11_slow_[ i ]; + std::tie( V_.P31_slow_[ i ], V_.P32_slow_[ i ] ) = + IAFPropagatorAlpha( P_.tau_syn_slow_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + + + V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_[ i ]; + V_.PSCInitialValues_slow_[ i ] = 1.0 * numerics::e / P_.tau_syn_slow_[ i ] * P_.amp_slow_[ i ]; + B_.spikes_[ i ].resize(); + } + + V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); +} + +/* ---------------------------------------------------------------- + * Update and spike handling functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::update( Time const& origin, const long from, const long to ) +{ + + double v_old = S_.U_; + + for ( long lag = from; lag < to; ++lag ) + { + + if ( S_.refractory_steps_ == 0 ) + { + // neuron not refractory, integrate voltage and currents + + // update threshold via exact solution of dynamics of spike component of + // threshold for glif2/4/5 models with "R" + if ( P_.has_theta_spike_ ) + { + S_.threshold_spike_ = S_.threshold_spike_ * V_.theta_spike_decay_rate_; + } + + // Calculate new ASCurrents value using exponential methods + S_.ASCurrents_sum_ = 0.0; + // for glif3/4/5 models with "ASC" + // take after spike current value at the beginning of the time to compute + // the exact mean ASC for the time step and sum the exact ASCs of all ports; + // and then update the current values to the value at the end of the time + // step, ready for the next time step + if ( P_.has_asc_ ) + { + for ( std::size_t a = 0; a < S_.ASCurrents_.size(); ++a ) + { + S_.ASCurrents_sum_ += ( V_.asc_stable_coeff_[ a ] * S_.ASCurrents_[ a ] ); + S_.ASCurrents_[ a ] = S_.ASCurrents_[ a ] * V_.asc_decay_rates_[ a ]; + } + } + + // voltage dynamics of membranes, linear exact to find next V_m value + S_.U_ = v_old * V_.P33_ + ( S_.I_ + S_.ASCurrents_sum_ ) * V_.P30_; + + // add synapse component for voltage dynamics + S_.I_syn_fast_ = 0.0; + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + S_.U_ += V_.P31_[ i ] * S_.y1_[ i ] + V_.P32_[ i ] * S_.y2_[ i ]; + S_.I_syn_fast_ += S_.y2_[ i ]; + } + // slow component + S_.I_syn_slow_ = 0.0; + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + // S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + (V_.P32_slow_[ i ] * S_.y2_slow_[ i ]) * P_.amp_slow_[ i ]; + S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + V_.P32_slow_[ i ] * S_.y2_slow_[ i ]; + // multiply amp_slow_ + S_.I_syn_slow_ += S_.y2_slow_[ i ]; + // print all components of P3x and y the slow current. + // std::cout << "P31_slow: " << V_.P31_slow_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << " y2_slow:" << + // S_.y2_slow_[ i ] << std::endl; + } + // add the fast and slow components + S_.I_syn_ = S_.I_syn_fast_ + S_.I_syn_slow_; + // print to diagnose the synaptic current + // std::cout << "I_syn_fast: " << S_.I_syn_fast_ << " I_syn_slow: " << S_.I_syn_slow_ << std::endl; + + // Calculate exact voltage component of the threshold for glif5 model with + // "A" + if ( P_.has_theta_voltage_ ) + { + const double beta = ( S_.I_ + S_.ASCurrents_sum_ ) / P_.G_; + S_.threshold_voltage_ = V_.phi * ( v_old - beta ) * V_.potential_decay_rate_ + + V_.theta_voltage_decay_rate_inverse_ + * ( S_.threshold_voltage_ - V_.phi * ( v_old - beta ) - V_.abpara_ratio_voltage_ * beta ) + + V_.abpara_ratio_voltage_ * beta; + } + + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + + // Check if there is an action potential + if ( S_.U_ > S_.threshold_ ) + { + // Marks that the neuron is in a refractory period + S_.refractory_steps_ = V_.RefractoryCounts_; + + // Reset ASC_currents for glif3/4/5 models with "ASC" + if ( P_.has_asc_ ) + { + for ( std::size_t a = 0; a < S_.ASCurrents_.size(); ++a ) + { + S_.ASCurrents_[ a ] = P_.asc_amps_[ a ] + S_.ASCurrents_[ a ] * V_.asc_refractory_decay_rates_[ a ]; + } + } + + // Reset voltage + if ( not P_.has_theta_spike_ ) + { + // Reset voltage for glif1/3 models without "R" + S_.U_ = P_.V_reset_; + } + else + { + // Reset voltage for glif2/4/5 models with "R" + S_.U_ = P_.voltage_reset_fraction_ * v_old + P_.voltage_reset_add_; + + // reset spike component of threshold + // (decay for refractory period and then add additive constant) + S_.threshold_spike_ = S_.threshold_spike_ * V_.theta_spike_refractory_decay_rate_ + P_.th_spike_add_; + + // rest the global threshold (voltage component of threshold: stay the + // same) + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + } + + set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); + SpikeEvent se; + kernel().event_delivery_manager.send( *this, se, lag ); + } + } + else + { + // neuron is absolute refractory + --S_.refractory_steps_; + + // While neuron is in refractory period count-down in time steps (since dt + // may change while in refractory) while holding the voltage at last peak. + S_.U_ = v_old; + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + } + + // alpha shape PSCs + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + double spike_value = B_.spikes_[ i ].get_value( lag ); + // Apply spikes delivered in this step: The spikes arriving at T+1 have an + // immediate effect on the state of the neuron + // store B_.spikes_[ i ].get_value( lag ); into a variable + S_.y2_[ i ] = V_.P21_[ i ] * S_.y1_[ i ] + V_.P22_[ i ] * S_.y2_[ i ]; + S_.y1_[ i ] *= V_.P11_[ i ]; + S_.y1_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; + + // slow component + S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; + S_.y1_slow_[ i ] *= V_.P11_slow_[ i ]; + S_.y1_slow_[ i ] += V_.PSCInitialValues_slow_[ i ] * spike_value; + // print y1 and y1 slow + // std::cout << "y1: " << S_.y1_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << std::endl; + // print PSCInitlaValues + // std::cout << "PSCInitialValues: " << V_.PSCInitialValues_[ i ] << " PSCInitialValues_slow: " << + // V_.PSCInitialValues_slow_[ i ] << std::endl; + } + + // Update any external currents + S_.I_ = B_.currents_.get_value( lag ); + + // Save voltage + B_.logger_.record_data( origin.get_steps() + lag ); + v_old = S_.U_; + } +} + +size_t +nest::glif_psc_double_alpha::handles_test_event( SpikeEvent&, size_t receptor_type ) +{ + if ( receptor_type <= 0 or receptor_type > P_.n_receptors_() ) + { + throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" ); + } + + P_.has_connections_ = true; + return receptor_type; +} + +void +nest::glif_psc_double_alpha::handle( SpikeEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + B_.spikes_[ e.get_rport() - 1 ].add_value( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), e.get_weight() * e.get_multiplicity() ); +} + +void +nest::glif_psc_double_alpha::handle( CurrentEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + B_.currents_.add_value( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), e.get_weight() * e.get_current() ); +} + +// Do not move this function as inline to h-file. It depends on +// universal_data_logger_impl.h being included here. +void +nest::glif_psc_double_alpha::handle( DataLoggingRequest& e ) +{ + B_.logger_.handle( e ); // the logger does this for us +} diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h new file mode 100644 index 0000000000..2b6336fae3 --- /dev/null +++ b/models/glif_psc_double_alpha.h @@ -0,0 +1,494 @@ +/* + * glif_psc_double_alpha.h + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#ifndef GLIF_PSC_DOUBLE_ALPHA_H +#define GLIF_PSC_DOUBLE_ALPHA_H + +#include "archiving_node.h" +#include "connection.h" +#include "event.h" +#include "nest_types.h" +#include "ring_buffer.h" +#include "universal_data_logger.h" + +#include "dictdatum.h" + +/* BeginUserDocs: integrate-and-fire, current-based + +Short description ++++++++++++++++++ + +Current-based generalized leaky integrate-and-fire (GLIF) models (from the Allen Institute) + +Description ++++++++++++ + +``glif_psc_double_alpha`` provides five generalized leaky integrate-and-fire +(GLIF) models [1]_ with alpha-function shaped synaptic currents. +Incoming spike events induce a postsynaptic change of current modeled +by sum of two alpha functions [2]_. The alpha function is normalized such that an event +of weight 1.0 results in a peak current of the fast component of the alpha function to +be 1 pA at :math:`t = tau_syn`. The peak current of the slow component of the alpha is +given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha has +a single synapse that is accessible through receptor_port 1. An arbitrary number of +synapses with different time constants can be configured by setting the desired time +constants as tau_syn array. The resulting synapses are addressed through receptor_port +1, 2, 3, .... + +The five GLIF models are: + +* **GLIF Model 1** - Traditional leaky integrate and fire (LIF) +* **GLIF Model 2** - Leaky integrate and fire with biologically defined reset rules + (LIF_R) +* **GLIF Model 3** - Leaky integrate and fire with after-spike currents (LIF_ASC) +* **GLIF Model 4** - Leaky integrate and fire with biologically defined reset rules + and after-spike currents (LIF_R_ASC) +* **GLIF Model 5** - Leaky integrate and fire with biologically defined reset rules, + after-spike currents and a voltage dependent threshold (LIF_R_ASC_A) + +GLIF model mechanism setting is based on three parameters +(``spike_dependent_threshold``, ``after_spike_currents``, ``adapting_threshold``). +The settings of these three parameters for the five GLIF models are listed +below. Other combinations of these parameters will not be supported. + ++--------+---------------------------+----------------------+--------------------+ +| Model | spike_dependent_threshold | after_spike_currents | adapting_threshold | ++========+===========================+======================+====================+ +| GLIF1 | False | False | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF2 | True | False | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF3 | False | True | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF4 | True | True | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF5 | True | True | True | ++--------+---------------------------+----------------------+--------------------+ + +Typical parameter setting of different levels of GLIF models for different cells +can be found and downloaded in the `Allen Cell Type Database +`_. For example, the default parameter setting of this +``glif_psc_double_alpha`` neuron model was from the parameter values of GLIF Model 5 of +Cell 490626718, which can be retrieved from the `Allen Brain Atlas +`_, with units being converted from SI units (i.e., V, S (1/Ohm), +F, s, A) to NEST used units (i.e., mV, nS (1/GOhm), pF, ms, pA) and values +being rounded to appropriate digits for simplification. + +For models with spike dependent threshold (i.e., GLIF2, GLIF4 and GLIF5), +parameter setting of voltage_reset_fraction and voltage_reset_add may lead to the +situation that voltage is bigger than threshold after reset. In this case, the neuron +will continue to spike until the end of the simulation regardless the stimulated inputs. +We recommend the setting of the parameters of these three models to follow the +condition of + +.. math:: + + E_L + \mathrm{voltage\_reset\_fraction} \cdot \left( V_\mathrm{th} - E_L \right) + + \mathrm{voltage\_reset\_add} < V_\mathrm{th} + \mathrm{th\_spike\_add} + +.. note:: + + If ``tau_m`` is very close to ``tau_syn_ex`` or ``tau_syn_in``, the model + will numerically behave as if ``tau_m`` is equal to ``tau_syn_ex`` or + ``tau_syn_in``, respectively, to avoid numerical instabilities. + + For implementation details see the + `IAF_neurons_singularity <../model_details/IAF_neurons_singularity.ipynb>`_ notebook. + +Parameters +++++++++++ + +The following parameters can be set in the status dictionary. + +========= ======== ============================================================ +**Membrane parameters** +------------------------------------------------------------------------------- +V_m double Membrane potential in mV (absolute value) +V_th double Instantaneous threshold in mV +g double Membrane conductance in nS +E_L double Resting membrane potential in mV +C_m double Capacitance of the membrane in pF +t_ref double Duration of refractory time in ms +V_reset double Reset potential of the membrane in mV (GLIF 1 or GLIF 3) +========= ======== ============================================================ + +========================= =============== ===================================== +**Spike adaptation and firing intensity parameters** +------------------------------------------------------------------------------- +th_spike_add double Threshold addition following spike + in mV (delta_theta_s in Equation (6) + in [1]_) +th_spike_decay double Spike-induced threshold time + constant in 1/ms (bs in Equation (2) + in [1]_) +voltage_reset_fraction double Voltage fraction coefficient + following spike (fv in Equation (5) + in [1]_) +voltage_reset_add double Voltage addition following spike in + mV (-delta_V (sign flipped) in + Equation (5) in [1]_) +asc_init double vector Initial values of after-spike + currents in pA +asc_decay double vector After-spike current time constants + in 1/ms (kj in Equation (3) in [1]_) +asc_amps double vector After-spike current amplitudes in + pA (deltaIj in Equation (7) in [1]_) +asc_r double vector Current fraction following spike + coefficients for fj in Equation (7) + in [1]_ +th_voltage_index double Adaptation index of threshold - A + 'leak-conductance' for the + voltage-dependent component of the + threshold in 1/ms (av in Equation + (4) in [1]_) +th_voltage_decay double Voltage-induced threshold time + constant - Inverse of which is the + time constant of the + voltage-dependent component of the + threshold in 1/ms (bv in Equation + (4) in [1]_) +tau_syn double vector Rise time constants of the synaptic + alpha function in ms +tau_syn_slow double vector Rise time constants of the slower + synaptic alpha function in ms +amp_slow double vector Relative amplitude of the slower + synaptic alpha function +E_rev double vector Reversal potential in mV +spike_dependent_threshold bool flag whether the neuron has + biologically defined reset rules + with a spike dependent threshold + component +after_spike_currents bool flag whether the neuron has after + spike currents +adapting_threshold bool flag whether the neuron has a + voltage dependent threshold component +========================= =============== ===================================== + +References +++++++++++ + +.. [1] Teeter C, Iyer R, Menon V, Gouwens N, Feng D, Berg J, Szafer A, + Cain N, Zeng H, Hawrylycz M, Koch C, & Mihalas S (2018) + Generalized leaky integrate-and-fire models classify multiple neuron + types. Nature Communications 9:709. +.. [2] Meffin, H., Burkitt, A. N., & Grayden, D. B. (2004). An analytical + model for the large, fluctuating synaptic conductance state typical of + neocortical neurons in vivo. J. Comput. Neurosci., 16, 159-175. + +See also +++++++++ + +gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp, +glif_psc + +EndUserDocs */ + +namespace nest +{ + +class glif_psc_double_alpha : public ArchivingNode +{ +public: + glif_psc_double_alpha(); + + glif_psc_double_alpha( const glif_psc_double_alpha& ); + + using nest::Node::handle; + using nest::Node::handles_test_event; + + size_t send_test_event( nest::Node&, size_t, nest::synindex, bool ) override; + + void handle( nest::SpikeEvent& ) override; + void handle( nest::CurrentEvent& ) override; + void handle( nest::DataLoggingRequest& ) override; + + size_t handles_test_event( nest::SpikeEvent&, size_t ) override; + size_t handles_test_event( nest::CurrentEvent&, size_t ) override; + size_t handles_test_event( nest::DataLoggingRequest&, size_t ) override; + + void get_status( DictionaryDatum& ) const override; + void set_status( const DictionaryDatum& ) override; + +private: + //! Reset internal buffers of neuron. + void init_buffers_() override; + + //! Initialize auxiliary quantities, leave parameters and state untouched. + void pre_run_hook() override; + + //! Take neuron through given time interval + void update( nest::Time const&, const long, const long ) override; + + // The next two classes need to be friends to access the State_ class/member + friend class nest::RecordablesMap< glif_psc_double_alpha >; + friend class nest::UniversalDataLogger< glif_psc_double_alpha >; + + struct Parameters_ + { + double G_; //!< membrane conductance in nS + double E_L_; //!< resting potential in mV + double th_inf_; //!< infinity threshold in mV + double C_m_; //!< capacitance in pF + double t_ref_; //!< refractory time in ms + double V_reset_; //!< Membrane voltage following spike in mV + double th_spike_add_; //!< threshold additive constant following reset in mV + double th_spike_decay_; //!< spike induced threshold in 1/ms + double voltage_reset_fraction_; //!< voltage fraction following reset coefficient + double voltage_reset_add_; //!< voltage additive constant following reset in mV + double th_voltage_index_; //!< a 'leak-conductance' for the voltage-dependent + //!< component of the threshold in 1/ms + double th_voltage_decay_; //!< inverse of which is the time constant of the + //!< voltage-dependent component of the threshold in 1/ms + std::vector< double > asc_init_; //!< initial values of ASCurrents_ in pA + std::vector< double > asc_decay_; //!< predefined time scale in 1/ms + std::vector< double > asc_amps_; //!< in pA + std::vector< double > asc_r_; //!< coefficient + std::vector< double > tau_syn_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms + std::vector< double > amp_slow_; //!< synaptic port time constants in ms + + //! boolean flag which indicates whether the neuron has connections + bool has_connections_; + + //! boolean flag which indicates whether the neuron has spike dependent threshold component + bool has_theta_spike_; + + //! boolean flag which indicates whether the neuron has after spike currents + bool has_asc_; + + //! boolean flag which indicates whether the neuron has voltage dependent threshold component + bool has_theta_voltage_; + + size_t n_receptors_() const; //!< Returns the size of tau_syn_ + + Parameters_(); + + void get( DictionaryDatum& ) const; + double set( const DictionaryDatum&, Node* ); + }; + + struct State_ + { + double U_; //!< relative membrane potential in mV + double threshold_; //!< total threshold in mV + double threshold_spike_; //!< spike component of threshold in mV + double threshold_voltage_; //!< voltage component of threshold in mV + double I_; //!< external current in pA + double I_syn_; //!< postsynaptic current in pA + double I_syn_fast_; //!< postsynaptic current in pA + double I_syn_slow_; //!< postsynaptic current in pA + std::vector< double > ASCurrents_; //!< after-spike currents in pA + double ASCurrents_sum_; //!< in pA + int refractory_steps_; //!< Number of refractory steps remaining + std::vector< double > y1_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA + + State_( const Parameters_& ); + + void get( DictionaryDatum&, const Parameters_& ) const; + void set( const DictionaryDatum&, const Parameters_&, double, Node* ); + }; + + + struct Buffers_ + { + Buffers_( glif_psc_double_alpha& ); + Buffers_( const Buffers_&, glif_psc_double_alpha& ); + + std::vector< nest::RingBuffer > spikes_; //!< Buffer incoming spikes through delay, as sum + nest::RingBuffer currents_; //!< Buffer incoming currents through delay, + + //! Logger for all analog data + nest::UniversalDataLogger< glif_psc_double_alpha > logger_; + }; + + struct Variables_ + { + int RefractoryCounts_; //!< counter during refractory period + double theta_spike_decay_rate_; //!< threshold spike component decay rate + double theta_spike_refractory_decay_rate_; //!< threshold spike component decay rate during refractory + double theta_voltage_decay_rate_inverse_; //!< inverse of threshold voltage component decay rate + double potential_decay_rate_; //!< membrane potential decay rate + double abpara_ratio_voltage_; //!< ratio of parameters of voltage threshold component av/bv + std::vector< double > asc_decay_rates_; //!< after spike current decay rates + std::vector< double > asc_stable_coeff_; //!< after spike current stable coefficient + std::vector< double > asc_refractory_decay_rates_; //!< after spike current decay rates during refractory + double phi; //!< threshold voltage component coefficient + + std::vector< double > P11_; //!< synaptic current evolution parameter + std::vector< double > P21_; //!< synaptic current evolution parameter + std::vector< double > P22_; //!< synaptic current evolution parameter + double P30_; //!< membrane current/voltage evolution parameter + double P33_; //!< membrane voltage evolution parameter + std::vector< double > P31_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_; //!< synaptic/membrane current evolution parameter + + // slow component + std::vector< double > P11_slow_; //!< synaptic current evolution parameter + std::vector< double > P21_slow_; //!< synaptic current evolution parameter + std::vector< double > P22_slow_; //!< synaptic current evolution parameter + double P30_slow_; //!< membrane current/voltage evolution parameter + double P33_slow_; //!< membrane voltage evolution parameter + std::vector< double > P31_slow_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_slow_; //!< synaptic/membrane current evolution parameter + + + /** Amplitude of the synaptic current. + This value is chosen such that a postsynaptic current with + weight one has an amplitude of 1 pA. + */ + std::vector< double > PSCInitialValues_; + std::vector< double > PSCInitialValues_slow_; + }; + + double + get_V_m_() const + { + return S_.U_ + P_.E_L_; + } + + double + get_ASCurrents_sum_() const + { + return S_.ASCurrents_sum_; + } + + double + get_I_() const + { + return S_.I_; + } + + double + get_I_syn_() const + { + return S_.I_syn_; + } + + // double + // get_I_syn_fast_() const + // { + // return S_.I_syn_fast_; + // } + + // double + // get_I_syn_slow_() const + // { + // return S_.I_syn_slow_; + // } + + double + get_threshold_() const + { + return S_.threshold_ + P_.E_L_; + } + + double + get_threshold_spike_() const + { + return S_.threshold_spike_; + } + + double + get_threshold_voltage_() const + { + return S_.threshold_voltage_; + } + + Parameters_ P_; + State_ S_; + Variables_ V_; + Buffers_ B_; + + //! Mapping of recordables names to access functions + static nest::RecordablesMap< glif_psc_double_alpha > recordablesMap_; +}; + + +inline size_t +nest::glif_psc_double_alpha::Parameters_::n_receptors_() const +{ + return tau_syn_.size(); +} + +inline size_t +nest::glif_psc_double_alpha::send_test_event( nest::Node& target, size_t receptor_type, nest::synindex, bool ) +{ + nest::SpikeEvent e; + e.set_sender( *this ); + return target.handles_test_event( e, receptor_type ); +} + +inline size_t +nest::glif_psc_double_alpha::handles_test_event( nest::CurrentEvent&, size_t receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw nest::UnknownReceptorType( receptor_type, get_name() ); + } + return 0; +} + +inline size_t +nest::glif_psc_double_alpha::handles_test_event( nest::DataLoggingRequest& dlr, size_t receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw nest::UnknownReceptorType( receptor_type, get_name() ); + } + return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); +} + +inline void +glif_psc_double_alpha::get_status( DictionaryDatum& d ) const +{ + // get our own parameter and state data + P_.get( d ); + S_.get( d, P_ ); + + // get information managed by parent class + ArchivingNode::get_status( d ); + + ( *d )[ nest::names::recordables ] = recordablesMap_.get_list(); +} + +inline void +glif_psc_double_alpha::set_status( const DictionaryDatum& d ) +{ + Parameters_ ptmp = P_; // temporary copy in case of errors + const double delta_EL = ptmp.set( d, this ); // throws if BadProperty + State_ stmp = S_; // temporary copy in case of errors + stmp.set( d, ptmp, delta_EL, this ); // throws if BadProperty + + ArchivingNode::set_status( d ); + + // if we get here, temporaries contain consistent set of properties + P_ = ptmp; + S_ = stmp; +} + +} // namespace nest + +#endif diff --git a/modelsets/full b/modelsets/full index 4f0e835f41..9f7ff8bbed 100644 --- a/modelsets/full +++ b/modelsets/full @@ -32,6 +32,7 @@ gif_pop_psc_exp ginzburg_neuron glif_cond glif_psc +glif_psc_double_alpha hh_cond_exp_traub hh_cond_beta_gap_traub hh_psc_alpha diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 1eed7b4866..782e1f1f9f 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -62,6 +62,7 @@ const Name allow_oversized_mask( "allow_oversized_mask" ); const Name alpha( "alpha" ); const Name alpha_1( "alpha_1" ); const Name alpha_2( "alpha_2" ); +const Name amp_slow( "amp_slow" ); const Name amplitude( "amplitude" ); const Name amplitude_times( "amplitude_times" ); const Name amplitude_values( "amplitude_values" ); @@ -502,6 +503,7 @@ const Name tau_stc( "tau_stc" ); const Name tau_syn( "tau_syn" ); const Name tau_syn_ex( "tau_syn_ex" ); const Name tau_syn_in( "tau_syn_in" ); +const Name tau_syn_slow( "tau_syn_slow" ); const Name tau_theta( "tau_theta" ); const Name tau_u_bar_bar( "tau_u_bar_bar" ); const Name tau_u_bar_minus( "tau_u_bar_minus" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index 03caa15f0c..8a4ca4f1d0 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -87,6 +87,7 @@ extern const Name allow_oversized_mask; extern const Name alpha; extern const Name alpha_1; extern const Name alpha_2; +extern const Name amp_slow; extern const Name amplitude; extern const Name amplitude_times; extern const Name amplitude_values; @@ -527,6 +528,7 @@ extern const Name tau_stc; extern const Name tau_syn; extern const Name tau_syn_ex; extern const Name tau_syn_in; +extern const Name tau_syn_slow; extern const Name tau_theta; extern const Name tau_u_bar_bar; extern const Name tau_u_bar_minus; From cd5698a1365b3ba622459f2adea3e54d5fedc830 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 1 Aug 2023 10:48:57 -0700 Subject: [PATCH 02/45] Cleaning up commented out codes (They were used during development.) --- models/glif_psc_double_alpha.cpp | 12 ------------ models/glif_psc_double_alpha.h | 12 ------------ 2 files changed, 24 deletions(-) diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index b48b75d7b4..10011da27b 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -564,18 +564,11 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.I_syn_slow_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { - // S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + (V_.P32_slow_[ i ] * S_.y2_slow_[ i ]) * P_.amp_slow_[ i ]; S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + V_.P32_slow_[ i ] * S_.y2_slow_[ i ]; - // multiply amp_slow_ S_.I_syn_slow_ += S_.y2_slow_[ i ]; - // print all components of P3x and y the slow current. - // std::cout << "P31_slow: " << V_.P31_slow_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << " y2_slow:" << - // S_.y2_slow_[ i ] << std::endl; } // add the fast and slow components S_.I_syn_ = S_.I_syn_fast_ + S_.I_syn_slow_; - // print to diagnose the synaptic current - // std::cout << "I_syn_fast: " << S_.I_syn_fast_ << " I_syn_slow: " << S_.I_syn_slow_ << std::endl; // Calculate exact voltage component of the threshold for glif5 model with // "A" @@ -656,11 +649,6 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; S_.y1_slow_[ i ] *= V_.P11_slow_[ i ]; S_.y1_slow_[ i ] += V_.PSCInitialValues_slow_[ i ] * spike_value; - // print y1 and y1 slow - // std::cout << "y1: " << S_.y1_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << std::endl; - // print PSCInitlaValues - // std::cout << "PSCInitialValues: " << V_.PSCInitialValues_[ i ] << " PSCInitialValues_slow: " << - // V_.PSCInitialValues_slow_[ i ] << std::endl; } // Update any external currents diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2b6336fae3..06bbc08376 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -387,18 +387,6 @@ class glif_psc_double_alpha : public ArchivingNode return S_.I_syn_; } - // double - // get_I_syn_fast_() const - // { - // return S_.I_syn_fast_; - // } - - // double - // get_I_syn_slow_() const - // { - // return S_.I_syn_slow_; - // } - double get_threshold_() const { From f5010ee5fcde878bf29c38fca8e64ba1bb2f48da Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 15 Aug 2023 17:30:04 -0700 Subject: [PATCH 03/45] Variable name correction and adding tests Some of the parameters that are associated with the fast synaptic components get '_fast' added to their names. Pytest script added to test the GLIF model with double alpha synapses. A couple of test cases require adding extra parameters for the new glif_psc_double_alpha model, so they are added. --- models/glif_psc_double_alpha.cpp | 48 +-- models/glif_psc_double_alpha.h | 39 ++- nestkernel/nest_names.cpp | 1 + nestkernel/nest_names.h | 1 + .../test_neurons_handle_multiplicity.py | 3 + .../test_multimeter_stepping.py | 1 + .../pytests/test_glif_psc_double_alpha.py | 305 ++++++++++++++++++ 7 files changed, 354 insertions(+), 44 deletions(-) create mode 100644 testsuite/pytests/test_glif_psc_double_alpha.py diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index 10011da27b..562154488e 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -80,7 +80,7 @@ nest::glif_psc_double_alpha::Parameters_::Parameters_() , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA , asc_r_( std::vector< double >( 2, 1.0 ) ) - , tau_syn_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless , has_connections_( false ) @@ -107,8 +107,8 @@ nest::glif_psc_double_alpha::State_::State_( const Parameters_& p ) { ASCurrents_sum_ += ASCurrents_[ a ]; } - y1_.clear(); - y2_.clear(); + y1_fast_.clear(); + y2_fast_.clear(); y1_slow_.clear(); y2_slow_.clear(); } @@ -140,8 +140,8 @@ nest::glif_psc_double_alpha::Parameters_::get( DictionaryDatum& d ) const def< std::vector< double > >( d, names::asc_amps, asc_amps_ ); def< std::vector< double > >( d, names::asc_r, asc_r_ ); - ArrayDatum tau_syn_ad( tau_syn_ ); - def< ArrayDatum >( d, names::tau_syn, tau_syn_ad ); + ArrayDatum tau_syn_fast_ad( tau_syn_fast_ ); + def< ArrayDatum >( d, names::tau_syn_fast, tau_syn_fast_ad ); ArrayDatum tau_syn_slow_ad( tau_syn_slow_ ); def< ArrayDatum >( d, names::tau_syn_slow, tau_syn_slow_ad ); @@ -289,7 +289,7 @@ nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* n } const size_t old_n_receptors = this->n_receptors_(); - if ( updateValue< std::vector< double > >( d, names::tau_syn, tau_syn_ ) ) + if ( updateValue< std::vector< double > >( d, names::tau_syn_fast, tau_syn_fast_ ) ) { if ( this->n_receptors_() != old_n_receptors and has_connections_ ) { @@ -297,9 +297,9 @@ nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* n "The neuron has connections, therefore the number of ports cannot be " "reduced." ); } - for ( size_t i = 0; i < tau_syn_.size(); ++i ) + for ( size_t i = 0; i < tau_syn_fast_.size(); ++i ) { - if ( tau_syn_[ i ] <= 0 ) + if ( tau_syn_fast_[ i ] <= 0 ) { throw BadProperty( "All synaptic time constants must be strictly positive." ); } @@ -458,11 +458,11 @@ nest::glif_psc_double_alpha::pre_run_hook() } // postsynaptic currents - V_.P11_.resize( P_.n_receptors_() ); - V_.P21_.resize( P_.n_receptors_() ); - V_.P22_.resize( P_.n_receptors_() ); - V_.P31_.resize( P_.n_receptors_() ); - V_.P32_.resize( P_.n_receptors_() ); + V_.P11_fast_.resize( P_.n_receptors_() ); + V_.P21_fast_.resize( P_.n_receptors_() ); + V_.P22_fast_.resize( P_.n_receptors_() ); + V_.P31_fast_.resize( P_.n_receptors_() ); + V_.P32_fast_.resize( P_.n_receptors_() ); // Additional variables for the slow component V_.P11_slow_.resize( P_.n_receptors_() ); @@ -471,8 +471,8 @@ nest::glif_psc_double_alpha::pre_run_hook() V_.P31_slow_.resize( P_.n_receptors_() ); V_.P32_slow_.resize( P_.n_receptors_() ); - S_.y1_.resize( P_.n_receptors_() ); - S_.y2_.resize( P_.n_receptors_() ); + S_.y1_fast_.resize( P_.n_receptors_() ); + S_.y2_fast_.resize( P_.n_receptors_() ); S_.y1_slow_.resize( P_.n_receptors_() ); // Slow component S_.y2_slow_.resize( P_.n_receptors_() ); // Slow component V_.PSCInitialValues_.resize( P_.n_receptors_() ); @@ -487,13 +487,13 @@ nest::glif_psc_double_alpha::pre_run_hook() for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { // these P are independent - V_.P11_[ i ] = V_.P22_[ i ] = std::exp( -h / P_.tau_syn_[ i ] ); + V_.P11_fast_[ i ] = V_.P22_fast_[ i ] = std::exp( -h / P_.tau_syn_fast_[ i ] ); - V_.P21_[ i ] = h * V_.P11_[ i ]; + V_.P21_fast_[ i ] = h * V_.P11_fast_[ i ]; // these are determined according to a numeric stability criterion // input time parameter shall be in ms, capacity in pF - std::tie( V_.P31_[ i ], V_.P32_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); // Slow component V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); @@ -502,7 +502,7 @@ nest::glif_psc_double_alpha::pre_run_hook() IAFPropagatorAlpha( P_.tau_syn_slow_[ i ], Tau_, P_.C_m_ ).evaluate( h ); - V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_[ i ]; + V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_fast_[ i ]; V_.PSCInitialValues_slow_[ i ] = 1.0 * numerics::e / P_.tau_syn_slow_[ i ] * P_.amp_slow_[ i ]; B_.spikes_[ i ].resize(); } @@ -557,8 +557,8 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.I_syn_fast_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { - S_.U_ += V_.P31_[ i ] * S_.y1_[ i ] + V_.P32_[ i ] * S_.y2_[ i ]; - S_.I_syn_fast_ += S_.y2_[ i ]; + S_.U_ += V_.P31_fast_[ i ] * S_.y1_fast_[ i ] + V_.P32_fast_[ i ] * S_.y2_fast_[ i ]; + S_.I_syn_fast_ += S_.y2_fast_[ i ]; } // slow component S_.I_syn_slow_ = 0.0; @@ -641,9 +641,9 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const // Apply spikes delivered in this step: The spikes arriving at T+1 have an // immediate effect on the state of the neuron // store B_.spikes_[ i ].get_value( lag ); into a variable - S_.y2_[ i ] = V_.P21_[ i ] * S_.y1_[ i ] + V_.P22_[ i ] * S_.y2_[ i ]; - S_.y1_[ i ] *= V_.P11_[ i ]; - S_.y1_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; + S_.y2_fast_[ i ] = V_.P21_fast_[ i ] * S_.y1_fast_[ i ] + V_.P22_fast_[ i ] * S_.y2_fast_[ i ]; + S_.y1_fast_[ i ] *= V_.P11_fast_[ i ]; + S_.y1_fast_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; // slow component S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 06bbc08376..cdc7449529 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -37,7 +37,8 @@ Short description +++++++++++++++++ -Current-based generalized leaky integrate-and-fire (GLIF) models (from the Allen Institute) +Current-based generalized leaky integrate-and-fire (GLIF) models +(from the Allen Institute) Description +++++++++++ @@ -47,12 +48,12 @@ Description Incoming spike events induce a postsynaptic change of current modeled by sum of two alpha functions [2]_. The alpha function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to -be 1 pA at :math:`t = tau_syn`. The peak current of the slow component of the alpha is -given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha has -a single synapse that is accessible through receptor_port 1. An arbitrary number of +be 1 pA at :math:`t = tau_syn_fast`. The peak current of the slow component of the alpha +is given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha +has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different time constants can be configured by setting the desired time -constants as tau_syn array. The resulting synapses are addressed through receptor_port -1, 2, 3, .... +constants as tau_syn_fast array. The resulting synapses are addressed through +receptor_port 1, 2, 3, .... The five GLIF models are: @@ -167,8 +168,8 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn double vector Rise time constants of the synaptic - alpha function in ms +tau_syn_fast double vector Rise time constants of the faster + synaptic alpha function in ms tau_syn_slow double vector Rise time constants of the slower synaptic alpha function in ms amp_slow double vector Relative amplitude of the slower @@ -263,7 +264,7 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_decay_; //!< predefined time scale in 1/ms std::vector< double > asc_amps_; //!< in pA std::vector< double > asc_r_; //!< coefficient - std::vector< double > tau_syn_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms std::vector< double > amp_slow_; //!< synaptic port time constants in ms @@ -279,7 +280,7 @@ class glif_psc_double_alpha : public ArchivingNode //! boolean flag which indicates whether the neuron has voltage dependent threshold component bool has_theta_voltage_; - size_t n_receptors_() const; //!< Returns the size of tau_syn_ + size_t n_receptors_() const; //!< Returns the size of tau_syn_fast_ Parameters_(); @@ -300,8 +301,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > ASCurrents_; //!< after-spike currents in pA double ASCurrents_sum_; //!< in pA int refractory_steps_; //!< Number of refractory steps remaining - std::vector< double > y1_; //!< synapse current evolution state 1 in pA - std::vector< double > y2_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA @@ -337,20 +338,18 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_refractory_decay_rates_; //!< after spike current decay rates during refractory double phi; //!< threshold voltage component coefficient - std::vector< double > P11_; //!< synaptic current evolution parameter - std::vector< double > P21_; //!< synaptic current evolution parameter - std::vector< double > P22_; //!< synaptic current evolution parameter + std::vector< double > P11_fast_; //!< synaptic current evolution parameter + std::vector< double > P21_fast_; //!< synaptic current evolution parameter + std::vector< double > P22_fast_; //!< synaptic current evolution parameter double P30_; //!< membrane current/voltage evolution parameter double P33_; //!< membrane voltage evolution parameter - std::vector< double > P31_; //!< synaptic/membrane current evolution parameter - std::vector< double > P32_; //!< synaptic/membrane current evolution parameter + std::vector< double > P31_fast_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_fast_; //!< synaptic/membrane current evolution parameter // slow component std::vector< double > P11_slow_; //!< synaptic current evolution parameter std::vector< double > P21_slow_; //!< synaptic current evolution parameter std::vector< double > P22_slow_; //!< synaptic current evolution parameter - double P30_slow_; //!< membrane current/voltage evolution parameter - double P33_slow_; //!< membrane voltage evolution parameter std::vector< double > P31_slow_; //!< synaptic/membrane current evolution parameter std::vector< double > P32_slow_; //!< synaptic/membrane current evolution parameter @@ -418,7 +417,7 @@ class glif_psc_double_alpha : public ArchivingNode inline size_t nest::glif_psc_double_alpha::Parameters_::n_receptors_() const { - return tau_syn_.size(); + return tau_syn_fast_.size(); } inline size_t diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 9fe6b9a993..df28532990 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -502,6 +502,7 @@ const Name tau_spike( "tau_spike" ); const Name tau_stc( "tau_stc" ); const Name tau_syn( "tau_syn" ); const Name tau_syn_ex( "tau_syn_ex" ); +const Name tau_syn_fast( "tau_syn_fast" ); const Name tau_syn_in( "tau_syn_in" ); const Name tau_syn_slow( "tau_syn_slow" ); const Name tau_theta( "tau_theta" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index dc34d576a0..d659ab2c1d 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -527,6 +527,7 @@ extern const Name tau_spike; extern const Name tau_stc; extern const Name tau_syn; extern const Name tau_syn_ex; +extern const Name tau_syn_fast; extern const Name tau_syn_in; extern const Name tau_syn_slow; extern const Name tau_theta; diff --git a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py index e5a73b69a3..b511d314ce 100644 --- a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py +++ b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py @@ -74,6 +74,9 @@ "gif_cond_exp_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "glif_cond": {"params": {"tau_syn": [1.0], "E_rev": [-85.0]}, "receptor_type": 1}, "glif_psc": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, + "glif_psc_double_alpha": { + "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, "receptor_type": 1, + }, "aeif_cond_alpha_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "aeif_cond_beta_multisynapse": { "params": {"E_rev": [0.0], "tau_rise": [1.0], "tau_decay": [1.0]}, diff --git a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py index 5f5451e1a6..13d9634913 100644 --- a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py +++ b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py @@ -64,6 +64,7 @@ "gif_psc_exp_multisynapse": {"receptor_type": 1}, "gif_cond_exp_multisynapse": {"receptor_type": 1}, "glif_psc": {"receptor_type": 1}, + "glif_psc_double_alpha": {"receptor_type": 1}, "iaf_cond_alpha_mc": {"receptor_type": 1}, "glif_cond": {"receptor_type": 1}, "ht_neuron": {"receptor_type": 1}, diff --git a/testsuite/pytests/test_glif_psc_double_alpha.py b/testsuite/pytests/test_glif_psc_double_alpha.py new file mode 100644 index 0000000000..3f4253e8dc --- /dev/null +++ b/testsuite/pytests/test_glif_psc_double_alpha.py @@ -0,0 +1,305 @@ +# -*- coding: utf-8 -*- +# +# test_glif_psc_double_alpha.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +import unittest +import nest + +try: + import scipy.stats + + HAVE_SCIPY = True +except ImportError: + HAVE_SCIPY = False + + +@unittest.skipIf(not HAVE_SCIPY, "SciPy package is not available") +class GLIFPSCTestCase(unittest.TestCase): + def setUp(self): + """ + Clean up and initialize NEST before each test. + """ + nest.ResetKernel() + nest.resolution = 0.01 + nest.rng_seed = 123456 + + def simulate_w_stim(self, model_params): + """ + Runs a one second simulation with different types of input stimulation. + Returns the time-steps, voltages, and collected spike times. + """ + + # Create glif model with double alpha synapses to simulate + # For this model, tau_syn_fast = tau_syn_slow = 2.0, and amp_slow = 1.0, + # so that this model behaves the same as glif_psc, when weight is halved + # of the original, for the connections that go into receptor_type == 1. + nrn = nest.Create("glif_psc_double_alpha", params=model_params) + nrn.set({"tau_syn_fast": [2.0], "tau_syn_slow": [2.0], "amp_slow": [1.0]}) + + # Create and connect inputs to glif model (weight is halved to account for double alpha synapse) + espikes = nest.Create( + "spike_generator", params={"spike_times": [10.0, 100.0, 400.0, 700.0], "spike_weights": [10.0] * 4} + ) + ispikes = nest.Create( + "spike_generator", params={"spike_times": [15.0, 99.0, 300.0, 800.0], "spike_weights": [-10.0] * 4} + ) + + cg = nest.Create( + "step_current_generator", + params={ + "amplitude_values": [ + 100.0, + ], + "amplitude_times": [ + 600.0, + ], + "start": 600.0, + "stop": 900.0, + }, + ) + pg = nest.Create("poisson_generator", params={"rate": 1000.0, "start": 200.0, "stop": 500.0}) + pn = nest.Create("parrot_neuron") + + nest.Connect(espikes, nrn, syn_spec={"receptor_type": 1}) + nest.Connect(ispikes, nrn, syn_spec={"receptor_type": 1}) + nest.Connect(cg, nrn, syn_spec={"weight": 3.0}) + nest.Connect(pg, pn) + nest.Connect(pn, nrn, syn_spec={"weight": 17.5, "receptor_type": 1}) # This is also halved + + # For recording spikes and voltage traces + sr = nest.Create("spike_recorder") + nest.Connect(nrn, sr) + + mm = nest.Create("multimeter", params={"record_from": ["V_m"]}) + nest.Connect(mm, nrn) + + nest.Simulate(1000.0) + + times = nest.GetStatus(mm, "events")[0]["times"] + V_m = nest.GetStatus(mm, "events")[0]["V_m"] + spikes = nest.GetStatus(sr, "events")[0]["times"] + + return times, V_m, spikes + + def ks_assert_spikes(self, spikes, reference_spikes): + """ + Runs a two-sided Kolmogorov-Smirnov statistic test on a set of spikes against a set of reference spikes. + """ + p_value_lim = 0.9 + d_lim = 0.2 + d, p_value = scipy.stats.ks_2samp(spikes, reference_spikes) + print(f"d={d}, p_value={p_value}") + self.assertGreater(p_value, p_value_lim) + self.assertLess(d, d_lim) + + def test_lif(self): + """ + Check LIF model + """ + lif_params = { + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_params) + spikes_expected = [ + 394.67, + 424.25, + 483.25, + 612.99, + 628.73, + 644.47, + 660.21, + 675.95, + 691.69, + 706.3, + 722.0, + 737.74, + 753.48, + 769.22, + 784.96, + 800.7, + 816.72, + 832.46, + 848.2, + 863.94, + 879.68, + 895.42, + ] + + self.ks_assert_spikes(spikes, spikes_expected) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r(self): + """ + Check LIF_R model + """ + lif_r_params = { + "spike_dependent_threshold": True, + "after_spike_currents": False, + "adapting_threshold": False, + "V_m": -78.85, + } + times, V_m, spikes = self.simulate_w_stim(lif_r_params) + expected_spikes = [ + 394.67, + 400.67, + 424.51, + 484.48, + 613.4, + 621.31, + 629.67, + 638.47, + 647.7, + 657.36, + 667.42, + 677.87, + 688.67, + 699.8, + 709.84, + 721.58, + 733.57, + 745.76, + 758.11, + 770.6, + 783.21, + 795.92, + 811.36, + 824.05, + 836.81, + 849.64, + 862.54, + 875.49, + 888.48, + ] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_asc(self): + """ + Check LIF_ASC model + """ + lif_asc_params = { + "spike_dependent_threshold": False, + "after_spike_currents": True, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_asc_params) + expected_spikes = [394.67, 484.81, 614.85, 648.28, 685.18, 725.06, 769.76, 821.69, 876.09] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r_asc(self): + """ + Check LIF_R_ASC model + """ + lif_r_asc_params = { + "spike_dependent_threshold": True, + "after_spike_currents": True, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_r_asc_params) + expected_spikes = [394.67, 485.07, 615.16, 649.52, 689.1, 734.18, 786.34, 845.85] + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r_asc_a(self): + """ + Check LIF_R_ASC_A model + """ + lif_r_asc_a_params = { + "spike_dependent_threshold": True, + "after_spike_currents": True, + "adapting_threshold": True, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_r_asc_a_params) + expected_spikes = [395.44, 615.27, 653.77, 701.42, 768.22, 859.82] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_double_alpha_synapse(self): + """ + Check double alpha synapse + """ + neuron = nest.Create("glif_psc_double_alpha") + + # create 3 spikes as inputs + input1 = nest.Create("spike_generator", 1, {"spike_times": [10.0]}) + input2 = nest.Create("spike_generator", 1, {"spike_times": [210.0]}) + input3 = nest.Create("spike_generator", 1, {"spike_times": [410.0]}) + + # input1: double alpha at the same location with same amplitude. should double the amplitude + # input2: double alpha at the same location with different amplitude. net effect shoulld 1.5 of single + # input3: double alpha at different locations. The peak of the combined input shifts to ther right. + neuron_settings = { + "tau_syn_fast": [2.0, 2.0, 2.0], + "tau_syn_slow": [2.0, 2.0, 4.0], + "amp_slow": [1.0, 0.5, 0.5], + } + + neuron.set(neuron_settings) + # record the current with 0.1 ms resolution to capture the peak values. + multimeter = nest.Create("multimeter", params={"record_from": ["I_syn"], "interval": 0.1}) + + # set the delay to 1 ms. and amplitude is 1. + # the expected peak time of the fast components are, 13 ms, 213 ms, 413 ms. + nest.Connect(input1, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 1}) + nest.Connect(input2, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 2}) + nest.Connect(input3, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 3}) + nest.Connect(multimeter, neuron) + + # Do simulation. + nest.Simulate(500.0) + times = nest.GetStatus(multimeter, "events")[0]["times"] + I_syn = nest.GetStatus(multimeter, "events")[0]["I_syn"] + + # Check the results. + # the peak timing should be 13 ms (index 129), 213 ms (index 2129), + # and slightly above 413 ms (413.4 ms, or index 4133). + + # for the first two we can estimate the exact peak values (2 and 1.5 pA) + self.assertAlmostEqual(I_syn[129], 2.0, delta=0.0001) + self.assertAlmostEqual(I_syn[2129], 1.5, delta=0.0001) + # for the last one, check if the the new peak is larger than peak of the fast component + self.assertGreater(I_syn[4133], I_syn[4129]) # peak is shifted to the right. + + +def suite(): + return unittest.makeSuite(GLIFPSCTestCase, "test") + + +def run(): + runner = unittest.TextTestRunner(verbosity=2) + runner.run(suite()) + + +if __name__ == "__main__": + run() From 4d1c60e0b89682d0afd141c9f18757bec7e46f49 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Fri, 25 Aug 2023 18:35:13 -0700 Subject: [PATCH 04/45] Making double-alpha description clearer --- models/glif_psc.h | 2 +- models/glif_psc_double_alpha.h | 26 +++++++++++++++----------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/models/glif_psc.h b/models/glif_psc.h index cb11214933..e85c48f182 100644 --- a/models/glif_psc.h +++ b/models/glif_psc.h @@ -165,7 +165,7 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn double vector Rise time constants of the synaptic +tau_syn double vector Time constants of the synaptic alpha function in ms E_rev double vector Reversal potential in mV spike_dependent_threshold bool flag whether the neuron has diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index cdc7449529..232929a079 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -44,16 +44,20 @@ Description +++++++++++ ``glif_psc_double_alpha`` provides five generalized leaky integrate-and-fire -(GLIF) models [1]_ with alpha-function shaped synaptic currents. +(GLIF) models [1]_ with double alpha-function shaped synaptic currents. Incoming spike events induce a postsynaptic change of current modeled -by sum of two alpha functions [2]_. The alpha function is normalized such that an event -of weight 1.0 results in a peak current of the fast component of the alpha function to -be 1 pA at :math:`t = tau_syn_fast`. The peak current of the slow component of the alpha -is given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha -has a single synapse that is accessible through receptor_port 1. An arbitrary number of -synapses with different time constants can be configured by setting the desired time -constants as tau_syn_fast array. The resulting synapses are addressed through -receptor_port 1, 2, 3, .... +by the sum of two alpha functions (fast and slow components) for each receptor [2]_. +This function is normalized such that an event of weight 1.0 results in a peak current +of the fast component of the alpha function to be 1 pA at :math:`t = tau_syn_fast`. +The relative peak current of the slow component is given as amp_slow, at +:math:`t = tau_syn_slow`. Namely, +:math:`I_{syn} = alpha_function(tau_syn=tau_syn_fast) + amp_slow * +alpha_function(tau_syn=tau_syn_slow)`. Therefore if amp_slow is not 0, the peak current +of the total synaptic current is larger than the specified weight. By default, +glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. +An arbitrary number of synapses with different time constants and amp_slow can be +configured by setting the desired parameters of tau_syn_fast, tau_syn_slow, and amp_slow +arrays. The resulting synapses are addressed through receptor_port 1, 2, 3, .... The five GLIF models are: @@ -168,9 +172,9 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn_fast double vector Rise time constants of the faster +tau_syn_fast double vector Time constants of the faster synaptic alpha function in ms -tau_syn_slow double vector Rise time constants of the slower +tau_syn_slow double vector Time constants of the slower synaptic alpha function in ms amp_slow double vector Relative amplitude of the slower synaptic alpha function From 4bb6edbf2cfbc13f5dc29c0405c1455a118b6811 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 29 Aug 2023 18:45:44 -0700 Subject: [PATCH 05/45] Making math format LaTeX compatible --- models/glif_psc_double_alpha.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 232929a079..f818ac80d3 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -48,16 +48,18 @@ Description Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current -of the fast component of the alpha function to be 1 pA at :math:`t = tau_syn_fast`. +of the fast component of the alpha function to be 1 pA at +:math:`t = \tau_\text{syn\_fast}`. The relative peak current of the slow component is given as amp_slow, at -:math:`t = tau_syn_slow`. Namely, -:math:`I_{syn} = alpha_function(tau_syn=tau_syn_fast) + amp_slow * -alpha_function(tau_syn=tau_syn_slow)`. Therefore if amp_slow is not 0, the peak current -of the total synaptic current is larger than the specified weight. By default, -glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. -An arbitrary number of synapses with different time constants and amp_slow can be -configured by setting the desired parameters of tau_syn_fast, tau_syn_slow, and amp_slow -arrays. The resulting synapses are addressed through receptor_port 1, 2, 3, .... +:math:`t = \tau_\text{syn\_slow}`. Namely, +:math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + +\text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. +Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger +than the specified weight. By default, glif_psc_double_alpha has a single synapse that +is accessible through receptor_port 1. An arbitrary number of synapses with different +time constants and amp_slow can be configured by setting the desired parameters of +tau_syn_fast, tau_syn_slow, and amp_slow arrays. The resulting synapses are addressed +through receptor_port 1, 2, 3, .... The five GLIF models are: From 5500b1591928c3a209dc0812b9100a81869a1090 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:30 -0700 Subject: [PATCH 06/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index f818ac80d3..202e650078 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -49,7 +49,7 @@ Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at -:math:`t = \tau_\text{syn\_fast}`. +:math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn\_slow}`. Namely, :math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + From 7222dd8a5852c4c2ba31b3649f54161d2f48fae9 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:45 -0700 Subject: [PATCH 07/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 202e650078..84c15b322f 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -51,7 +51,7 @@ This function is normalized such that an event of weight 1.0 results in a peak c of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at -:math:`t = \tau_\text{syn\_slow}`. Namely, +:math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + \text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger From e2a3aacb35c9e234da90fc21f8e58144191a4b43 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:55 -0700 Subject: [PATCH 08/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 84c15b322f..93cb4e5766 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -52,7 +52,7 @@ of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, -:math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + +:math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + \text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that From c71e76c2cbc1cf0ae49b7da548d63b9f4cd06f9b Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:29:06 -0700 Subject: [PATCH 09/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 93cb4e5766..749105cfe9 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -53,7 +53,7 @@ of the fast component of the alpha function to be 1 pA at The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + -\text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. +\text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different From 851f072f3219c662bed40b0ce68c2092ce43f52a Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Wed, 30 Aug 2023 13:40:51 -0700 Subject: [PATCH 10/45] black applied to a modified file --- .../pytests/sli2py_neurons/test_neurons_handle_multiplicity.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py index b511d314ce..e595490aff 100644 --- a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py +++ b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py @@ -75,7 +75,8 @@ "glif_cond": {"params": {"tau_syn": [1.0], "E_rev": [-85.0]}, "receptor_type": 1}, "glif_psc": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "glif_psc_double_alpha": { - "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, "receptor_type": 1, + "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, + "receptor_type": 1, }, "aeif_cond_alpha_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "aeif_cond_beta_multisynapse": { From 222a6176e4b10067cc6a7940d7218a7d21c0ed80 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Wed, 30 Aug 2023 13:45:30 -0700 Subject: [PATCH 11/45] clang-format -i applied to committed files --- models/glif_psc_double_alpha.cpp | 5 +++-- models/glif_psc_double_alpha.h | 10 +++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index 562154488e..6245501a88 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -80,7 +80,7 @@ nest::glif_psc_double_alpha::Parameters_::Parameters_() , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA , asc_r_( std::vector< double >( 2, 1.0 ) ) - , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless , has_connections_( false ) @@ -493,7 +493,8 @@ nest::glif_psc_double_alpha::pre_run_hook() // these are determined according to a numeric stability criterion // input time parameter shall be in ms, capacity in pF - std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = + IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); // Slow component V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 749105cfe9..2be27965e7 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -270,7 +270,7 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_decay_; //!< predefined time scale in 1/ms std::vector< double > asc_amps_; //!< in pA std::vector< double > asc_r_; //!< coefficient - std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms std::vector< double > amp_slow_; //!< synaptic port time constants in ms @@ -307,8 +307,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > ASCurrents_; //!< after-spike currents in pA double ASCurrents_sum_; //!< in pA int refractory_steps_; //!< Number of refractory steps remaining - std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA - std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA @@ -347,8 +347,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > P11_fast_; //!< synaptic current evolution parameter std::vector< double > P21_fast_; //!< synaptic current evolution parameter std::vector< double > P22_fast_; //!< synaptic current evolution parameter - double P30_; //!< membrane current/voltage evolution parameter - double P33_; //!< membrane voltage evolution parameter + double P30_; //!< membrane current/voltage evolution parameter + double P33_; //!< membrane voltage evolution parameter std::vector< double > P31_fast_; //!< synaptic/membrane current evolution parameter std::vector< double > P32_fast_; //!< synaptic/membrane current evolution parameter From b2940b5bd56479e4eee3a4c72c5f525f3d3a8153 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:28:54 -0700 Subject: [PATCH 12/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2be27965e7..a252d7e3c8 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -37,7 +37,7 @@ Short description +++++++++++++++++ -Current-based generalized leaky integrate-and-fire (GLIF) models +Current-based generalized leaky integrate-and-fire (GLIF) models with double alpha-function (from the Allen Institute) Description From ecd4522d0188df34ed833b867045ea7d367a8a49 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:30:13 -0700 Subject: [PATCH 13/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index a252d7e3c8..2298289b57 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -49,7 +49,7 @@ Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at -:math:`t = \tau_\text{syn_fast}`. +:math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + From 097ecf62fdc4d61f038901e591846c084295bef2 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:30:29 -0700 Subject: [PATCH 14/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2298289b57..4ca0b54dee 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -50,7 +50,7 @@ by the sum of two alpha functions (fast and slow components) for each receptor [ This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. -The relative peak current of the slow component is given as amp_slow, at +The relative peak current of the slow component is given as ``amp_slow``, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + \text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. From 58b489d6fcc52a2ebb753d02818ac3389b2445f7 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:00 -0700 Subject: [PATCH 15/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 4ca0b54dee..80c21a2cbc 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -52,8 +52,11 @@ of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as ``amp_slow``, at :math:`t = \tau_\text{syn_slow}`. Namely, -:math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + -\text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. + +.. math:: + + I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). + Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different From 5084585d8fc5b778eae6fabed67b474b10a5f05e Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:15 -0700 Subject: [PATCH 16/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 80c21a2cbc..920a2bd326 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -57,12 +57,12 @@ The relative peak current of the slow component is given as ``amp_slow``, at I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). -Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger -than the specified weight. By default, glif_psc_double_alpha has a single synapse that -is accessible through receptor_port 1. An arbitrary number of synapses with different -time constants and amp_slow can be configured by setting the desired parameters of -tau_syn_fast, tau_syn_slow, and amp_slow arrays. The resulting synapses are addressed -through receptor_port 1, 2, 3, .... +Therefore if ``amp_slow`` is not 0, the peak current of the total synaptic current is larger +than the specified weight. By default, ``glif_psc_double_alpha`` has a single synapse that +is accessible through ``receptor_port`` 1. An arbitrary number of synapses with different +time constants and ``amp_slow`` can be configured by setting the desired parameters of +``tau_syn_fast``, ``tau_syn_slow``, and ``amp_slow`` arrays. The resulting synapses are addressed +through ``receptor_port`` 1, 2, 3, .... The five GLIF models are: From 0d752cdb35e8db0a3d90235bec7e3cd0593ad96b Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:30 -0700 Subject: [PATCH 17/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 920a2bd326..b4bcbde9f3 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -105,7 +105,7 @@ F, s, A) to NEST used units (i.e., mV, nS (1/GOhm), pF, ms, pA) and values being rounded to appropriate digits for simplification. For models with spike dependent threshold (i.e., GLIF2, GLIF4 and GLIF5), -parameter setting of voltage_reset_fraction and voltage_reset_add may lead to the +parameter setting of ``voltage_reset_fraction`` and ``voltage_reset_add`` may lead to the situation that voltage is bigger than threshold after reset. In this case, the neuron will continue to spike until the end of the simulation regardless the stimulated inputs. We recommend the setting of the parameters of these three models to follow the From 56397f3e9646674b2d848f551a580a305f42d74c Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Thu, 31 Aug 2023 10:34:49 -0700 Subject: [PATCH 18/45] Adding double_alpha reference to glif_psc.h --- models/glif_psc.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/models/glif_psc.h b/models/glif_psc.h index e85c48f182..d133c9f6b4 100644 --- a/models/glif_psc.h +++ b/models/glif_psc.h @@ -192,7 +192,8 @@ References See also ++++++++ -gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp +gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp, +glif_psc_double_alpha EndUserDocs */ From 35dbf5e3d1d3c59fa63d98ba75c54d05553f5982 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Thu, 31 Aug 2023 10:38:32 -0700 Subject: [PATCH 19/45] One more change for math notation --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index b4bcbde9f3..d0ac80de1f 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -51,7 +51,7 @@ This function is normalized such that an event of weight 1.0 results in a peak c of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as ``amp_slow``, at -:math:`t = \tau_\text{syn_slow}`. Namely, +:math:`t = \tau_\text{syn, slow}`. Namely, .. math:: From 2dc50d82ebcc4e6796ea78a54754bbe587cef026 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Fri, 1 Sep 2023 11:27:11 -0700 Subject: [PATCH 20/45] Adding example, and reapplying clang-format --- doc/htmldoc/examples/index.rst | 2 + models/glif_psc_double_alpha.h | 3 +- .../examples/glif_psc_double_alpha_neuron.py | 214 ++++++++++++++++++ 3 files changed, 218 insertions(+), 1 deletion(-) create mode 100644 pynest/examples/glif_psc_double_alpha_neuron.py diff --git a/doc/htmldoc/examples/index.rst b/doc/htmldoc/examples/index.rst index 549e6024cc..d37f8b4102 100644 --- a/doc/htmldoc/examples/index.rst +++ b/doc/htmldoc/examples/index.rst @@ -76,6 +76,7 @@ PyNEST examples * :doc:`../auto_examples/glif_cond_neuron` * :doc:`../auto_examples/glif_psc_neuron` + * :doc:`../auto_examples/glif_psc_double_alpha_neuron` @@ -260,6 +261,7 @@ PyNEST examples ../auto_examples/mc_neuron ../auto_examples/glif_cond_neuron ../auto_examples/glif_psc_neuron + ../auto_examples/glif_psc_double_alpha_neuron ../auto_examples/precise_spiking ../auto_examples/CampbellSiegert ../auto_examples/vinit_example diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index d0ac80de1f..b0589d03b2 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -55,7 +55,8 @@ The relative peak current of the slow component is given as ``amp_slow``, at .. math:: - I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). + I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot +\text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). Therefore if ``amp_slow`` is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, ``glif_psc_double_alpha`` has a single synapse that diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py new file mode 100644 index 0000000000..e448382bea --- /dev/null +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -0,0 +1,214 @@ +# -*- coding: utf-8 -*- +# +# glif_psc_double_alpha_neuron.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +# %% +""" +Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha +synaptic function example +------------------------------------------------------------------------ + +Simple example of how to use the ``glif_psc_double_alpha`` neuron model that illustrates +difference from the ``glif_psc`` neuron model. + +The behavior of the ``glif_psc_double_alpha`` neuron model the same as the ``glif_psc`` +neuron model, except that the synaptic currents are modeled as a double alpha function. +Therefore, in this example, we only compare the difference in the synaptic currents +between the two models. compared to the single alpha function, the double alpha function +has much more control over the shape of the tail of the synaptic current. + +Simple synaptic inputs are applied to the neuron and the resulting voltage traces, +current traces are shown for the two models. + +""" + +############################################################################## +# First, we import all necessary modules to simulate, analyze and plot this +# example. + +import nest +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec + +############################################################################## +# We initialize NEST and set the simulation resolution. + +nest.ResetKernel() +resolution = 0.05 +nest.resolution = resolution + +############################################################################## +# We also pre-define the synapse time constant arrays. +# In contrast to ``glif_psc`` models, ``glif_psc_double_alpha`` models have +# two components of synaptic currents, one for the fast component and the other +# for the slow component. The relative amplitude also need to be set, so there +# are three parameters to define per port. For this example, we keep the +# ``tau_syn_fast`` to 2 ms for simplicity, and vary the ``tau_syn_slow`` and +# ``amp_slow`` to illustrate how the parameters affect the shape of the synaptic +# currents. + +tau_syn_glif_psc = [2.0, 2.0, 2.0] # value for the ``glif_psc`` model + +tau_syn_fast = [2.0, 2.0, 2.0] # common between 'timing' and 'amp' manipulations + +# for the slow component timing manipuation +tau_syn_slow_timing = [4.0, 6.0, 8.0] +amp_slow_timing = [0.5, 0.5, 0.5] + +# for the slow component amplitude manipulation +tau_syn_slow_amp = [6.0, 6.0, 6.0] +amp_slow_amp = [0.2, 0.5, 0.8] + +############################################################################### +# Now we create three neurons: ``glif_psc``, ``glif_psc_double_alpha_timing``, +# and ``glif_psc_double_alpha_amp``. The parameters for the ``glif_psc`` neuron +# are set as default. The parameters for the ``glif_psc_double_alpha_timing`` +# neuron are set to have the same ``tau_syn_fast`` as the ``glif_psc`` neuron, +# and the ``tau_syn_slow`` and ``amp_slow`` are set to the values defined above +# for the timing manipulation. + +n_glif_psc = nest.Create( + "glif_psc", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn": tau_syn_glif_psc, + }, +) + +n_glif_psc_double_alpha_timing = nest.Create( + "glif_psc_double_alpha", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn_fast": tau_syn_fast, + "tau_syn_slow": tau_syn_slow_timing, + "amp_slow": amp_slow_timing, + }, +) + +n_glif_psc_double_alpha_amp = nest.Create( + "glif_psc_double_alpha", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn_fast": tau_syn_fast, + "tau_syn_slow": tau_syn_slow_amp, + "amp_slow": amp_slow_amp, + }, +) + +neurons = n_glif_psc + n_glif_psc_double_alpha_timing + n_glif_psc_double_alpha_amp + +############################################################################### +# For the stimulation input to the glif_psc neurons, we create three excitation +# spike generators, each one with a single spike. + +espike1 = nest.Create("spike_generator", params={"spike_times": [10.0], "spike_weights": [20.0]}) +espike2 = nest.Create("spike_generator", params={"spike_times": [110.0], "spike_weights": [20.0]}) +espike3 = nest.Create("spike_generator", params={"spike_times": [210.0], "spike_weights": [20.0]}) + +############################################################################### +# The generators are then connected to the neurons. Specification of +# the ``receptor_type`` uniquely defines the target receptor. +# We connect each of the spikes generator to a different receptor that have different +# parameters. + +nest.Connect(espike1, neurons, syn_spec={"delay": resolution, "receptor_type": 1}) +nest.Connect(espike2, neurons, syn_spec={"delay": resolution, "receptor_type": 2}) +nest.Connect(espike3, neurons, syn_spec={"delay": resolution, "receptor_type": 3}) + +############################################################################### +# A ``multimeter`` is created and connected to the neurons. The parameters +# specified for the multimeter include the list of quantities that should be +# recorded and the time interval at which quantities are measured. + +mm = nest.Create( + "multimeter", + params={ + "interval": resolution, + "record_from": ["V_m", "I_syn"], + }, +) +nest.Connect(mm, neurons) + +############################################################################### +# Run the simulation for 300 ms and retrieve recorded data from +# the multimeter and spike recorder. + +nest.Simulate(300.0) +data = mm.events + +############################################################################### +# We plot the time traces of the synaptic current and the membrane potential. +# Each input current is annotated with the corresponding parameter value of the +# receptor. The blue line is the synaptic current of the ``glif_psc`` neuron, and +# the red line is the synaptic current of the ``glif_psc_double_alpha`` neuron. + + +# defining the figure property for each parameter variation type, +variation_types = ["timing", "amp"] +annotate_variable = ["tau_syn_slow", "amp_slow"] +annotate_values = [tau_syn_slow_timing, amp_slow_amp] +fig_titles = [ + "Variation of tau_syn_slow: tau_syn_fast = 2.0, amp_slow = 0.5", + "Variation of amp_slow: tau_syn_fast = 2.0, tau_syn_slow = 6.0", +] + + +senders = data["senders"] +t = data["times"][senders == 1] + +for i, variation_type in enumerate(variation_types): + plt.figure(variation_type, figsize=(10, 5)) + gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1]) + data_types = ["I_syn", "V_m"] + data_types_names = ["Synaptic current (pA)", "Membrane potential (mV)"] + for j, data_type in enumerate(data_types): + d = data[data_type] + ax = plt.subplot(gs[j]) + ax.plot(t, d[senders == 1], "b", label="glif_psc (tau_syn=2.0)") + ax.plot(t, d[senders == 2 + i], "r", label="glif_psc_double_alpha") + if j == 0: + # place legend outside the plot + ax.legend(bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0) + else: + ax.set_xlabel("time (ms)") + + ax.set_ylabel(data_types_names[j]) + + # now let's annotate each of the input with the corresponding parameter. + spike_timings = [10.0, 110.0, 210.0] + ax = plt.subplot(gs[0]) + for j, spike_timing in enumerate(spike_timings): + ax.annotate( + f"{annotate_variable[i]}={annotate_values[i][j]}", + xy=(spike_timing + 10, 20), + xytext=(spike_timing + 10, 25), + arrowprops=dict(arrowstyle="->", connectionstyle="arc3"), + ) + plt.title(fig_titles[i]) + plt.tight_layout() + + +plt.show() From f24d0de0ff684c9e836877c02fc4ac62214d8969 Mon Sep 17 00:00:00 2001 From: Dennis Terhorst Date: Mon, 4 Sep 2023 11:26:07 +0200 Subject: [PATCH 21/45] Apply suggestions from code review Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- pynest/examples/glif_psc_double_alpha_neuron.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py index e448382bea..50895ccb70 100644 --- a/pynest/examples/glif_psc_double_alpha_neuron.py +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -19,7 +19,6 @@ # You should have received a copy of the GNU General Public License # along with NEST. If not, see . -# %% """ Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha synaptic function example @@ -28,15 +27,14 @@ Simple example of how to use the ``glif_psc_double_alpha`` neuron model that illustrates difference from the ``glif_psc`` neuron model. -The behavior of the ``glif_psc_double_alpha`` neuron model the same as the ``glif_psc`` +The behavior of the ``glif_psc_double_alpha`` neuron model is the same as the ``glif_psc`` neuron model, except that the synaptic currents are modeled as a double alpha function. Therefore, in this example, we only compare the difference in the synaptic currents -between the two models. compared to the single alpha function, the double alpha function +between the two models. Compared to the single alpha function, the double alpha function has much more control over the shape of the tail of the synaptic current. -Simple synaptic inputs are applied to the neuron and the resulting voltage traces, +Simple synaptic inputs are applied to the neuron and the resulting voltage and current traces are shown for the two models. - """ ############################################################################## @@ -58,8 +56,8 @@ # We also pre-define the synapse time constant arrays. # In contrast to ``glif_psc`` models, ``glif_psc_double_alpha`` models have # two components of synaptic currents, one for the fast component and the other -# for the slow component. The relative amplitude also need to be set, so there -# are three parameters to define per port. For this example, we keep the +# for the slow component. The relative amplitude also needs to be set, so there +# are three parameters to define per receptor port. For this example, we keep the # ``tau_syn_fast`` to 2 ms for simplicity, and vary the ``tau_syn_slow`` and # ``amp_slow`` to illustrate how the parameters affect the shape of the synaptic # currents. @@ -121,7 +119,7 @@ neurons = n_glif_psc + n_glif_psc_double_alpha_timing + n_glif_psc_double_alpha_amp ############################################################################### -# For the stimulation input to the glif_psc neurons, we create three excitation +# For the stimulation input to the ``glif_psc`` neurons, we create three excitation # spike generators, each one with a single spike. espike1 = nest.Create("spike_generator", params={"spike_times": [10.0], "spike_weights": [20.0]}) From bf917690f219aad59cc86b8992a2b2a95b36d04f Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 18 Jul 2023 15:51:06 -0700 Subject: [PATCH 22/45] Adding GLIF with double alpha psc. --- models/glif_psc_double_alpha.cpp | 711 +++++++++++++++++++++++++++++++ models/glif_psc_double_alpha.h | 494 +++++++++++++++++++++ modelsets/full | 1 + nestkernel/nest_names.cpp | 2 + nestkernel/nest_names.h | 2 + 5 files changed, 1210 insertions(+) create mode 100644 models/glif_psc_double_alpha.cpp create mode 100644 models/glif_psc_double_alpha.h diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp new file mode 100644 index 0000000000..b48b75d7b4 --- /dev/null +++ b/models/glif_psc_double_alpha.cpp @@ -0,0 +1,711 @@ +/* + * glif_psc_double_alpha.cpp + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#include "glif_psc_double_alpha.h" + +// C++ includes: +#include +#include + +// Includes from libnestutil: +#include "dict_util.h" +#include "exceptions.h" +#include "iaf_propagator.h" +#include "kernel_manager.h" +#include "universal_data_logger_impl.h" + +// Includes from sli: +#include "dict.h" +#include "dictutils.h" + +using namespace nest; + +nest::RecordablesMap< nest::glif_psc_double_alpha > nest::glif_psc_double_alpha::recordablesMap_; + +namespace nest +{ +// Override the create() method with one call to RecordablesMap::insert_() +// for each quantity to be recorded. +template <> +void +RecordablesMap< nest::glif_psc_double_alpha >::create() +{ + insert_( names::V_m, &nest::glif_psc_double_alpha::get_V_m_ ); + insert_( names::ASCurrents_sum, &nest::glif_psc_double_alpha::get_ASCurrents_sum_ ); + insert_( names::I, &nest::glif_psc_double_alpha::get_I_ ); + insert_( names::I_syn, &nest::glif_psc_double_alpha::get_I_syn_ ); + insert_( names::threshold, &nest::glif_psc_double_alpha::get_threshold_ ); + insert_( names::threshold_spike, &nest::glif_psc_double_alpha::get_threshold_spike_ ); + insert_( names::threshold_voltage, &nest::glif_psc_double_alpha::get_threshold_voltage_ ); +} +} + +/* ---------------------------------------------------------------- + * Default constructors defining default parameters and state + * ---------------------------------------------------------------- */ + +nest::glif_psc_double_alpha::Parameters_::Parameters_() + : G_( 9.43 ) // in nS + , E_L_( -78.85 ) // in mV + , th_inf_( -51.68 - E_L_ ) // in mv, rel to E_L_, - 51.68 - E_L_, i.e., 27.17 + , C_m_( 58.72 ) // in pF + , t_ref_( 3.75 ) // in ms + , V_reset_( 0.0 ) // in mV, rel to E_L_, -78.85 - E_L_ + , th_spike_add_( 0.37 ) // in mV + , th_spike_decay_( 0.009 ) // in 1/ms + , voltage_reset_fraction_( 0.20 ) + , voltage_reset_add_( 18.51 ) // in mV + , th_voltage_index_( 0.005 ) // in 1/ms + , th_voltage_decay_( 0.09 ) // in 1/ms + , asc_init_( std::vector< double >( 2, 0.0 ) ) // in pA + , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms + , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA + , asc_r_( std::vector< double >( 2, 1.0 ) ) + , tau_syn_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms + , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless + , has_connections_( false ) + , has_theta_spike_( false ) + , has_asc_( false ) + , has_theta_voltage_( false ) +{ +} + +nest::glif_psc_double_alpha::State_::State_( const Parameters_& p ) + : U_( 0.0 ) // in mV + , threshold_( p.th_inf_ ) // in mV + , threshold_spike_( 0.0 ) // in mV + , threshold_voltage_( 0.0 ) // in mV + , I_( 0.0 ) // in pA + , I_syn_( 0.0 ) // in pA + , I_syn_fast_( 0.0 ) // in pA + , I_syn_slow_( 0.0 ) // in pA + , ASCurrents_( p.asc_init_ ) // in pA + , ASCurrents_sum_( 0.0 ) // in pA + , refractory_steps_( 0 ) +{ + for ( std::size_t a = 0; a < p.asc_init_.size(); ++a ) + { + ASCurrents_sum_ += ASCurrents_[ a ]; + } + y1_.clear(); + y2_.clear(); + y1_slow_.clear(); + y2_slow_.clear(); +} + +/* ---------------------------------------------------------------- + * Parameter and state extractions and manipulation functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::Parameters_::get( DictionaryDatum& d ) const +{ + def< double >( d, names::V_th, th_inf_ + E_L_ ); + def< double >( d, names::g, G_ ); + def< double >( d, names::E_L, E_L_ ); + def< double >( d, names::C_m, C_m_ ); + def< double >( d, names::t_ref, t_ref_ ); + def< double >( d, names::V_reset, V_reset_ + E_L_ ); + + def< double >( d, names::th_spike_add, th_spike_add_ ); + def< double >( d, names::th_spike_decay, th_spike_decay_ ); + def< double >( d, names::voltage_reset_fraction, voltage_reset_fraction_ ); + def< double >( d, names::voltage_reset_add, voltage_reset_add_ ); + + def< double >( d, names::th_voltage_index, th_voltage_index_ ); + def< double >( d, names::th_voltage_decay, th_voltage_decay_ ); + + def< std::vector< double > >( d, names::asc_init, asc_init_ ); + def< std::vector< double > >( d, names::asc_decay, asc_decay_ ); + def< std::vector< double > >( d, names::asc_amps, asc_amps_ ); + def< std::vector< double > >( d, names::asc_r, asc_r_ ); + + ArrayDatum tau_syn_ad( tau_syn_ ); + def< ArrayDatum >( d, names::tau_syn, tau_syn_ad ); + + ArrayDatum tau_syn_slow_ad( tau_syn_slow_ ); + def< ArrayDatum >( d, names::tau_syn_slow, tau_syn_slow_ad ); + + ArrayDatum amp_slow_ad( amp_slow_ ); + def< ArrayDatum >( d, names::amp_slow, amp_slow_ad ); + + def< bool >( d, names::has_connections, has_connections_ ); + def< bool >( d, names::spike_dependent_threshold, has_theta_spike_ ); + def< bool >( d, names::after_spike_currents, has_asc_ ); + def< bool >( d, names::adapting_threshold, has_theta_voltage_ ); +} + +double +nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* node ) +{ + // if E_L_ is changed, we need to adjust all variables defined relative to + // E_L_ + const double ELold = E_L_; + updateValueParam< double >( d, names::E_L, E_L_, node ); + const double delta_EL = E_L_ - ELold; + + if ( updateValueParam< double >( d, names::V_reset, V_reset_, node ) ) + { + V_reset_ -= E_L_; + } + else + { + V_reset_ -= delta_EL; + } + + if ( updateValueParam< double >( d, names::V_th, th_inf_, node ) ) + { + th_inf_ -= E_L_; + } + else + { + th_inf_ -= delta_EL; + } + + updateValueParam< double >( d, names::g, G_, node ); + updateValueParam< double >( d, names::C_m, C_m_, node ); + updateValueParam< double >( d, names::t_ref, t_ref_, node ); + + updateValueParam< double >( d, names::th_spike_add, th_spike_add_, node ); + updateValueParam< double >( d, names::th_spike_decay, th_spike_decay_, node ); + updateValueParam< double >( d, names::voltage_reset_fraction, voltage_reset_fraction_, node ); + updateValueParam< double >( d, names::voltage_reset_add, voltage_reset_add_, node ); + + updateValueParam< double >( d, names::th_voltage_index, th_voltage_index_, node ); + updateValueParam< double >( d, names::th_voltage_decay, th_voltage_decay_, node ); + + updateValue< std::vector< double > >( d, names::asc_init, asc_init_ ); + updateValue< std::vector< double > >( d, names::asc_decay, asc_decay_ ); + updateValue< std::vector< double > >( d, names::asc_amps, asc_amps_ ); + updateValue< std::vector< double > >( d, names::asc_r, asc_r_ ); + + // set model mechanisms + updateValueParam< bool >( d, names::spike_dependent_threshold, has_theta_spike_, node ); + updateValueParam< bool >( d, names::after_spike_currents, has_asc_, node ); + updateValueParam< bool >( d, names::adapting_threshold, has_theta_voltage_, node ); + + // check model mechanisms parameter + if ( not( ( not has_theta_spike_ and not has_asc_ and not has_theta_voltage_ ) or // glif1 + ( has_theta_spike_ and not has_asc_ and not has_theta_voltage_ ) or // glif2 + ( not has_theta_spike_ and has_asc_ and not has_theta_voltage_ ) or // glif3 + ( has_theta_spike_ and has_asc_ and not has_theta_voltage_ ) or // glif4 + ( has_theta_spike_ and has_asc_ and has_theta_voltage_ ) // glif5 + ) ) + { + throw BadProperty( + "Incorrect model mechanism combination setting." + "See documentation for setting of model mechanism parameters:" + "spike_dependent_threshold, after_spike_currents, adapting_threshold." ); + } + + // check after ASC parameters' sizes and values + if ( has_asc_ ) + { + // check size + const size_t asc_size = asc_decay_.size(); + if ( not( + ( asc_init_.size() == asc_size ) and ( asc_amps_.size() == asc_size ) and ( asc_r_.size() == asc_size ) ) ) + { + throw BadProperty( + "All after spike current parameters (i.e., asc_init, k, asc_amps, r) " + "must have the same size." ); + } + + // check values + for ( std::size_t a = 0; a < asc_decay_.size(); ++a ) + { + if ( asc_decay_[ a ] <= 0.0 ) + { + throw BadProperty( "After-spike current time constant must be strictly positive." ); + } + + if ( asc_r_[ a ] < 0.0 or asc_r_[ a ] > 1.0 ) + { + throw BadProperty( "After spike current fraction following spike coefficients r must be within [0.0, 1.0]." ); + } + } + } + + if ( V_reset_ >= th_inf_ ) + { + throw BadProperty( "Reset potential must be smaller than threshold." ); + } + + if ( C_m_ <= 0.0 ) + { + throw BadProperty( "Capacitance must be strictly positive." ); + } + + if ( G_ <= 0.0 ) + { + throw BadProperty( "Membrane conductance must be strictly positive." ); + } + + if ( t_ref_ <= 0.0 ) + { + throw BadProperty( "Refractory time constant must be strictly positive." ); + } + + if ( has_theta_voltage_ ) + { + if ( th_voltage_decay_ <= 0.0 ) + { + throw BadProperty( "Voltage-induced threshold time constant must be strictly positive." ); + } + } + + // check spike component parameters + if ( has_theta_spike_ ) + { + if ( th_spike_decay_ <= 0.0 ) + { + throw BadProperty( "Spike induced threshold time constant must be strictly positive." ); + } + + if ( voltage_reset_fraction_ < 0.0 or voltage_reset_fraction_ > 1.0 ) + { + throw BadProperty( "Voltage fraction coefficient following spike must be within [0.0, 1.0]." ); + } + } + + const size_t old_n_receptors = this->n_receptors_(); + if ( updateValue< std::vector< double > >( d, names::tau_syn, tau_syn_ ) ) + { + if ( this->n_receptors_() != old_n_receptors and has_connections_ ) + { + throw BadProperty( + "The neuron has connections, therefore the number of ports cannot be " + "reduced." ); + } + for ( size_t i = 0; i < tau_syn_.size(); ++i ) + { + if ( tau_syn_[ i ] <= 0 ) + { + throw BadProperty( "All synaptic time constants must be strictly positive." ); + } + } + } + if ( updateValue< std::vector< double > >( d, names::tau_syn_slow, tau_syn_slow_ ) ) + { + for ( size_t i = 0; i < tau_syn_slow_.size(); ++i ) + { + if ( tau_syn_slow_[ i ] <= 0 ) + { + throw BadProperty( "All slow synaptic time constants must be strictly positive." ); + } + } + } + // for amp_slow + if ( updateValue< std::vector< double > >( d, names::amp_slow, amp_slow_ ) ) + { + for ( size_t i = 0; i < amp_slow_.size(); ++i ) + { + if ( amp_slow_[ i ] <= 0 ) + { + throw BadProperty( "All slow synaptic amplitudes must be strictly positive." ); + } + } + } + + return delta_EL; +} + +void +nest::glif_psc_double_alpha::State_::get( DictionaryDatum& d, const Parameters_& p ) const +{ + def< double >( d, names::V_m, U_ + p.E_L_ ); + def< std::vector< double > >( d, names::ASCurrents, ASCurrents_ ); + def< double >( d, names::threshold_spike, threshold_spike_ ); + def< double >( d, names::threshold_voltage, threshold_voltage_ ); +} + +void +nest::glif_psc_double_alpha::State_::set( const DictionaryDatum& d, const Parameters_& p, double delta_EL, Node* node ) +{ + if ( updateValueParam< double >( d, names::V_m, U_, node ) ) + { + U_ -= p.E_L_; + } + else + { + U_ -= delta_EL; + } + + bool asc_flag = updateValue< std::vector< double > >( d, names::ASCurrents, ASCurrents_ ); + if ( asc_flag and not p.has_asc_ ) + { + throw BadProperty( "After spike currents are not supported or settable in the current model mechanisms." ); + } + + const size_t asc_size = p.asc_decay_.size(); + if ( asc_flag ) + { + if ( ASCurrents_.size() != asc_size ) + { + throw BadProperty( "After spike current values must have have the same size (" + std::to_string( asc_size ) + + ") of its parameters (i.e., asc_init, k, asc_amps, r)." ); + } + } + + if ( updateValueParam< double >( d, names::threshold_spike, threshold_spike_, node ) and not p.has_theta_spike_ ) + { + throw BadProperty( "Threshold spike component is not supported or settable in the current model mechanisms." ); + } + + if ( updateValueParam< double >( d, names::threshold_voltage, threshold_voltage_, node ) + and not p.has_theta_voltage_ ) + { + throw BadProperty( "Threshold voltage component is not supported or settable in the current model mechanisms." ); + } +} + +nest::glif_psc_double_alpha::Buffers_::Buffers_( glif_psc_double_alpha& n ) + : logger_( n ) +{ +} + +nest::glif_psc_double_alpha::Buffers_::Buffers_( const Buffers_&, glif_psc_double_alpha& n ) + : logger_( n ) +{ +} + +/* ---------------------------------------------------------------- + * Default and copy constructor for node + * ---------------------------------------------------------------- */ + +nest::glif_psc_double_alpha::glif_psc_double_alpha() + : ArchivingNode() + , P_() + , S_( P_ ) + , B_( *this ) +{ + recordablesMap_.create(); +} + +nest::glif_psc_double_alpha::glif_psc_double_alpha( const glif_psc_double_alpha& n ) + : ArchivingNode( n ) + , P_( n.P_ ) + , S_( n.S_ ) + , B_( n.B_, *this ) +{ +} + +/* ---------------------------------------------------------------- + * Node initialization functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::init_buffers_() +{ + B_.spikes_.clear(); // includes resize + B_.currents_.clear(); // include resize + B_.logger_.reset(); // includes resize +} + +void +nest::glif_psc_double_alpha::pre_run_hook() +{ + B_.logger_.init(); + + const double h = Time::get_resolution().get_ms(); // in ms + + // pre-computing of decay parameters + if ( P_.has_theta_spike_ ) + { + V_.theta_spike_decay_rate_ = std::exp( -P_.th_spike_decay_ * h ); + V_.theta_spike_refractory_decay_rate_ = std::exp( -P_.th_spike_decay_ * P_.t_ref_ ); + } + + if ( P_.has_asc_ ) + { + V_.asc_decay_rates_.resize( P_.asc_decay_.size() ); + V_.asc_stable_coeff_.resize( P_.asc_decay_.size() ); + V_.asc_refractory_decay_rates_.resize( P_.asc_decay_.size() ); + for ( std::size_t a = 0; a < P_.asc_decay_.size(); ++a ) + { + V_.asc_decay_rates_[ a ] = std::exp( -P_.asc_decay_[ a ] * h ); + V_.asc_stable_coeff_[ a ] = ( ( 1.0 / P_.asc_decay_[ a ] ) / h ) * ( 1.0 - V_.asc_decay_rates_[ a ] ); + V_.asc_refractory_decay_rates_[ a ] = P_.asc_r_[ a ] * std::exp( -P_.asc_decay_[ a ] * P_.t_ref_ ); + } + } + + if ( P_.has_theta_voltage_ ) + { + V_.potential_decay_rate_ = std::exp( -P_.G_ * h / P_.C_m_ ); + V_.theta_voltage_decay_rate_inverse_ = 1 / ( std::exp( P_.th_voltage_decay_ * h ) ); + V_.phi = P_.th_voltage_index_ / ( P_.th_voltage_decay_ - P_.G_ / P_.C_m_ ); + V_.abpara_ratio_voltage_ = P_.th_voltage_index_ / P_.th_voltage_decay_; + } + + // postsynaptic currents + V_.P11_.resize( P_.n_receptors_() ); + V_.P21_.resize( P_.n_receptors_() ); + V_.P22_.resize( P_.n_receptors_() ); + V_.P31_.resize( P_.n_receptors_() ); + V_.P32_.resize( P_.n_receptors_() ); + + // Additional variables for the slow component + V_.P11_slow_.resize( P_.n_receptors_() ); + V_.P21_slow_.resize( P_.n_receptors_() ); + V_.P22_slow_.resize( P_.n_receptors_() ); + V_.P31_slow_.resize( P_.n_receptors_() ); + V_.P32_slow_.resize( P_.n_receptors_() ); + + S_.y1_.resize( P_.n_receptors_() ); + S_.y2_.resize( P_.n_receptors_() ); + S_.y1_slow_.resize( P_.n_receptors_() ); // Slow component + S_.y2_slow_.resize( P_.n_receptors_() ); // Slow component + V_.PSCInitialValues_.resize( P_.n_receptors_() ); + V_.PSCInitialValues_slow_.resize( P_.n_receptors_() ); // Slow component + + B_.spikes_.resize( P_.n_receptors_() ); + + double Tau_ = P_.C_m_ / P_.G_; // in ms + V_.P33_ = std::exp( -h / Tau_ ); + V_.P30_ = 1 / P_.C_m_ * ( 1 - V_.P33_ ) * Tau_; + + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + // these P are independent + V_.P11_[ i ] = V_.P22_[ i ] = std::exp( -h / P_.tau_syn_[ i ] ); + + V_.P21_[ i ] = h * V_.P11_[ i ]; + + // these are determined according to a numeric stability criterion + // input time parameter shall be in ms, capacity in pF + std::tie( V_.P31_[ i ], V_.P32_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + + // Slow component + V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); + V_.P21_slow_[ i ] = h * V_.P11_slow_[ i ]; + std::tie( V_.P31_slow_[ i ], V_.P32_slow_[ i ] ) = + IAFPropagatorAlpha( P_.tau_syn_slow_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + + + V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_[ i ]; + V_.PSCInitialValues_slow_[ i ] = 1.0 * numerics::e / P_.tau_syn_slow_[ i ] * P_.amp_slow_[ i ]; + B_.spikes_[ i ].resize(); + } + + V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); +} + +/* ---------------------------------------------------------------- + * Update and spike handling functions + * ---------------------------------------------------------------- */ + +void +nest::glif_psc_double_alpha::update( Time const& origin, const long from, const long to ) +{ + + double v_old = S_.U_; + + for ( long lag = from; lag < to; ++lag ) + { + + if ( S_.refractory_steps_ == 0 ) + { + // neuron not refractory, integrate voltage and currents + + // update threshold via exact solution of dynamics of spike component of + // threshold for glif2/4/5 models with "R" + if ( P_.has_theta_spike_ ) + { + S_.threshold_spike_ = S_.threshold_spike_ * V_.theta_spike_decay_rate_; + } + + // Calculate new ASCurrents value using exponential methods + S_.ASCurrents_sum_ = 0.0; + // for glif3/4/5 models with "ASC" + // take after spike current value at the beginning of the time to compute + // the exact mean ASC for the time step and sum the exact ASCs of all ports; + // and then update the current values to the value at the end of the time + // step, ready for the next time step + if ( P_.has_asc_ ) + { + for ( std::size_t a = 0; a < S_.ASCurrents_.size(); ++a ) + { + S_.ASCurrents_sum_ += ( V_.asc_stable_coeff_[ a ] * S_.ASCurrents_[ a ] ); + S_.ASCurrents_[ a ] = S_.ASCurrents_[ a ] * V_.asc_decay_rates_[ a ]; + } + } + + // voltage dynamics of membranes, linear exact to find next V_m value + S_.U_ = v_old * V_.P33_ + ( S_.I_ + S_.ASCurrents_sum_ ) * V_.P30_; + + // add synapse component for voltage dynamics + S_.I_syn_fast_ = 0.0; + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + S_.U_ += V_.P31_[ i ] * S_.y1_[ i ] + V_.P32_[ i ] * S_.y2_[ i ]; + S_.I_syn_fast_ += S_.y2_[ i ]; + } + // slow component + S_.I_syn_slow_ = 0.0; + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + // S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + (V_.P32_slow_[ i ] * S_.y2_slow_[ i ]) * P_.amp_slow_[ i ]; + S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + V_.P32_slow_[ i ] * S_.y2_slow_[ i ]; + // multiply amp_slow_ + S_.I_syn_slow_ += S_.y2_slow_[ i ]; + // print all components of P3x and y the slow current. + // std::cout << "P31_slow: " << V_.P31_slow_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << " y2_slow:" << + // S_.y2_slow_[ i ] << std::endl; + } + // add the fast and slow components + S_.I_syn_ = S_.I_syn_fast_ + S_.I_syn_slow_; + // print to diagnose the synaptic current + // std::cout << "I_syn_fast: " << S_.I_syn_fast_ << " I_syn_slow: " << S_.I_syn_slow_ << std::endl; + + // Calculate exact voltage component of the threshold for glif5 model with + // "A" + if ( P_.has_theta_voltage_ ) + { + const double beta = ( S_.I_ + S_.ASCurrents_sum_ ) / P_.G_; + S_.threshold_voltage_ = V_.phi * ( v_old - beta ) * V_.potential_decay_rate_ + + V_.theta_voltage_decay_rate_inverse_ + * ( S_.threshold_voltage_ - V_.phi * ( v_old - beta ) - V_.abpara_ratio_voltage_ * beta ) + + V_.abpara_ratio_voltage_ * beta; + } + + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + + // Check if there is an action potential + if ( S_.U_ > S_.threshold_ ) + { + // Marks that the neuron is in a refractory period + S_.refractory_steps_ = V_.RefractoryCounts_; + + // Reset ASC_currents for glif3/4/5 models with "ASC" + if ( P_.has_asc_ ) + { + for ( std::size_t a = 0; a < S_.ASCurrents_.size(); ++a ) + { + S_.ASCurrents_[ a ] = P_.asc_amps_[ a ] + S_.ASCurrents_[ a ] * V_.asc_refractory_decay_rates_[ a ]; + } + } + + // Reset voltage + if ( not P_.has_theta_spike_ ) + { + // Reset voltage for glif1/3 models without "R" + S_.U_ = P_.V_reset_; + } + else + { + // Reset voltage for glif2/4/5 models with "R" + S_.U_ = P_.voltage_reset_fraction_ * v_old + P_.voltage_reset_add_; + + // reset spike component of threshold + // (decay for refractory period and then add additive constant) + S_.threshold_spike_ = S_.threshold_spike_ * V_.theta_spike_refractory_decay_rate_ + P_.th_spike_add_; + + // rest the global threshold (voltage component of threshold: stay the + // same) + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + } + + set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); + SpikeEvent se; + kernel().event_delivery_manager.send( *this, se, lag ); + } + } + else + { + // neuron is absolute refractory + --S_.refractory_steps_; + + // While neuron is in refractory period count-down in time steps (since dt + // may change while in refractory) while holding the voltage at last peak. + S_.U_ = v_old; + S_.threshold_ = S_.threshold_spike_ + S_.threshold_voltage_ + P_.th_inf_; + } + + // alpha shape PSCs + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) + { + double spike_value = B_.spikes_[ i ].get_value( lag ); + // Apply spikes delivered in this step: The spikes arriving at T+1 have an + // immediate effect on the state of the neuron + // store B_.spikes_[ i ].get_value( lag ); into a variable + S_.y2_[ i ] = V_.P21_[ i ] * S_.y1_[ i ] + V_.P22_[ i ] * S_.y2_[ i ]; + S_.y1_[ i ] *= V_.P11_[ i ]; + S_.y1_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; + + // slow component + S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; + S_.y1_slow_[ i ] *= V_.P11_slow_[ i ]; + S_.y1_slow_[ i ] += V_.PSCInitialValues_slow_[ i ] * spike_value; + // print y1 and y1 slow + // std::cout << "y1: " << S_.y1_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << std::endl; + // print PSCInitlaValues + // std::cout << "PSCInitialValues: " << V_.PSCInitialValues_[ i ] << " PSCInitialValues_slow: " << + // V_.PSCInitialValues_slow_[ i ] << std::endl; + } + + // Update any external currents + S_.I_ = B_.currents_.get_value( lag ); + + // Save voltage + B_.logger_.record_data( origin.get_steps() + lag ); + v_old = S_.U_; + } +} + +size_t +nest::glif_psc_double_alpha::handles_test_event( SpikeEvent&, size_t receptor_type ) +{ + if ( receptor_type <= 0 or receptor_type > P_.n_receptors_() ) + { + throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" ); + } + + P_.has_connections_ = true; + return receptor_type; +} + +void +nest::glif_psc_double_alpha::handle( SpikeEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + B_.spikes_[ e.get_rport() - 1 ].add_value( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), e.get_weight() * e.get_multiplicity() ); +} + +void +nest::glif_psc_double_alpha::handle( CurrentEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + B_.currents_.add_value( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), e.get_weight() * e.get_current() ); +} + +// Do not move this function as inline to h-file. It depends on +// universal_data_logger_impl.h being included here. +void +nest::glif_psc_double_alpha::handle( DataLoggingRequest& e ) +{ + B_.logger_.handle( e ); // the logger does this for us +} diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h new file mode 100644 index 0000000000..2b6336fae3 --- /dev/null +++ b/models/glif_psc_double_alpha.h @@ -0,0 +1,494 @@ +/* + * glif_psc_double_alpha.h + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#ifndef GLIF_PSC_DOUBLE_ALPHA_H +#define GLIF_PSC_DOUBLE_ALPHA_H + +#include "archiving_node.h" +#include "connection.h" +#include "event.h" +#include "nest_types.h" +#include "ring_buffer.h" +#include "universal_data_logger.h" + +#include "dictdatum.h" + +/* BeginUserDocs: integrate-and-fire, current-based + +Short description ++++++++++++++++++ + +Current-based generalized leaky integrate-and-fire (GLIF) models (from the Allen Institute) + +Description ++++++++++++ + +``glif_psc_double_alpha`` provides five generalized leaky integrate-and-fire +(GLIF) models [1]_ with alpha-function shaped synaptic currents. +Incoming spike events induce a postsynaptic change of current modeled +by sum of two alpha functions [2]_. The alpha function is normalized such that an event +of weight 1.0 results in a peak current of the fast component of the alpha function to +be 1 pA at :math:`t = tau_syn`. The peak current of the slow component of the alpha is +given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha has +a single synapse that is accessible through receptor_port 1. An arbitrary number of +synapses with different time constants can be configured by setting the desired time +constants as tau_syn array. The resulting synapses are addressed through receptor_port +1, 2, 3, .... + +The five GLIF models are: + +* **GLIF Model 1** - Traditional leaky integrate and fire (LIF) +* **GLIF Model 2** - Leaky integrate and fire with biologically defined reset rules + (LIF_R) +* **GLIF Model 3** - Leaky integrate and fire with after-spike currents (LIF_ASC) +* **GLIF Model 4** - Leaky integrate and fire with biologically defined reset rules + and after-spike currents (LIF_R_ASC) +* **GLIF Model 5** - Leaky integrate and fire with biologically defined reset rules, + after-spike currents and a voltage dependent threshold (LIF_R_ASC_A) + +GLIF model mechanism setting is based on three parameters +(``spike_dependent_threshold``, ``after_spike_currents``, ``adapting_threshold``). +The settings of these three parameters for the five GLIF models are listed +below. Other combinations of these parameters will not be supported. + ++--------+---------------------------+----------------------+--------------------+ +| Model | spike_dependent_threshold | after_spike_currents | adapting_threshold | ++========+===========================+======================+====================+ +| GLIF1 | False | False | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF2 | True | False | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF3 | False | True | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF4 | True | True | False | ++--------+---------------------------+----------------------+--------------------+ +| GLIF5 | True | True | True | ++--------+---------------------------+----------------------+--------------------+ + +Typical parameter setting of different levels of GLIF models for different cells +can be found and downloaded in the `Allen Cell Type Database +`_. For example, the default parameter setting of this +``glif_psc_double_alpha`` neuron model was from the parameter values of GLIF Model 5 of +Cell 490626718, which can be retrieved from the `Allen Brain Atlas +`_, with units being converted from SI units (i.e., V, S (1/Ohm), +F, s, A) to NEST used units (i.e., mV, nS (1/GOhm), pF, ms, pA) and values +being rounded to appropriate digits for simplification. + +For models with spike dependent threshold (i.e., GLIF2, GLIF4 and GLIF5), +parameter setting of voltage_reset_fraction and voltage_reset_add may lead to the +situation that voltage is bigger than threshold after reset. In this case, the neuron +will continue to spike until the end of the simulation regardless the stimulated inputs. +We recommend the setting of the parameters of these three models to follow the +condition of + +.. math:: + + E_L + \mathrm{voltage\_reset\_fraction} \cdot \left( V_\mathrm{th} - E_L \right) + + \mathrm{voltage\_reset\_add} < V_\mathrm{th} + \mathrm{th\_spike\_add} + +.. note:: + + If ``tau_m`` is very close to ``tau_syn_ex`` or ``tau_syn_in``, the model + will numerically behave as if ``tau_m`` is equal to ``tau_syn_ex`` or + ``tau_syn_in``, respectively, to avoid numerical instabilities. + + For implementation details see the + `IAF_neurons_singularity <../model_details/IAF_neurons_singularity.ipynb>`_ notebook. + +Parameters +++++++++++ + +The following parameters can be set in the status dictionary. + +========= ======== ============================================================ +**Membrane parameters** +------------------------------------------------------------------------------- +V_m double Membrane potential in mV (absolute value) +V_th double Instantaneous threshold in mV +g double Membrane conductance in nS +E_L double Resting membrane potential in mV +C_m double Capacitance of the membrane in pF +t_ref double Duration of refractory time in ms +V_reset double Reset potential of the membrane in mV (GLIF 1 or GLIF 3) +========= ======== ============================================================ + +========================= =============== ===================================== +**Spike adaptation and firing intensity parameters** +------------------------------------------------------------------------------- +th_spike_add double Threshold addition following spike + in mV (delta_theta_s in Equation (6) + in [1]_) +th_spike_decay double Spike-induced threshold time + constant in 1/ms (bs in Equation (2) + in [1]_) +voltage_reset_fraction double Voltage fraction coefficient + following spike (fv in Equation (5) + in [1]_) +voltage_reset_add double Voltage addition following spike in + mV (-delta_V (sign flipped) in + Equation (5) in [1]_) +asc_init double vector Initial values of after-spike + currents in pA +asc_decay double vector After-spike current time constants + in 1/ms (kj in Equation (3) in [1]_) +asc_amps double vector After-spike current amplitudes in + pA (deltaIj in Equation (7) in [1]_) +asc_r double vector Current fraction following spike + coefficients for fj in Equation (7) + in [1]_ +th_voltage_index double Adaptation index of threshold - A + 'leak-conductance' for the + voltage-dependent component of the + threshold in 1/ms (av in Equation + (4) in [1]_) +th_voltage_decay double Voltage-induced threshold time + constant - Inverse of which is the + time constant of the + voltage-dependent component of the + threshold in 1/ms (bv in Equation + (4) in [1]_) +tau_syn double vector Rise time constants of the synaptic + alpha function in ms +tau_syn_slow double vector Rise time constants of the slower + synaptic alpha function in ms +amp_slow double vector Relative amplitude of the slower + synaptic alpha function +E_rev double vector Reversal potential in mV +spike_dependent_threshold bool flag whether the neuron has + biologically defined reset rules + with a spike dependent threshold + component +after_spike_currents bool flag whether the neuron has after + spike currents +adapting_threshold bool flag whether the neuron has a + voltage dependent threshold component +========================= =============== ===================================== + +References +++++++++++ + +.. [1] Teeter C, Iyer R, Menon V, Gouwens N, Feng D, Berg J, Szafer A, + Cain N, Zeng H, Hawrylycz M, Koch C, & Mihalas S (2018) + Generalized leaky integrate-and-fire models classify multiple neuron + types. Nature Communications 9:709. +.. [2] Meffin, H., Burkitt, A. N., & Grayden, D. B. (2004). An analytical + model for the large, fluctuating synaptic conductance state typical of + neocortical neurons in vivo. J. Comput. Neurosci., 16, 159-175. + +See also +++++++++ + +gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp, +glif_psc + +EndUserDocs */ + +namespace nest +{ + +class glif_psc_double_alpha : public ArchivingNode +{ +public: + glif_psc_double_alpha(); + + glif_psc_double_alpha( const glif_psc_double_alpha& ); + + using nest::Node::handle; + using nest::Node::handles_test_event; + + size_t send_test_event( nest::Node&, size_t, nest::synindex, bool ) override; + + void handle( nest::SpikeEvent& ) override; + void handle( nest::CurrentEvent& ) override; + void handle( nest::DataLoggingRequest& ) override; + + size_t handles_test_event( nest::SpikeEvent&, size_t ) override; + size_t handles_test_event( nest::CurrentEvent&, size_t ) override; + size_t handles_test_event( nest::DataLoggingRequest&, size_t ) override; + + void get_status( DictionaryDatum& ) const override; + void set_status( const DictionaryDatum& ) override; + +private: + //! Reset internal buffers of neuron. + void init_buffers_() override; + + //! Initialize auxiliary quantities, leave parameters and state untouched. + void pre_run_hook() override; + + //! Take neuron through given time interval + void update( nest::Time const&, const long, const long ) override; + + // The next two classes need to be friends to access the State_ class/member + friend class nest::RecordablesMap< glif_psc_double_alpha >; + friend class nest::UniversalDataLogger< glif_psc_double_alpha >; + + struct Parameters_ + { + double G_; //!< membrane conductance in nS + double E_L_; //!< resting potential in mV + double th_inf_; //!< infinity threshold in mV + double C_m_; //!< capacitance in pF + double t_ref_; //!< refractory time in ms + double V_reset_; //!< Membrane voltage following spike in mV + double th_spike_add_; //!< threshold additive constant following reset in mV + double th_spike_decay_; //!< spike induced threshold in 1/ms + double voltage_reset_fraction_; //!< voltage fraction following reset coefficient + double voltage_reset_add_; //!< voltage additive constant following reset in mV + double th_voltage_index_; //!< a 'leak-conductance' for the voltage-dependent + //!< component of the threshold in 1/ms + double th_voltage_decay_; //!< inverse of which is the time constant of the + //!< voltage-dependent component of the threshold in 1/ms + std::vector< double > asc_init_; //!< initial values of ASCurrents_ in pA + std::vector< double > asc_decay_; //!< predefined time scale in 1/ms + std::vector< double > asc_amps_; //!< in pA + std::vector< double > asc_r_; //!< coefficient + std::vector< double > tau_syn_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms + std::vector< double > amp_slow_; //!< synaptic port time constants in ms + + //! boolean flag which indicates whether the neuron has connections + bool has_connections_; + + //! boolean flag which indicates whether the neuron has spike dependent threshold component + bool has_theta_spike_; + + //! boolean flag which indicates whether the neuron has after spike currents + bool has_asc_; + + //! boolean flag which indicates whether the neuron has voltage dependent threshold component + bool has_theta_voltage_; + + size_t n_receptors_() const; //!< Returns the size of tau_syn_ + + Parameters_(); + + void get( DictionaryDatum& ) const; + double set( const DictionaryDatum&, Node* ); + }; + + struct State_ + { + double U_; //!< relative membrane potential in mV + double threshold_; //!< total threshold in mV + double threshold_spike_; //!< spike component of threshold in mV + double threshold_voltage_; //!< voltage component of threshold in mV + double I_; //!< external current in pA + double I_syn_; //!< postsynaptic current in pA + double I_syn_fast_; //!< postsynaptic current in pA + double I_syn_slow_; //!< postsynaptic current in pA + std::vector< double > ASCurrents_; //!< after-spike currents in pA + double ASCurrents_sum_; //!< in pA + int refractory_steps_; //!< Number of refractory steps remaining + std::vector< double > y1_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA + + State_( const Parameters_& ); + + void get( DictionaryDatum&, const Parameters_& ) const; + void set( const DictionaryDatum&, const Parameters_&, double, Node* ); + }; + + + struct Buffers_ + { + Buffers_( glif_psc_double_alpha& ); + Buffers_( const Buffers_&, glif_psc_double_alpha& ); + + std::vector< nest::RingBuffer > spikes_; //!< Buffer incoming spikes through delay, as sum + nest::RingBuffer currents_; //!< Buffer incoming currents through delay, + + //! Logger for all analog data + nest::UniversalDataLogger< glif_psc_double_alpha > logger_; + }; + + struct Variables_ + { + int RefractoryCounts_; //!< counter during refractory period + double theta_spike_decay_rate_; //!< threshold spike component decay rate + double theta_spike_refractory_decay_rate_; //!< threshold spike component decay rate during refractory + double theta_voltage_decay_rate_inverse_; //!< inverse of threshold voltage component decay rate + double potential_decay_rate_; //!< membrane potential decay rate + double abpara_ratio_voltage_; //!< ratio of parameters of voltage threshold component av/bv + std::vector< double > asc_decay_rates_; //!< after spike current decay rates + std::vector< double > asc_stable_coeff_; //!< after spike current stable coefficient + std::vector< double > asc_refractory_decay_rates_; //!< after spike current decay rates during refractory + double phi; //!< threshold voltage component coefficient + + std::vector< double > P11_; //!< synaptic current evolution parameter + std::vector< double > P21_; //!< synaptic current evolution parameter + std::vector< double > P22_; //!< synaptic current evolution parameter + double P30_; //!< membrane current/voltage evolution parameter + double P33_; //!< membrane voltage evolution parameter + std::vector< double > P31_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_; //!< synaptic/membrane current evolution parameter + + // slow component + std::vector< double > P11_slow_; //!< synaptic current evolution parameter + std::vector< double > P21_slow_; //!< synaptic current evolution parameter + std::vector< double > P22_slow_; //!< synaptic current evolution parameter + double P30_slow_; //!< membrane current/voltage evolution parameter + double P33_slow_; //!< membrane voltage evolution parameter + std::vector< double > P31_slow_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_slow_; //!< synaptic/membrane current evolution parameter + + + /** Amplitude of the synaptic current. + This value is chosen such that a postsynaptic current with + weight one has an amplitude of 1 pA. + */ + std::vector< double > PSCInitialValues_; + std::vector< double > PSCInitialValues_slow_; + }; + + double + get_V_m_() const + { + return S_.U_ + P_.E_L_; + } + + double + get_ASCurrents_sum_() const + { + return S_.ASCurrents_sum_; + } + + double + get_I_() const + { + return S_.I_; + } + + double + get_I_syn_() const + { + return S_.I_syn_; + } + + // double + // get_I_syn_fast_() const + // { + // return S_.I_syn_fast_; + // } + + // double + // get_I_syn_slow_() const + // { + // return S_.I_syn_slow_; + // } + + double + get_threshold_() const + { + return S_.threshold_ + P_.E_L_; + } + + double + get_threshold_spike_() const + { + return S_.threshold_spike_; + } + + double + get_threshold_voltage_() const + { + return S_.threshold_voltage_; + } + + Parameters_ P_; + State_ S_; + Variables_ V_; + Buffers_ B_; + + //! Mapping of recordables names to access functions + static nest::RecordablesMap< glif_psc_double_alpha > recordablesMap_; +}; + + +inline size_t +nest::glif_psc_double_alpha::Parameters_::n_receptors_() const +{ + return tau_syn_.size(); +} + +inline size_t +nest::glif_psc_double_alpha::send_test_event( nest::Node& target, size_t receptor_type, nest::synindex, bool ) +{ + nest::SpikeEvent e; + e.set_sender( *this ); + return target.handles_test_event( e, receptor_type ); +} + +inline size_t +nest::glif_psc_double_alpha::handles_test_event( nest::CurrentEvent&, size_t receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw nest::UnknownReceptorType( receptor_type, get_name() ); + } + return 0; +} + +inline size_t +nest::glif_psc_double_alpha::handles_test_event( nest::DataLoggingRequest& dlr, size_t receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw nest::UnknownReceptorType( receptor_type, get_name() ); + } + return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); +} + +inline void +glif_psc_double_alpha::get_status( DictionaryDatum& d ) const +{ + // get our own parameter and state data + P_.get( d ); + S_.get( d, P_ ); + + // get information managed by parent class + ArchivingNode::get_status( d ); + + ( *d )[ nest::names::recordables ] = recordablesMap_.get_list(); +} + +inline void +glif_psc_double_alpha::set_status( const DictionaryDatum& d ) +{ + Parameters_ ptmp = P_; // temporary copy in case of errors + const double delta_EL = ptmp.set( d, this ); // throws if BadProperty + State_ stmp = S_; // temporary copy in case of errors + stmp.set( d, ptmp, delta_EL, this ); // throws if BadProperty + + ArchivingNode::set_status( d ); + + // if we get here, temporaries contain consistent set of properties + P_ = ptmp; + S_ = stmp; +} + +} // namespace nest + +#endif diff --git a/modelsets/full b/modelsets/full index 4f0e835f41..9f7ff8bbed 100644 --- a/modelsets/full +++ b/modelsets/full @@ -32,6 +32,7 @@ gif_pop_psc_exp ginzburg_neuron glif_cond glif_psc +glif_psc_double_alpha hh_cond_exp_traub hh_cond_beta_gap_traub hh_psc_alpha diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 7559f5279b..9fe6b9a993 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -62,6 +62,7 @@ const Name allow_oversized_mask( "allow_oversized_mask" ); const Name alpha( "alpha" ); const Name alpha_1( "alpha_1" ); const Name alpha_2( "alpha_2" ); +const Name amp_slow( "amp_slow" ); const Name amplitude( "amplitude" ); const Name amplitude_times( "amplitude_times" ); const Name amplitude_values( "amplitude_values" ); @@ -502,6 +503,7 @@ const Name tau_stc( "tau_stc" ); const Name tau_syn( "tau_syn" ); const Name tau_syn_ex( "tau_syn_ex" ); const Name tau_syn_in( "tau_syn_in" ); +const Name tau_syn_slow( "tau_syn_slow" ); const Name tau_theta( "tau_theta" ); const Name tau_u_bar_bar( "tau_u_bar_bar" ); const Name tau_u_bar_minus( "tau_u_bar_minus" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index cc2aa1d516..dc34d576a0 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -87,6 +87,7 @@ extern const Name allow_oversized_mask; extern const Name alpha; extern const Name alpha_1; extern const Name alpha_2; +extern const Name amp_slow; extern const Name amplitude; extern const Name amplitude_times; extern const Name amplitude_values; @@ -527,6 +528,7 @@ extern const Name tau_stc; extern const Name tau_syn; extern const Name tau_syn_ex; extern const Name tau_syn_in; +extern const Name tau_syn_slow; extern const Name tau_theta; extern const Name tau_u_bar_bar; extern const Name tau_u_bar_minus; From b09fdda7281293a898b8337e2808ad153139cdf4 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 1 Aug 2023 10:48:57 -0700 Subject: [PATCH 23/45] Cleaning up commented out codes (They were used during development.) --- models/glif_psc_double_alpha.cpp | 12 ------------ models/glif_psc_double_alpha.h | 12 ------------ 2 files changed, 24 deletions(-) diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index b48b75d7b4..10011da27b 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -564,18 +564,11 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.I_syn_slow_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { - // S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + (V_.P32_slow_[ i ] * S_.y2_slow_[ i ]) * P_.amp_slow_[ i ]; S_.U_ += V_.P31_slow_[ i ] * S_.y1_slow_[ i ] + V_.P32_slow_[ i ] * S_.y2_slow_[ i ]; - // multiply amp_slow_ S_.I_syn_slow_ += S_.y2_slow_[ i ]; - // print all components of P3x and y the slow current. - // std::cout << "P31_slow: " << V_.P31_slow_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << " y2_slow:" << - // S_.y2_slow_[ i ] << std::endl; } // add the fast and slow components S_.I_syn_ = S_.I_syn_fast_ + S_.I_syn_slow_; - // print to diagnose the synaptic current - // std::cout << "I_syn_fast: " << S_.I_syn_fast_ << " I_syn_slow: " << S_.I_syn_slow_ << std::endl; // Calculate exact voltage component of the threshold for glif5 model with // "A" @@ -656,11 +649,6 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; S_.y1_slow_[ i ] *= V_.P11_slow_[ i ]; S_.y1_slow_[ i ] += V_.PSCInitialValues_slow_[ i ] * spike_value; - // print y1 and y1 slow - // std::cout << "y1: " << S_.y1_[ i ] << " y1_slow: " << S_.y1_slow_[ i ] << std::endl; - // print PSCInitlaValues - // std::cout << "PSCInitialValues: " << V_.PSCInitialValues_[ i ] << " PSCInitialValues_slow: " << - // V_.PSCInitialValues_slow_[ i ] << std::endl; } // Update any external currents diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2b6336fae3..06bbc08376 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -387,18 +387,6 @@ class glif_psc_double_alpha : public ArchivingNode return S_.I_syn_; } - // double - // get_I_syn_fast_() const - // { - // return S_.I_syn_fast_; - // } - - // double - // get_I_syn_slow_() const - // { - // return S_.I_syn_slow_; - // } - double get_threshold_() const { From ae508728d6279d0d7d5d6f781a3eb02f3b4f631d Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 15 Aug 2023 17:30:04 -0700 Subject: [PATCH 24/45] Variable name correction and adding tests Some of the parameters that are associated with the fast synaptic components get '_fast' added to their names. Pytest script added to test the GLIF model with double alpha synapses. A couple of test cases require adding extra parameters for the new glif_psc_double_alpha model, so they are added. --- models/glif_psc_double_alpha.cpp | 48 +-- models/glif_psc_double_alpha.h | 39 ++- nestkernel/nest_names.cpp | 1 + nestkernel/nest_names.h | 1 + .../test_neurons_handle_multiplicity.py | 3 + .../test_multimeter_stepping.py | 1 + .../pytests/test_glif_psc_double_alpha.py | 305 ++++++++++++++++++ 7 files changed, 354 insertions(+), 44 deletions(-) create mode 100644 testsuite/pytests/test_glif_psc_double_alpha.py diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index 10011da27b..562154488e 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -80,7 +80,7 @@ nest::glif_psc_double_alpha::Parameters_::Parameters_() , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA , asc_r_( std::vector< double >( 2, 1.0 ) ) - , tau_syn_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless , has_connections_( false ) @@ -107,8 +107,8 @@ nest::glif_psc_double_alpha::State_::State_( const Parameters_& p ) { ASCurrents_sum_ += ASCurrents_[ a ]; } - y1_.clear(); - y2_.clear(); + y1_fast_.clear(); + y2_fast_.clear(); y1_slow_.clear(); y2_slow_.clear(); } @@ -140,8 +140,8 @@ nest::glif_psc_double_alpha::Parameters_::get( DictionaryDatum& d ) const def< std::vector< double > >( d, names::asc_amps, asc_amps_ ); def< std::vector< double > >( d, names::asc_r, asc_r_ ); - ArrayDatum tau_syn_ad( tau_syn_ ); - def< ArrayDatum >( d, names::tau_syn, tau_syn_ad ); + ArrayDatum tau_syn_fast_ad( tau_syn_fast_ ); + def< ArrayDatum >( d, names::tau_syn_fast, tau_syn_fast_ad ); ArrayDatum tau_syn_slow_ad( tau_syn_slow_ ); def< ArrayDatum >( d, names::tau_syn_slow, tau_syn_slow_ad ); @@ -289,7 +289,7 @@ nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* n } const size_t old_n_receptors = this->n_receptors_(); - if ( updateValue< std::vector< double > >( d, names::tau_syn, tau_syn_ ) ) + if ( updateValue< std::vector< double > >( d, names::tau_syn_fast, tau_syn_fast_ ) ) { if ( this->n_receptors_() != old_n_receptors and has_connections_ ) { @@ -297,9 +297,9 @@ nest::glif_psc_double_alpha::Parameters_::set( const DictionaryDatum& d, Node* n "The neuron has connections, therefore the number of ports cannot be " "reduced." ); } - for ( size_t i = 0; i < tau_syn_.size(); ++i ) + for ( size_t i = 0; i < tau_syn_fast_.size(); ++i ) { - if ( tau_syn_[ i ] <= 0 ) + if ( tau_syn_fast_[ i ] <= 0 ) { throw BadProperty( "All synaptic time constants must be strictly positive." ); } @@ -458,11 +458,11 @@ nest::glif_psc_double_alpha::pre_run_hook() } // postsynaptic currents - V_.P11_.resize( P_.n_receptors_() ); - V_.P21_.resize( P_.n_receptors_() ); - V_.P22_.resize( P_.n_receptors_() ); - V_.P31_.resize( P_.n_receptors_() ); - V_.P32_.resize( P_.n_receptors_() ); + V_.P11_fast_.resize( P_.n_receptors_() ); + V_.P21_fast_.resize( P_.n_receptors_() ); + V_.P22_fast_.resize( P_.n_receptors_() ); + V_.P31_fast_.resize( P_.n_receptors_() ); + V_.P32_fast_.resize( P_.n_receptors_() ); // Additional variables for the slow component V_.P11_slow_.resize( P_.n_receptors_() ); @@ -471,8 +471,8 @@ nest::glif_psc_double_alpha::pre_run_hook() V_.P31_slow_.resize( P_.n_receptors_() ); V_.P32_slow_.resize( P_.n_receptors_() ); - S_.y1_.resize( P_.n_receptors_() ); - S_.y2_.resize( P_.n_receptors_() ); + S_.y1_fast_.resize( P_.n_receptors_() ); + S_.y2_fast_.resize( P_.n_receptors_() ); S_.y1_slow_.resize( P_.n_receptors_() ); // Slow component S_.y2_slow_.resize( P_.n_receptors_() ); // Slow component V_.PSCInitialValues_.resize( P_.n_receptors_() ); @@ -487,13 +487,13 @@ nest::glif_psc_double_alpha::pre_run_hook() for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { // these P are independent - V_.P11_[ i ] = V_.P22_[ i ] = std::exp( -h / P_.tau_syn_[ i ] ); + V_.P11_fast_[ i ] = V_.P22_fast_[ i ] = std::exp( -h / P_.tau_syn_fast_[ i ] ); - V_.P21_[ i ] = h * V_.P11_[ i ]; + V_.P21_fast_[ i ] = h * V_.P11_fast_[ i ]; // these are determined according to a numeric stability criterion // input time parameter shall be in ms, capacity in pF - std::tie( V_.P31_[ i ], V_.P32_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); // Slow component V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); @@ -502,7 +502,7 @@ nest::glif_psc_double_alpha::pre_run_hook() IAFPropagatorAlpha( P_.tau_syn_slow_[ i ], Tau_, P_.C_m_ ).evaluate( h ); - V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_[ i ]; + V_.PSCInitialValues_[ i ] = 1.0 * numerics::e / P_.tau_syn_fast_[ i ]; V_.PSCInitialValues_slow_[ i ] = 1.0 * numerics::e / P_.tau_syn_slow_[ i ] * P_.amp_slow_[ i ]; B_.spikes_[ i ].resize(); } @@ -557,8 +557,8 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const S_.I_syn_fast_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { - S_.U_ += V_.P31_[ i ] * S_.y1_[ i ] + V_.P32_[ i ] * S_.y2_[ i ]; - S_.I_syn_fast_ += S_.y2_[ i ]; + S_.U_ += V_.P31_fast_[ i ] * S_.y1_fast_[ i ] + V_.P32_fast_[ i ] * S_.y2_fast_[ i ]; + S_.I_syn_fast_ += S_.y2_fast_[ i ]; } // slow component S_.I_syn_slow_ = 0.0; @@ -641,9 +641,9 @@ nest::glif_psc_double_alpha::update( Time const& origin, const long from, const // Apply spikes delivered in this step: The spikes arriving at T+1 have an // immediate effect on the state of the neuron // store B_.spikes_[ i ].get_value( lag ); into a variable - S_.y2_[ i ] = V_.P21_[ i ] * S_.y1_[ i ] + V_.P22_[ i ] * S_.y2_[ i ]; - S_.y1_[ i ] *= V_.P11_[ i ]; - S_.y1_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; + S_.y2_fast_[ i ] = V_.P21_fast_[ i ] * S_.y1_fast_[ i ] + V_.P22_fast_[ i ] * S_.y2_fast_[ i ]; + S_.y1_fast_[ i ] *= V_.P11_fast_[ i ]; + S_.y1_fast_[ i ] += V_.PSCInitialValues_[ i ] * spike_value; // slow component S_.y2_slow_[ i ] = V_.P21_slow_[ i ] * S_.y1_slow_[ i ] + V_.P22_slow_[ i ] * S_.y2_slow_[ i ]; diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 06bbc08376..cdc7449529 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -37,7 +37,8 @@ Short description +++++++++++++++++ -Current-based generalized leaky integrate-and-fire (GLIF) models (from the Allen Institute) +Current-based generalized leaky integrate-and-fire (GLIF) models +(from the Allen Institute) Description +++++++++++ @@ -47,12 +48,12 @@ Description Incoming spike events induce a postsynaptic change of current modeled by sum of two alpha functions [2]_. The alpha function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to -be 1 pA at :math:`t = tau_syn`. The peak current of the slow component of the alpha is -given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha has -a single synapse that is accessible through receptor_port 1. An arbitrary number of +be 1 pA at :math:`t = tau_syn_fast`. The peak current of the slow component of the alpha +is given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha +has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different time constants can be configured by setting the desired time -constants as tau_syn array. The resulting synapses are addressed through receptor_port -1, 2, 3, .... +constants as tau_syn_fast array. The resulting synapses are addressed through +receptor_port 1, 2, 3, .... The five GLIF models are: @@ -167,8 +168,8 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn double vector Rise time constants of the synaptic - alpha function in ms +tau_syn_fast double vector Rise time constants of the faster + synaptic alpha function in ms tau_syn_slow double vector Rise time constants of the slower synaptic alpha function in ms amp_slow double vector Relative amplitude of the slower @@ -263,7 +264,7 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_decay_; //!< predefined time scale in 1/ms std::vector< double > asc_amps_; //!< in pA std::vector< double > asc_r_; //!< coefficient - std::vector< double > tau_syn_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms std::vector< double > amp_slow_; //!< synaptic port time constants in ms @@ -279,7 +280,7 @@ class glif_psc_double_alpha : public ArchivingNode //! boolean flag which indicates whether the neuron has voltage dependent threshold component bool has_theta_voltage_; - size_t n_receptors_() const; //!< Returns the size of tau_syn_ + size_t n_receptors_() const; //!< Returns the size of tau_syn_fast_ Parameters_(); @@ -300,8 +301,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > ASCurrents_; //!< after-spike currents in pA double ASCurrents_sum_; //!< in pA int refractory_steps_; //!< Number of refractory steps remaining - std::vector< double > y1_; //!< synapse current evolution state 1 in pA - std::vector< double > y2_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA @@ -337,20 +338,18 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_refractory_decay_rates_; //!< after spike current decay rates during refractory double phi; //!< threshold voltage component coefficient - std::vector< double > P11_; //!< synaptic current evolution parameter - std::vector< double > P21_; //!< synaptic current evolution parameter - std::vector< double > P22_; //!< synaptic current evolution parameter + std::vector< double > P11_fast_; //!< synaptic current evolution parameter + std::vector< double > P21_fast_; //!< synaptic current evolution parameter + std::vector< double > P22_fast_; //!< synaptic current evolution parameter double P30_; //!< membrane current/voltage evolution parameter double P33_; //!< membrane voltage evolution parameter - std::vector< double > P31_; //!< synaptic/membrane current evolution parameter - std::vector< double > P32_; //!< synaptic/membrane current evolution parameter + std::vector< double > P31_fast_; //!< synaptic/membrane current evolution parameter + std::vector< double > P32_fast_; //!< synaptic/membrane current evolution parameter // slow component std::vector< double > P11_slow_; //!< synaptic current evolution parameter std::vector< double > P21_slow_; //!< synaptic current evolution parameter std::vector< double > P22_slow_; //!< synaptic current evolution parameter - double P30_slow_; //!< membrane current/voltage evolution parameter - double P33_slow_; //!< membrane voltage evolution parameter std::vector< double > P31_slow_; //!< synaptic/membrane current evolution parameter std::vector< double > P32_slow_; //!< synaptic/membrane current evolution parameter @@ -418,7 +417,7 @@ class glif_psc_double_alpha : public ArchivingNode inline size_t nest::glif_psc_double_alpha::Parameters_::n_receptors_() const { - return tau_syn_.size(); + return tau_syn_fast_.size(); } inline size_t diff --git a/nestkernel/nest_names.cpp b/nestkernel/nest_names.cpp index 9fe6b9a993..df28532990 100644 --- a/nestkernel/nest_names.cpp +++ b/nestkernel/nest_names.cpp @@ -502,6 +502,7 @@ const Name tau_spike( "tau_spike" ); const Name tau_stc( "tau_stc" ); const Name tau_syn( "tau_syn" ); const Name tau_syn_ex( "tau_syn_ex" ); +const Name tau_syn_fast( "tau_syn_fast" ); const Name tau_syn_in( "tau_syn_in" ); const Name tau_syn_slow( "tau_syn_slow" ); const Name tau_theta( "tau_theta" ); diff --git a/nestkernel/nest_names.h b/nestkernel/nest_names.h index dc34d576a0..d659ab2c1d 100644 --- a/nestkernel/nest_names.h +++ b/nestkernel/nest_names.h @@ -527,6 +527,7 @@ extern const Name tau_spike; extern const Name tau_stc; extern const Name tau_syn; extern const Name tau_syn_ex; +extern const Name tau_syn_fast; extern const Name tau_syn_in; extern const Name tau_syn_slow; extern const Name tau_theta; diff --git a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py index 3b15223925..6a3fd7ee60 100644 --- a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py +++ b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py @@ -76,6 +76,9 @@ "gif_cond_exp_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "glif_cond": {"params": {"tau_syn": [1.0], "E_rev": [-85.0]}, "receptor_type": 1}, "glif_psc": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, + "glif_psc_double_alpha": { + "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, "receptor_type": 1, + }, "aeif_cond_alpha_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "aeif_cond_beta_multisynapse": { "params": {"E_rev": [0.0], "tau_rise": [1.0], "tau_decay": [1.0]}, diff --git a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py index 313b30afcc..188f28157f 100644 --- a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py +++ b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py @@ -66,6 +66,7 @@ "gif_psc_exp_multisynapse": {"receptor_type": 1}, "gif_cond_exp_multisynapse": {"receptor_type": 1}, "glif_psc": {"receptor_type": 1}, + "glif_psc_double_alpha": {"receptor_type": 1}, "iaf_cond_alpha_mc": {"receptor_type": 1}, "glif_cond": {"receptor_type": 1}, "ht_neuron": {"receptor_type": 1}, diff --git a/testsuite/pytests/test_glif_psc_double_alpha.py b/testsuite/pytests/test_glif_psc_double_alpha.py new file mode 100644 index 0000000000..3f4253e8dc --- /dev/null +++ b/testsuite/pytests/test_glif_psc_double_alpha.py @@ -0,0 +1,305 @@ +# -*- coding: utf-8 -*- +# +# test_glif_psc_double_alpha.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +import unittest +import nest + +try: + import scipy.stats + + HAVE_SCIPY = True +except ImportError: + HAVE_SCIPY = False + + +@unittest.skipIf(not HAVE_SCIPY, "SciPy package is not available") +class GLIFPSCTestCase(unittest.TestCase): + def setUp(self): + """ + Clean up and initialize NEST before each test. + """ + nest.ResetKernel() + nest.resolution = 0.01 + nest.rng_seed = 123456 + + def simulate_w_stim(self, model_params): + """ + Runs a one second simulation with different types of input stimulation. + Returns the time-steps, voltages, and collected spike times. + """ + + # Create glif model with double alpha synapses to simulate + # For this model, tau_syn_fast = tau_syn_slow = 2.0, and amp_slow = 1.0, + # so that this model behaves the same as glif_psc, when weight is halved + # of the original, for the connections that go into receptor_type == 1. + nrn = nest.Create("glif_psc_double_alpha", params=model_params) + nrn.set({"tau_syn_fast": [2.0], "tau_syn_slow": [2.0], "amp_slow": [1.0]}) + + # Create and connect inputs to glif model (weight is halved to account for double alpha synapse) + espikes = nest.Create( + "spike_generator", params={"spike_times": [10.0, 100.0, 400.0, 700.0], "spike_weights": [10.0] * 4} + ) + ispikes = nest.Create( + "spike_generator", params={"spike_times": [15.0, 99.0, 300.0, 800.0], "spike_weights": [-10.0] * 4} + ) + + cg = nest.Create( + "step_current_generator", + params={ + "amplitude_values": [ + 100.0, + ], + "amplitude_times": [ + 600.0, + ], + "start": 600.0, + "stop": 900.0, + }, + ) + pg = nest.Create("poisson_generator", params={"rate": 1000.0, "start": 200.0, "stop": 500.0}) + pn = nest.Create("parrot_neuron") + + nest.Connect(espikes, nrn, syn_spec={"receptor_type": 1}) + nest.Connect(ispikes, nrn, syn_spec={"receptor_type": 1}) + nest.Connect(cg, nrn, syn_spec={"weight": 3.0}) + nest.Connect(pg, pn) + nest.Connect(pn, nrn, syn_spec={"weight": 17.5, "receptor_type": 1}) # This is also halved + + # For recording spikes and voltage traces + sr = nest.Create("spike_recorder") + nest.Connect(nrn, sr) + + mm = nest.Create("multimeter", params={"record_from": ["V_m"]}) + nest.Connect(mm, nrn) + + nest.Simulate(1000.0) + + times = nest.GetStatus(mm, "events")[0]["times"] + V_m = nest.GetStatus(mm, "events")[0]["V_m"] + spikes = nest.GetStatus(sr, "events")[0]["times"] + + return times, V_m, spikes + + def ks_assert_spikes(self, spikes, reference_spikes): + """ + Runs a two-sided Kolmogorov-Smirnov statistic test on a set of spikes against a set of reference spikes. + """ + p_value_lim = 0.9 + d_lim = 0.2 + d, p_value = scipy.stats.ks_2samp(spikes, reference_spikes) + print(f"d={d}, p_value={p_value}") + self.assertGreater(p_value, p_value_lim) + self.assertLess(d, d_lim) + + def test_lif(self): + """ + Check LIF model + """ + lif_params = { + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_params) + spikes_expected = [ + 394.67, + 424.25, + 483.25, + 612.99, + 628.73, + 644.47, + 660.21, + 675.95, + 691.69, + 706.3, + 722.0, + 737.74, + 753.48, + 769.22, + 784.96, + 800.7, + 816.72, + 832.46, + 848.2, + 863.94, + 879.68, + 895.42, + ] + + self.ks_assert_spikes(spikes, spikes_expected) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r(self): + """ + Check LIF_R model + """ + lif_r_params = { + "spike_dependent_threshold": True, + "after_spike_currents": False, + "adapting_threshold": False, + "V_m": -78.85, + } + times, V_m, spikes = self.simulate_w_stim(lif_r_params) + expected_spikes = [ + 394.67, + 400.67, + 424.51, + 484.48, + 613.4, + 621.31, + 629.67, + 638.47, + 647.7, + 657.36, + 667.42, + 677.87, + 688.67, + 699.8, + 709.84, + 721.58, + 733.57, + 745.76, + 758.11, + 770.6, + 783.21, + 795.92, + 811.36, + 824.05, + 836.81, + 849.64, + 862.54, + 875.49, + 888.48, + ] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_asc(self): + """ + Check LIF_ASC model + """ + lif_asc_params = { + "spike_dependent_threshold": False, + "after_spike_currents": True, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_asc_params) + expected_spikes = [394.67, 484.81, 614.85, 648.28, 685.18, 725.06, 769.76, 821.69, 876.09] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r_asc(self): + """ + Check LIF_R_ASC model + """ + lif_r_asc_params = { + "spike_dependent_threshold": True, + "after_spike_currents": True, + "adapting_threshold": False, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_r_asc_params) + expected_spikes = [394.67, 485.07, 615.16, 649.52, 689.1, 734.18, 786.34, 845.85] + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_lif_r_asc_a(self): + """ + Check LIF_R_ASC_A model + """ + lif_r_asc_a_params = { + "spike_dependent_threshold": True, + "after_spike_currents": True, + "adapting_threshold": True, + "V_m": -78.85, + } + + times, V_m, spikes = self.simulate_w_stim(lif_r_asc_a_params) + expected_spikes = [395.44, 615.27, 653.77, 701.42, 768.22, 859.82] + + self.ks_assert_spikes(spikes, expected_spikes) + self.assertAlmostEqual(V_m[0], -78.85) + + def test_double_alpha_synapse(self): + """ + Check double alpha synapse + """ + neuron = nest.Create("glif_psc_double_alpha") + + # create 3 spikes as inputs + input1 = nest.Create("spike_generator", 1, {"spike_times": [10.0]}) + input2 = nest.Create("spike_generator", 1, {"spike_times": [210.0]}) + input3 = nest.Create("spike_generator", 1, {"spike_times": [410.0]}) + + # input1: double alpha at the same location with same amplitude. should double the amplitude + # input2: double alpha at the same location with different amplitude. net effect shoulld 1.5 of single + # input3: double alpha at different locations. The peak of the combined input shifts to ther right. + neuron_settings = { + "tau_syn_fast": [2.0, 2.0, 2.0], + "tau_syn_slow": [2.0, 2.0, 4.0], + "amp_slow": [1.0, 0.5, 0.5], + } + + neuron.set(neuron_settings) + # record the current with 0.1 ms resolution to capture the peak values. + multimeter = nest.Create("multimeter", params={"record_from": ["I_syn"], "interval": 0.1}) + + # set the delay to 1 ms. and amplitude is 1. + # the expected peak time of the fast components are, 13 ms, 213 ms, 413 ms. + nest.Connect(input1, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 1}) + nest.Connect(input2, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 2}) + nest.Connect(input3, neuron, syn_spec={"weight": 1.0, "delay": 1.0, "receptor_type": 3}) + nest.Connect(multimeter, neuron) + + # Do simulation. + nest.Simulate(500.0) + times = nest.GetStatus(multimeter, "events")[0]["times"] + I_syn = nest.GetStatus(multimeter, "events")[0]["I_syn"] + + # Check the results. + # the peak timing should be 13 ms (index 129), 213 ms (index 2129), + # and slightly above 413 ms (413.4 ms, or index 4133). + + # for the first two we can estimate the exact peak values (2 and 1.5 pA) + self.assertAlmostEqual(I_syn[129], 2.0, delta=0.0001) + self.assertAlmostEqual(I_syn[2129], 1.5, delta=0.0001) + # for the last one, check if the the new peak is larger than peak of the fast component + self.assertGreater(I_syn[4133], I_syn[4129]) # peak is shifted to the right. + + +def suite(): + return unittest.makeSuite(GLIFPSCTestCase, "test") + + +def run(): + runner = unittest.TextTestRunner(verbosity=2) + runner.run(suite()) + + +if __name__ == "__main__": + run() From 571dbc3b4d1030373f6bad766f1b6f5c7f19d7f2 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Fri, 25 Aug 2023 18:35:13 -0700 Subject: [PATCH 25/45] Making double-alpha description clearer --- models/glif_psc.h | 2 +- models/glif_psc_double_alpha.h | 26 +++++++++++++++----------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/models/glif_psc.h b/models/glif_psc.h index cb11214933..e85c48f182 100644 --- a/models/glif_psc.h +++ b/models/glif_psc.h @@ -165,7 +165,7 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn double vector Rise time constants of the synaptic +tau_syn double vector Time constants of the synaptic alpha function in ms E_rev double vector Reversal potential in mV spike_dependent_threshold bool flag whether the neuron has diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index cdc7449529..232929a079 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -44,16 +44,20 @@ Description +++++++++++ ``glif_psc_double_alpha`` provides five generalized leaky integrate-and-fire -(GLIF) models [1]_ with alpha-function shaped synaptic currents. +(GLIF) models [1]_ with double alpha-function shaped synaptic currents. Incoming spike events induce a postsynaptic change of current modeled -by sum of two alpha functions [2]_. The alpha function is normalized such that an event -of weight 1.0 results in a peak current of the fast component of the alpha function to -be 1 pA at :math:`t = tau_syn_fast`. The peak current of the slow component of the alpha -is given as amp_slow pA, at :math:`t = tau_syn_slow`. By default, glif_psc_double_alpha -has a single synapse that is accessible through receptor_port 1. An arbitrary number of -synapses with different time constants can be configured by setting the desired time -constants as tau_syn_fast array. The resulting synapses are addressed through -receptor_port 1, 2, 3, .... +by the sum of two alpha functions (fast and slow components) for each receptor [2]_. +This function is normalized such that an event of weight 1.0 results in a peak current +of the fast component of the alpha function to be 1 pA at :math:`t = tau_syn_fast`. +The relative peak current of the slow component is given as amp_slow, at +:math:`t = tau_syn_slow`. Namely, +:math:`I_{syn} = alpha_function(tau_syn=tau_syn_fast) + amp_slow * +alpha_function(tau_syn=tau_syn_slow)`. Therefore if amp_slow is not 0, the peak current +of the total synaptic current is larger than the specified weight. By default, +glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. +An arbitrary number of synapses with different time constants and amp_slow can be +configured by setting the desired parameters of tau_syn_fast, tau_syn_slow, and amp_slow +arrays. The resulting synapses are addressed through receptor_port 1, 2, 3, .... The five GLIF models are: @@ -168,9 +172,9 @@ th_voltage_decay double Voltage-induced threshold time voltage-dependent component of the threshold in 1/ms (bv in Equation (4) in [1]_) -tau_syn_fast double vector Rise time constants of the faster +tau_syn_fast double vector Time constants of the faster synaptic alpha function in ms -tau_syn_slow double vector Rise time constants of the slower +tau_syn_slow double vector Time constants of the slower synaptic alpha function in ms amp_slow double vector Relative amplitude of the slower synaptic alpha function From 2edb4482b750ff27a81204c99bbca09608177860 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Tue, 29 Aug 2023 18:45:44 -0700 Subject: [PATCH 26/45] Making math format LaTeX compatible --- models/glif_psc_double_alpha.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 232929a079..f818ac80d3 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -48,16 +48,18 @@ Description Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current -of the fast component of the alpha function to be 1 pA at :math:`t = tau_syn_fast`. +of the fast component of the alpha function to be 1 pA at +:math:`t = \tau_\text{syn\_fast}`. The relative peak current of the slow component is given as amp_slow, at -:math:`t = tau_syn_slow`. Namely, -:math:`I_{syn} = alpha_function(tau_syn=tau_syn_fast) + amp_slow * -alpha_function(tau_syn=tau_syn_slow)`. Therefore if amp_slow is not 0, the peak current -of the total synaptic current is larger than the specified weight. By default, -glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. -An arbitrary number of synapses with different time constants and amp_slow can be -configured by setting the desired parameters of tau_syn_fast, tau_syn_slow, and amp_slow -arrays. The resulting synapses are addressed through receptor_port 1, 2, 3, .... +:math:`t = \tau_\text{syn\_slow}`. Namely, +:math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + +\text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. +Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger +than the specified weight. By default, glif_psc_double_alpha has a single synapse that +is accessible through receptor_port 1. An arbitrary number of synapses with different +time constants and amp_slow can be configured by setting the desired parameters of +tau_syn_fast, tau_syn_slow, and amp_slow arrays. The resulting synapses are addressed +through receptor_port 1, 2, 3, .... The five GLIF models are: From 45b872d4bcb9f4687e87c081ce2f8dd24159c5cf Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Wed, 30 Aug 2023 13:40:51 -0700 Subject: [PATCH 27/45] black applied to a modified file --- .../pytests/sli2py_neurons/test_neurons_handle_multiplicity.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py index 6a3fd7ee60..1812745e2c 100644 --- a/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py +++ b/testsuite/pytests/sli2py_neurons/test_neurons_handle_multiplicity.py @@ -77,7 +77,8 @@ "glif_cond": {"params": {"tau_syn": [1.0], "E_rev": [-85.0]}, "receptor_type": 1}, "glif_psc": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "glif_psc_double_alpha": { - "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, "receptor_type": 1, + "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, + "receptor_type": 1, }, "aeif_cond_alpha_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "aeif_cond_beta_multisynapse": { From e44264a7671d4f74b2320b119240f8d9e4c86be3 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:30 -0700 Subject: [PATCH 28/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index f818ac80d3..202e650078 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -49,7 +49,7 @@ Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at -:math:`t = \tau_\text{syn\_fast}`. +:math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn\_slow}`. Namely, :math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + From 089db1149a3150199a4d8df0e6f9b16377213af8 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:45 -0700 Subject: [PATCH 29/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 202e650078..84c15b322f 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -51,7 +51,7 @@ This function is normalized such that an event of weight 1.0 results in a peak c of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at -:math:`t = \tau_\text{syn\_slow}`. Namely, +:math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + \text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger From 40537ee4ec66df1f590043965e62ebcfb02138f1 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:28:55 -0700 Subject: [PATCH 30/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 84c15b322f..93cb4e5766 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -52,7 +52,7 @@ of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn_fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, -:math:`I_\text{syn} = \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_fast}) + +:math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + \text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that From 4309902685b0b8cc00514392e700e9bea041dd50 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 30 Aug 2023 13:29:06 -0700 Subject: [PATCH 31/45] Update models/glif_psc_double_alpha.h Co-authored-by: clinssen --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 93cb4e5766..749105cfe9 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -53,7 +53,7 @@ of the fast component of the alpha function to be 1 pA at The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + -\text{amp_slow} * \text{alpha\_function}(\tau_\text{syn} = \tau_\text{syn\_slow})`. +\text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different From db0be04980aeaa973f4a20882c04e39b0c453447 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Wed, 30 Aug 2023 13:45:30 -0700 Subject: [PATCH 32/45] clang-format -i applied to committed files --- models/glif_psc_double_alpha.cpp | 5 +++-- models/glif_psc_double_alpha.h | 10 +++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/models/glif_psc_double_alpha.cpp b/models/glif_psc_double_alpha.cpp index 562154488e..6245501a88 100644 --- a/models/glif_psc_double_alpha.cpp +++ b/models/glif_psc_double_alpha.cpp @@ -80,7 +80,7 @@ nest::glif_psc_double_alpha::Parameters_::Parameters_() , asc_decay_( std::vector< double > { 0.003, 0.1 } ) // in 1/ms , asc_amps_( std::vector< double > { -9.18, -198.94 } ) // in pA , asc_r_( std::vector< double >( 2, 1.0 ) ) - , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms + , tau_syn_fast_( std::vector< double >( 1, 2.0 ) ) // in ms , tau_syn_slow_( std::vector< double >( 1, 6.0 ) ) // in ms , amp_slow_( std::vector< double >( 1, 0.3 ) ) // amplitude relative to the fast component, unitless , has_connections_( false ) @@ -493,7 +493,8 @@ nest::glif_psc_double_alpha::pre_run_hook() // these are determined according to a numeric stability criterion // input time parameter shall be in ms, capacity in pF - std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); + std::tie( V_.P31_fast_[ i ], V_.P32_fast_[ i ] ) = + IAFPropagatorAlpha( P_.tau_syn_fast_[ i ], Tau_, P_.C_m_ ).evaluate( h ); // Slow component V_.P11_slow_[ i ] = V_.P22_slow_[ i ] = std::exp( -h / P_.tau_syn_slow_[ i ] ); diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 749105cfe9..2be27965e7 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -270,7 +270,7 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > asc_decay_; //!< predefined time scale in 1/ms std::vector< double > asc_amps_; //!< in pA std::vector< double > asc_r_; //!< coefficient - std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms + std::vector< double > tau_syn_fast_; //!< synaptic port time constants in ms std::vector< double > tau_syn_slow_; //!< synaptic port time constants in ms std::vector< double > amp_slow_; //!< synaptic port time constants in ms @@ -307,8 +307,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > ASCurrents_; //!< after-spike currents in pA double ASCurrents_sum_; //!< in pA int refractory_steps_; //!< Number of refractory steps remaining - std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA - std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA + std::vector< double > y1_fast_; //!< synapse current evolution state 1 in pA + std::vector< double > y2_fast_; //!< synapse current evolution state 2 in pA std::vector< double > y1_slow_; //!< synapse current evolution state 1 in pA std::vector< double > y2_slow_; //!< synapse current evolution state 2 in pA @@ -347,8 +347,8 @@ class glif_psc_double_alpha : public ArchivingNode std::vector< double > P11_fast_; //!< synaptic current evolution parameter std::vector< double > P21_fast_; //!< synaptic current evolution parameter std::vector< double > P22_fast_; //!< synaptic current evolution parameter - double P30_; //!< membrane current/voltage evolution parameter - double P33_; //!< membrane voltage evolution parameter + double P30_; //!< membrane current/voltage evolution parameter + double P33_; //!< membrane voltage evolution parameter std::vector< double > P31_fast_; //!< synaptic/membrane current evolution parameter std::vector< double > P32_fast_; //!< synaptic/membrane current evolution parameter From 4c354f2670557c41e8cb2ee1fd45d8901423edf7 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Thu, 31 Aug 2023 10:34:49 -0700 Subject: [PATCH 33/45] Adding double_alpha reference to glif_psc.h --- models/glif_psc.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/models/glif_psc.h b/models/glif_psc.h index e85c48f182..d133c9f6b4 100644 --- a/models/glif_psc.h +++ b/models/glif_psc.h @@ -192,7 +192,8 @@ References See also ++++++++ -gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp +gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp, +glif_psc_double_alpha EndUserDocs */ From 983998435eb51aa055600a5b4a65981394764fdb Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:28:54 -0700 Subject: [PATCH 34/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2be27965e7..a252d7e3c8 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -37,7 +37,7 @@ Short description +++++++++++++++++ -Current-based generalized leaky integrate-and-fire (GLIF) models +Current-based generalized leaky integrate-and-fire (GLIF) models with double alpha-function (from the Allen Institute) Description From c553f21c277a3689592c36ef2974706f03b5a0b4 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:30:13 -0700 Subject: [PATCH 35/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index a252d7e3c8..2298289b57 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -49,7 +49,7 @@ Incoming spike events induce a postsynaptic change of current modeled by the sum of two alpha functions (fast and slow components) for each receptor [2]_. This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at -:math:`t = \tau_\text{syn_fast}`. +:math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as amp_slow, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + From ad016b3e77ab16b7d9432d5853940d10de2ff97c Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:30:29 -0700 Subject: [PATCH 36/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 2298289b57..4ca0b54dee 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -50,7 +50,7 @@ by the sum of two alpha functions (fast and slow components) for each receptor [ This function is normalized such that an event of weight 1.0 results in a peak current of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. -The relative peak current of the slow component is given as amp_slow, at +The relative peak current of the slow component is given as ``amp_slow``, at :math:`t = \tau_\text{syn_slow}`. Namely, :math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + \text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. From 75b3ccbdb3195bda137bc9d209d4212bc9926249 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:00 -0700 Subject: [PATCH 37/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 4ca0b54dee..80c21a2cbc 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -52,8 +52,11 @@ of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as ``amp_slow``, at :math:`t = \tau_\text{syn_slow}`. Namely, -:math:`I_\text{syn} = \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_fast}) + -\text{amp_slow} * \text{alpha_function}(\tau_\text{syn} = \tau_\text{syn_slow})`. + +.. math:: + + I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). + Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, glif_psc_double_alpha has a single synapse that is accessible through receptor_port 1. An arbitrary number of synapses with different From 7ed402dca8789a7d9abffa58025823381f0e2eb0 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:15 -0700 Subject: [PATCH 38/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 80c21a2cbc..920a2bd326 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -57,12 +57,12 @@ The relative peak current of the slow component is given as ``amp_slow``, at I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). -Therefore if amp_slow is not 0, the peak current of the total synaptic current is larger -than the specified weight. By default, glif_psc_double_alpha has a single synapse that -is accessible through receptor_port 1. An arbitrary number of synapses with different -time constants and amp_slow can be configured by setting the desired parameters of -tau_syn_fast, tau_syn_slow, and amp_slow arrays. The resulting synapses are addressed -through receptor_port 1, 2, 3, .... +Therefore if ``amp_slow`` is not 0, the peak current of the total synaptic current is larger +than the specified weight. By default, ``glif_psc_double_alpha`` has a single synapse that +is accessible through ``receptor_port`` 1. An arbitrary number of synapses with different +time constants and ``amp_slow`` can be configured by setting the desired parameters of +``tau_syn_fast``, ``tau_syn_slow``, and ``amp_slow`` arrays. The resulting synapses are addressed +through ``receptor_port`` 1, 2, 3, .... The five GLIF models are: From e9232a9fb94e9723c6e3b2b69241ebeea618d62b Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Thu, 31 Aug 2023 10:31:30 -0700 Subject: [PATCH 39/45] Update models/glif_psc_double_alpha.h Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index 920a2bd326..b4bcbde9f3 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -105,7 +105,7 @@ F, s, A) to NEST used units (i.e., mV, nS (1/GOhm), pF, ms, pA) and values being rounded to appropriate digits for simplification. For models with spike dependent threshold (i.e., GLIF2, GLIF4 and GLIF5), -parameter setting of voltage_reset_fraction and voltage_reset_add may lead to the +parameter setting of ``voltage_reset_fraction`` and ``voltage_reset_add`` may lead to the situation that voltage is bigger than threshold after reset. In this case, the neuron will continue to spike until the end of the simulation regardless the stimulated inputs. We recommend the setting of the parameters of these three models to follow the From 4f7ee9b0295c8a42534349905193b72863b23753 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Thu, 31 Aug 2023 10:38:32 -0700 Subject: [PATCH 40/45] One more change for math notation --- models/glif_psc_double_alpha.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index b4bcbde9f3..d0ac80de1f 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -51,7 +51,7 @@ This function is normalized such that an event of weight 1.0 results in a peak c of the fast component of the alpha function to be 1 pA at :math:`t = \tau_\text{syn, fast}`. The relative peak current of the slow component is given as ``amp_slow``, at -:math:`t = \tau_\text{syn_slow}`. Namely, +:math:`t = \tau_\text{syn, slow}`. Namely, .. math:: From 6d6e9bd147e7ad71c5225700991dfd6b13c375b0 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Fri, 1 Sep 2023 11:27:11 -0700 Subject: [PATCH 41/45] Adding example, and reapplying clang-format --- doc/htmldoc/examples/index.rst | 2 + models/glif_psc_double_alpha.h | 3 +- .../examples/glif_psc_double_alpha_neuron.py | 214 ++++++++++++++++++ 3 files changed, 218 insertions(+), 1 deletion(-) create mode 100644 pynest/examples/glif_psc_double_alpha_neuron.py diff --git a/doc/htmldoc/examples/index.rst b/doc/htmldoc/examples/index.rst index 549e6024cc..d37f8b4102 100644 --- a/doc/htmldoc/examples/index.rst +++ b/doc/htmldoc/examples/index.rst @@ -76,6 +76,7 @@ PyNEST examples * :doc:`../auto_examples/glif_cond_neuron` * :doc:`../auto_examples/glif_psc_neuron` + * :doc:`../auto_examples/glif_psc_double_alpha_neuron` @@ -260,6 +261,7 @@ PyNEST examples ../auto_examples/mc_neuron ../auto_examples/glif_cond_neuron ../auto_examples/glif_psc_neuron + ../auto_examples/glif_psc_double_alpha_neuron ../auto_examples/precise_spiking ../auto_examples/CampbellSiegert ../auto_examples/vinit_example diff --git a/models/glif_psc_double_alpha.h b/models/glif_psc_double_alpha.h index d0ac80de1f..b0589d03b2 100644 --- a/models/glif_psc_double_alpha.h +++ b/models/glif_psc_double_alpha.h @@ -55,7 +55,8 @@ The relative peak current of the slow component is given as ``amp_slow``, at .. math:: - I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). + I_\text{syn} = \text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, fast} \right) + \text{amp_slow} \cdot +\text{alpha_function} \left( \tau_\text{syn} = \tau_\text{syn, slow} \right). Therefore if ``amp_slow`` is not 0, the peak current of the total synaptic current is larger than the specified weight. By default, ``glif_psc_double_alpha`` has a single synapse that diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py new file mode 100644 index 0000000000..e448382bea --- /dev/null +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -0,0 +1,214 @@ +# -*- coding: utf-8 -*- +# +# glif_psc_double_alpha_neuron.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +# %% +""" +Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha +synaptic function example +------------------------------------------------------------------------ + +Simple example of how to use the ``glif_psc_double_alpha`` neuron model that illustrates +difference from the ``glif_psc`` neuron model. + +The behavior of the ``glif_psc_double_alpha`` neuron model the same as the ``glif_psc`` +neuron model, except that the synaptic currents are modeled as a double alpha function. +Therefore, in this example, we only compare the difference in the synaptic currents +between the two models. compared to the single alpha function, the double alpha function +has much more control over the shape of the tail of the synaptic current. + +Simple synaptic inputs are applied to the neuron and the resulting voltage traces, +current traces are shown for the two models. + +""" + +############################################################################## +# First, we import all necessary modules to simulate, analyze and plot this +# example. + +import nest +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec + +############################################################################## +# We initialize NEST and set the simulation resolution. + +nest.ResetKernel() +resolution = 0.05 +nest.resolution = resolution + +############################################################################## +# We also pre-define the synapse time constant arrays. +# In contrast to ``glif_psc`` models, ``glif_psc_double_alpha`` models have +# two components of synaptic currents, one for the fast component and the other +# for the slow component. The relative amplitude also need to be set, so there +# are three parameters to define per port. For this example, we keep the +# ``tau_syn_fast`` to 2 ms for simplicity, and vary the ``tau_syn_slow`` and +# ``amp_slow`` to illustrate how the parameters affect the shape of the synaptic +# currents. + +tau_syn_glif_psc = [2.0, 2.0, 2.0] # value for the ``glif_psc`` model + +tau_syn_fast = [2.0, 2.0, 2.0] # common between 'timing' and 'amp' manipulations + +# for the slow component timing manipuation +tau_syn_slow_timing = [4.0, 6.0, 8.0] +amp_slow_timing = [0.5, 0.5, 0.5] + +# for the slow component amplitude manipulation +tau_syn_slow_amp = [6.0, 6.0, 6.0] +amp_slow_amp = [0.2, 0.5, 0.8] + +############################################################################### +# Now we create three neurons: ``glif_psc``, ``glif_psc_double_alpha_timing``, +# and ``glif_psc_double_alpha_amp``. The parameters for the ``glif_psc`` neuron +# are set as default. The parameters for the ``glif_psc_double_alpha_timing`` +# neuron are set to have the same ``tau_syn_fast`` as the ``glif_psc`` neuron, +# and the ``tau_syn_slow`` and ``amp_slow`` are set to the values defined above +# for the timing manipulation. + +n_glif_psc = nest.Create( + "glif_psc", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn": tau_syn_glif_psc, + }, +) + +n_glif_psc_double_alpha_timing = nest.Create( + "glif_psc_double_alpha", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn_fast": tau_syn_fast, + "tau_syn_slow": tau_syn_slow_timing, + "amp_slow": amp_slow_timing, + }, +) + +n_glif_psc_double_alpha_amp = nest.Create( + "glif_psc_double_alpha", + params={ + "spike_dependent_threshold": False, + "after_spike_currents": False, + "adapting_threshold": False, + "tau_syn_fast": tau_syn_fast, + "tau_syn_slow": tau_syn_slow_amp, + "amp_slow": amp_slow_amp, + }, +) + +neurons = n_glif_psc + n_glif_psc_double_alpha_timing + n_glif_psc_double_alpha_amp + +############################################################################### +# For the stimulation input to the glif_psc neurons, we create three excitation +# spike generators, each one with a single spike. + +espike1 = nest.Create("spike_generator", params={"spike_times": [10.0], "spike_weights": [20.0]}) +espike2 = nest.Create("spike_generator", params={"spike_times": [110.0], "spike_weights": [20.0]}) +espike3 = nest.Create("spike_generator", params={"spike_times": [210.0], "spike_weights": [20.0]}) + +############################################################################### +# The generators are then connected to the neurons. Specification of +# the ``receptor_type`` uniquely defines the target receptor. +# We connect each of the spikes generator to a different receptor that have different +# parameters. + +nest.Connect(espike1, neurons, syn_spec={"delay": resolution, "receptor_type": 1}) +nest.Connect(espike2, neurons, syn_spec={"delay": resolution, "receptor_type": 2}) +nest.Connect(espike3, neurons, syn_spec={"delay": resolution, "receptor_type": 3}) + +############################################################################### +# A ``multimeter`` is created and connected to the neurons. The parameters +# specified for the multimeter include the list of quantities that should be +# recorded and the time interval at which quantities are measured. + +mm = nest.Create( + "multimeter", + params={ + "interval": resolution, + "record_from": ["V_m", "I_syn"], + }, +) +nest.Connect(mm, neurons) + +############################################################################### +# Run the simulation for 300 ms and retrieve recorded data from +# the multimeter and spike recorder. + +nest.Simulate(300.0) +data = mm.events + +############################################################################### +# We plot the time traces of the synaptic current and the membrane potential. +# Each input current is annotated with the corresponding parameter value of the +# receptor. The blue line is the synaptic current of the ``glif_psc`` neuron, and +# the red line is the synaptic current of the ``glif_psc_double_alpha`` neuron. + + +# defining the figure property for each parameter variation type, +variation_types = ["timing", "amp"] +annotate_variable = ["tau_syn_slow", "amp_slow"] +annotate_values = [tau_syn_slow_timing, amp_slow_amp] +fig_titles = [ + "Variation of tau_syn_slow: tau_syn_fast = 2.0, amp_slow = 0.5", + "Variation of amp_slow: tau_syn_fast = 2.0, tau_syn_slow = 6.0", +] + + +senders = data["senders"] +t = data["times"][senders == 1] + +for i, variation_type in enumerate(variation_types): + plt.figure(variation_type, figsize=(10, 5)) + gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1]) + data_types = ["I_syn", "V_m"] + data_types_names = ["Synaptic current (pA)", "Membrane potential (mV)"] + for j, data_type in enumerate(data_types): + d = data[data_type] + ax = plt.subplot(gs[j]) + ax.plot(t, d[senders == 1], "b", label="glif_psc (tau_syn=2.0)") + ax.plot(t, d[senders == 2 + i], "r", label="glif_psc_double_alpha") + if j == 0: + # place legend outside the plot + ax.legend(bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0) + else: + ax.set_xlabel("time (ms)") + + ax.set_ylabel(data_types_names[j]) + + # now let's annotate each of the input with the corresponding parameter. + spike_timings = [10.0, 110.0, 210.0] + ax = plt.subplot(gs[0]) + for j, spike_timing in enumerate(spike_timings): + ax.annotate( + f"{annotate_variable[i]}={annotate_values[i][j]}", + xy=(spike_timing + 10, 20), + xytext=(spike_timing + 10, 25), + arrowprops=dict(arrowstyle="->", connectionstyle="arc3"), + ) + plt.title(fig_titles[i]) + plt.tight_layout() + + +plt.show() From 2152cbd42ba62216359aff9ab36187825ddbef45 Mon Sep 17 00:00:00 2001 From: Dennis Terhorst Date: Mon, 4 Sep 2023 11:26:07 +0200 Subject: [PATCH 42/45] Apply suggestions from code review Co-authored-by: Nicolai Haug <39106781+nicolossus@users.noreply.github.com> --- pynest/examples/glif_psc_double_alpha_neuron.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py index e448382bea..50895ccb70 100644 --- a/pynest/examples/glif_psc_double_alpha_neuron.py +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -19,7 +19,6 @@ # You should have received a copy of the GNU General Public License # along with NEST. If not, see . -# %% """ Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha synaptic function example @@ -28,15 +27,14 @@ Simple example of how to use the ``glif_psc_double_alpha`` neuron model that illustrates difference from the ``glif_psc`` neuron model. -The behavior of the ``glif_psc_double_alpha`` neuron model the same as the ``glif_psc`` +The behavior of the ``glif_psc_double_alpha`` neuron model is the same as the ``glif_psc`` neuron model, except that the synaptic currents are modeled as a double alpha function. Therefore, in this example, we only compare the difference in the synaptic currents -between the two models. compared to the single alpha function, the double alpha function +between the two models. Compared to the single alpha function, the double alpha function has much more control over the shape of the tail of the synaptic current. -Simple synaptic inputs are applied to the neuron and the resulting voltage traces, +Simple synaptic inputs are applied to the neuron and the resulting voltage and current traces are shown for the two models. - """ ############################################################################## @@ -58,8 +56,8 @@ # We also pre-define the synapse time constant arrays. # In contrast to ``glif_psc`` models, ``glif_psc_double_alpha`` models have # two components of synaptic currents, one for the fast component and the other -# for the slow component. The relative amplitude also need to be set, so there -# are three parameters to define per port. For this example, we keep the +# for the slow component. The relative amplitude also needs to be set, so there +# are three parameters to define per receptor port. For this example, we keep the # ``tau_syn_fast`` to 2 ms for simplicity, and vary the ``tau_syn_slow`` and # ``amp_slow`` to illustrate how the parameters affect the shape of the synaptic # currents. @@ -121,7 +119,7 @@ neurons = n_glif_psc + n_glif_psc_double_alpha_timing + n_glif_psc_double_alpha_amp ############################################################################### -# For the stimulation input to the glif_psc neurons, we create three excitation +# For the stimulation input to the ``glif_psc`` neurons, we create three excitation # spike generators, each one with a single spike. espike1 = nest.Create("spike_generator", params={"spike_times": [10.0], "spike_weights": [20.0]}) From 8500f3edd0e17d33c80f3da6ad5628069d6f49d7 Mon Sep 17 00:00:00 2001 From: "shinya.ito" Date: Thu, 7 Sep 2023 12:42:51 -0700 Subject: [PATCH 43/45] Making the new test code compatible with glif_psc_double_alpha. --- testsuite/pytests/sli2py_regressions/test_issue_77.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/testsuite/pytests/sli2py_regressions/test_issue_77.py b/testsuite/pytests/sli2py_regressions/test_issue_77.py index fb93153584..020c827590 100644 --- a/testsuite/pytests/sli2py_regressions/test_issue_77.py +++ b/testsuite/pytests/sli2py_regressions/test_issue_77.py @@ -79,6 +79,10 @@ "gif_psc_exp_multisynapse": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, "glif_cond": {"params": {"tau_syn": [0.2], "E_rev": [0.0]}, "receptor_type": 1}, "glif_psc": {"params": {"tau_syn": [1.0]}, "receptor_type": 1}, + "glif_psc_double_alpha": { + "params": {"tau_syn_fast": [1.0], "tau_syn_slow": [2.0], "amp_slow": [0.5]}, + "receptor_type": 1, + }, "ht_neuron": {"receptor_type": 1}, "pp_cond_exp_mc_urbanczik": {"receptor_type": 1}, } From 5df9403e919d796153bf32a89bb9ba70c3e825e8 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 13 Sep 2023 17:43:04 -0700 Subject: [PATCH 44/45] Update pynest/examples/glif_psc_double_alpha_neuron.py Co-authored-by: jessica-mitchell --- pynest/examples/glif_psc_double_alpha_neuron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py index 50895ccb70..6a7a700f1a 100644 --- a/pynest/examples/glif_psc_double_alpha_neuron.py +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -20,7 +20,7 @@ # along with NEST. If not, see . """ -Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha +Current-based generalized leaky integrate and fire (GLIF) neuron with double alpha \ synaptic function example ------------------------------------------------------------------------ From 1e9e65efc20bd470b449b8f4ad69db706e911e80 Mon Sep 17 00:00:00 2001 From: Shinya Ito Date: Wed, 13 Sep 2023 17:43:13 -0700 Subject: [PATCH 45/45] Update pynest/examples/glif_psc_double_alpha_neuron.py Co-authored-by: jessica-mitchell --- pynest/examples/glif_psc_double_alpha_neuron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynest/examples/glif_psc_double_alpha_neuron.py b/pynest/examples/glif_psc_double_alpha_neuron.py index 6a7a700f1a..44787f01e2 100644 --- a/pynest/examples/glif_psc_double_alpha_neuron.py +++ b/pynest/examples/glif_psc_double_alpha_neuron.py @@ -25,7 +25,7 @@ ------------------------------------------------------------------------ Simple example of how to use the ``glif_psc_double_alpha`` neuron model that illustrates -difference from the ``glif_psc`` neuron model. +differences from the ``glif_psc`` neuron model. The behavior of the ``glif_psc_double_alpha`` neuron model is the same as the ``glif_psc`` neuron model, except that the synaptic currents are modeled as a double alpha function.