gif_pop_psc_exp.h

Go to the documentation of this file.

/*
 *  gif_pop_psc_exp.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 <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef PP_POP_PSC_BETA_H
#define PP_POP_PSC_BETA_H
 
 
// Includes from nestkernel:
#include "nest.h"
#include "node.h"
#include "random_generators.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
#ifdef HAVE_GSL
 
namespace nest
{
 
class Network;
 
/* BeginUserDocs: neuron, integrate-and-fire, current-based, adaptation, stochastic
 
Short description
+++++++++++++++++
 
Population of generalized integrate-and-fire neurons (GIF) with exponential
postsynaptic currents and adaptation (from the Gerstner lab)
 
 
Description
+++++++++++
 
This model simulates a population of spike-response model neurons with
multi-timescale adaptation and exponential postsynaptic currents, as
described by Schwalger et al. (2017) :footcite:p:`Schwalger2017`.
 
The single neuron model is defined by the hazard function
 
.. math::
 
 h(t) = \lambda_0  \exp\frac{V_m(t) - E_{\text{sfa}}(t)}{\Delta_V}
 
After each spike, the membrane potential :math:`V_m` is reset to
:math:`V_{\text{reset}}`. Spike frequency
adaptation is implemented by a set of exponentially decaying traces, the
sum of which is :math:`E_{\text{sfa}}`. Upon a spike, each of the adaptation traces is
incremented by the respective :math:`q_{\text{sfa}}` and decays with the respective time constant
:math:`\tau_{\text{sfa}}`.
 
The corresponding single neuron model is available in NEST as ``gif_psc_exp``.
The default parameters, although some are named slightly different, are not
matched in both models for historical reasons. See below for the parameter
translation.
 
Connecting two population models corresponds to full connectivity of every
neuron in each population. An approximation of random connectivity can be
implemented by connecting populations using a ``bernoulli_synapse``.
 
 
Parameters
++++++++++
 
The following parameters can be set in the status Dictionary.
 
 
=========== ============= =====================================================
 V_reset    mV            Membrane potential is reset to this value after
                          a spike
 V_T_star   mV            Threshold level of the membrane potential
 E_L        mV            Resting potential
 Delta_V    mV            Noise level of escape rate
 C_m        pF            Capacitance of the membrane
 tau_m      ms            Membrane time constant
 t_ref      ms            Duration of refractory period
 I_e        pA            Constant input current
 N          integer       Number of neurons in the population
 len_kernel integer       Refractory effects are accounted for up to len_kernel
                          time steps
 lambda_0   1/s           Firing rate at threshold
 tau_syn_ex ms            Time constant for excitatory synaptic currents
 tau_syn_in ms            Time constant for inhibitory synaptic currents
 tau_sfa    list of ms    vector Adaptation time constants
 q_sfa      list of ms    Adaptation kernel amplitudes
 BinoRand   boolean       If True, binomial random numbers are used, otherwise
                          we use Poisson distributed spike counts
=========== ============= =====================================================
 
 
=============== ============  =============================
**Parameter translation to gif_psc_exp**
-----------------------------------------------------------
gif_pop_psc_exp  gif_psc_exp  relation
tau_m            g_L          tau_m = C_m / g_L
N                ---          use N gif_psc_exp neurons
=============== ============  =============================
 
 
References
++++++++++
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
See also
++++++++
 
gif_psc_exp, bernoulli_synapse
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: gif_pop_psc_exp
 
EndUserDocs */
 
 
void register_gif_pop_psc_exp( const std::string& name );
 
class gif_pop_psc_exp : public Node
{
 
public:
  gif_pop_psc_exp();
  gif_pop_psc_exp( const gif_pop_psc_exp& );
 
  using Node::handle;
  using Node::handles_test_event;
 
  size_t send_test_event( Node&, size_t, synindex, bool ) override;
 
  void handle( SpikeEvent& ) override;
  void handle( CurrentEvent& ) override;
  void handle( DataLoggingRequest& ) override;
 
  size_t handles_test_event( SpikeEvent&, size_t ) override;
  size_t handles_test_event( CurrentEvent&, size_t ) override;
  size_t handles_test_event( DataLoggingRequest&, size_t ) override;
 
  void get_status( Dictionary& ) const override;
  void set_status( const Dictionary& ) override;
 
private:
  void init_buffers_() override;
  void pre_run_hook() override;
 
  void update( Time const&, const long, const long ) override;
 
  double escrate( const double );
  long draw_poisson( const double n_expect_ );
  long draw_binomial( const double n_expect_ );
  double adaptation_kernel( const int k );
  int get_history_size();
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< gif_pop_psc_exp >;
  friend class UniversalDataLogger< gif_pop_psc_exp >;
 
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    long N_;
 
    double tau_m_;
 
    double c_m_;
 
    double t_ref_;
 
    double lambda_0_;
 
    double Delta_V_;
 
    long len_kernel_;
 
    double I_e_;
 
    double V_reset_;
 
    double V_T_star_;
 
    double E_L_;
 
    double tau_syn_ex_;
    double tau_syn_in_;
 
    std::vector< double > tau_sfa_;
 
    std::vector< double > q_sfa_;
 
    bool BinoRand_;
 
    Parameters_();                              
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
  };
 
  // ----------------------------------------------------------------
 
  struct State_
  {
    double y0_;         // DC input current
    double I_syn_ex_;   // synaptic current
    double I_syn_in_;   // synaptic current
    double V_m_;        // membrane potential
    double n_expect_;   // expected spike number
    double theta_hat_;  // adapting threshold for non-refractory neurons
    long n_spikes_;     // number of spikes
 
    // internal switch signaling that state vectors are initialized
    bool initialized_;
 
    State_();  
 
    void get( Dictionary&, const Parameters_& ) const;
    void set( const Dictionary&, const Parameters_&, Node* );
  };
 
  // ----------------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( gif_pop_psc_exp& );
    Buffers_( const Buffers_&, gif_pop_psc_exp& );
 
    RingBuffer ex_spikes_;
    RingBuffer in_spikes_;
    RingBuffer currents_;
 
    UniversalDataLogger< gif_pop_psc_exp > logger_;
  };
 
  // ----------------------------------------------------------------
 
  struct Variables_
  {
 
    double R_;       // membrane resistance
    double P20_;     // membrane integration constant
    double P22_;     // membrane integration constant
    double P11_ex_;  // synaptic integration constant
    double P11_in_;  // synaptic integration constant
    int k_ref_;      // length of refractory period in time steps
 
    std::vector< double > Q30_;        // QR adaptation integration constant
    std::vector< double > Q30K_;       // QR adaptation integration constant
    std::vector< double > theta_;      // adaptation kernel
    std::vector< double > theta_tld_;  // QR adaptation kernel
 
    double h_;  // simulation time step in ms
    double min_double_;
 
    RngPtr rng_;  // random number generator of own thread
 
    poisson_distribution poisson_dist_;  
    binomial_distribution bino_dist_;    
 
    double x_;                      // internal variable of population dynamics
    double z_;                      // internal variable of population dynamics
    double lambda_free_;            // hazard rate for non-refractory neurons
    std::vector< double > m_;       // survival buffer
    std::vector< double > n_;       // population activity buffer
    std::vector< double > u_;       // mean of survivals
    std::vector< double > v_;       // variance of survivals
    std::vector< double > lambda_;  // escape rates buffer
    std::vector< double > g_;       // adaptation variables
 
    int k0_;  // rotating index of history buffers
  };
 
  // Access functions for UniversalDataLogger -----------------------
 
  double
  get_V_m_() const
  {
    return S_.V_m_;
  }
 
  double
  get_n_events_() const
  {
    return S_.n_spikes_;
  }
 
  double
  get_E_sfa_() const
  {
    return S_.theta_hat_;
  }
 
  double
  get_mean_() const
  {
    return S_.n_expect_;
  }
 
  double
  get_I_syn_ex_() const
  {
    return S_.I_syn_ex_;
  }
  double
  get_I_syn_in_() const
  {
    return S_.I_syn_in_;
  }
 
 
  // ----------------------------------------------------------------
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
  static RecordablesMap< gif_pop_psc_exp > recordablesMap_;
};
 
inline size_t
gif_pop_psc_exp::send_test_event( Node& target, size_t receptor_type, synindex, bool )
{
  SpikeEvent e;
  e.set_sender( *this );
 
  return target.handles_test_event( e, receptor_type );
}
 
inline size_t
gif_pop_psc_exp::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
gif_pop_psc_exp::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
gif_pop_psc_exp::handles_test_event( DataLoggingRequest& dlr, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}
 
inline void
gif_pop_psc_exp::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d, P_ );
  // MoD: In models derived from ArchivingNode, here get_status of the
  // parent class is called. Since this model derives from Node, and
  // not from ArchivingNode, this call has been disabled here
  // (Node does not have a comparable method).
  //  ArchivingNode::get_status(d);
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
gif_pop_psc_exp::set_status( const Dictionary& d )
{
  Parameters_ ptmp = P_;      // temporary copy in case of errors
  ptmp.set( d, this );        // throws if BadProperty
  State_ stmp = S_;           // temporary copy in case of errors
  stmp.set( d, ptmp, this );  // throws if BadProperty
 
  // We now know that (ptmp, stmp) are consistent. We do not
  // write them back to (P_, S_) before we are also sure that
  // the properties to be set in the parent class are internally
  // consistent.
 
  // MoD: In models derived from ArchivingNode, here set_status of the
  // parent class is called. Since this model derives from Node, and
  // not from ArchivingNode, this call has been disabled here
  // (Node does not have a comparable method).
  //  ArchivingNode::set_status(d);
 
  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace
 
 
#endif /* HAVE_GSL */
#endif /* #ifndef PP_POP_PSC_BETA_H */