iaf_cond_beta.h

Go to the documentation of this file.

/*
 *  iaf_cond_beta.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 IAF_COND_BETA_H
#define IAF_COND_BETA_H
 
// Generated includes:
#include "config.h"
 
#ifdef HAVE_GSL
 
// C includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
 
// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
namespace nest
{
extern "C" int iaf_cond_beta_dynamics( double, const double*, double*, void* );
 
// Disable clang-formatting for documentation due to over-wide tables.
// clang-format off
/* BeginUserDocs: neuron, integrate-and-fire, conductance-based, hard threshold
 
Short description
+++++++++++++++++
 
Leaky integrate-and-fire neuron model with beta-shaped synaptic conductances
 
Description
+++++++++++
 
``iaf_cond_beta`` is a leaky integrate-and-fire neuron model with
 
* a hard threshold,
* a fixed refractory period,
* beta-shaped synaptic conductances (a difference of two exponentials) with separate rise and decay
  times, normalized such that an event of weight 1.0 results in a peak conductance of 1 nS.
 
The model follows the conductance-based integrate-and-fire framework discussed in
:footcite:p:`Meffin2004`, :footcite:p:`Bernander1991`, :footcite:p:`Kuhn2004`, :footcite:p:`Rotter1999`,
:footcite:p:`Roth2010`.
 
Membrane potential evolution, spike emission, and refractoriness
................................................................
 
The membrane potential evolves according to
 
.. math::
 
   \frac{dV_\text{m}}{dt} = \frac{ -g_{\text{L}} (V_{\text{m}} - E_{\text{L}}) - I_{\text{syn}} + I_\text{e} } {C_{\text{m}}}
 
where the synaptic input current :math:`I_{\text{syn}}(t)` is discussed below and :math:`I_\text{e}` is
a constant input current set as a model parameter.
 
A spike is emitted at time step :math:`t^*=t^{k+1}` if
 
.. math::
 
   V_\text{m}(t_k) < V_{th} \quad\text{and}\quad V_\text{m}(t^{k+1})\geq V_\text{th} \;.
 
Subsequently,
 
.. math::
 
   V_\text{m}(t) = V_{\text{reset}} \quad\text{for}\quad t^* \leq t \leq t^* + t_{\text{ref}} \;,
 
that is, the membrane potential is clamped to :math:`V_{\text{reset}}` during the refractory period.
 
Synaptic input
..............
 
The synaptic input current has an excitatory and an inhibitory component
 
.. math::
 
   I_{\text{syn}}(t) = I_{\text{syn, ex}}(t) + I_{\text{syn, in}}(t)
 
where
 
.. math::
 
   I_{\text{syn, X}}(t) = (V_{\text{m}}(t) - E_{\text{syn, X}}) \sum_{j} w_j \sum_k g_{\text{X}}(t-t_j^k-d_j) \;,
 
where :math:`j` indexes either excitatory (:math:`\text{X} = \text{ex}`)
or inhibitory (:math:`\text{X} = \text{in}`) presynaptic neurons,
:math:`k` indexes the spike times of neuron :math:`j`, and :math:`d_j`
is the delay from neuron :math:`j`.
 
The individual synaptic conductances are beta-shaped, that is, given by the difference of two
exponentials with separate rise and decay time constants,
 
.. math::
 
   g_{\text{X}}(t) =
   \frac{ e^{-\frac{t}{\tau_{\text{decay, X}}}} - e^{-\frac{t}{\tau_{\text{rise, X}}}} }
        { e^{-\frac{t_{\text{peak, X}}}{\tau_{\text{decay, X}}}} - e^{-\frac{t_{\text{peak, X}}}{\tau_{\text{rise, X}}}} }
   \, \Theta(t)
 
where :math:`\Theta(x)` is the Heaviside step function and
 
.. math::
 
   t_{\text{peak, X}} = \frac{ \tau_{\text{decay, X}} \, \tau_{\text{rise, X}} }
                             { \tau_{\text{decay, X}} - \tau_{\text{rise, X}} }
                        \ln \frac{ \tau_{\text{decay, X}} }{ \tau_{\text{rise, X}} }
 
is the time at which the conductance attains its maximum. The conductances are normalized to
unit peak, that is,
 
.. math::
 
   g_{\text{X}}(t = t_{\text{peak, X}}) = 1 \;,
 
as the weight of the connection :math:`w_{j}` from neuron :math:`j` (excitatory if
:math:`w_{j} > 0` or inhibitory if :math:`w_{j} < 0`) is written explicitly in the equation for the total
synaptic current.
 
For :math:`\tau_{\text{rise, X}} = \tau_{\text{decay, X}}`, which is the case for the default
parameters, the beta function reduces to the alpha function of ``iaf_cond_alpha``, and the model
uses that form.
 
.. note::
 
   This class distinguishes between excitatory and inhibitory inputs by the sign of the
   synaptic weight. For input separation by ``receptor_type``, see ``iaf_cond_alpha_mc``.
   This model integrates the dynamics with an adaptive-step-size Runge-Kutta-Fehlberg
   solver from the GNU Scientific Library.
 
Parameters
++++++++++
 
The following parameters can be set in the status dictionary.
 
================ =========== =============================== ===================================================
**Parameter**    **Default** **Math equivalent**             **Description**
================ =========== =============================== ===================================================
``E_L``          -70 mV      :math:`E_\text{L}`              Leak reversal potential
``C_m``          250 pF      :math:`C_{\text{m}}`            Capacity of the membrane
``t_ref``        2 ms        :math:`t_{\text{ref}}`          Duration of refractory period
``V_th``         -55 mV      :math:`V_{\text{th}}`           Spike threshold
``V_reset``      -60 mV      :math:`V_{\text{reset}}`        Reset potential of the membrane
``E_ex``         0 mV        :math:`E_\text{ex}`             Excitatory reversal potential
``E_in``         -85 mV      :math:`E_\text{in}`             Inhibitory reversal potential
``g_L``          16.6667 nS  :math:`g_\text{L}`              Leak conductance
``tau_rise_ex``  0.2 ms      :math:`\tau_{\text{rise, ex}}`  Rise time of the excitatory synaptic conductance
``tau_decay_ex`` 0.2 ms      :math:`\tau_{\text{decay, ex}}` Decay time of the excitatory synaptic conductance
``tau_rise_in``  2.0 ms      :math:`\tau_{\text{rise, in}}`  Rise time of the inhibitory synaptic conductance
``tau_decay_in`` 2.0 ms      :math:`\tau_{\text{decay, in}}` Decay time of the inhibitory synaptic conductance
``I_e``          0 pA        :math:`I_\text{e}`              Constant input current
================ =========== =============================== ===================================================
 
The following state variables evolve during simulation and are available either as neuron properties or as recordables.
 
================== ================= ========================== =================================
**State variable** **Initial value** **Math equivalent**        **Description**
================== ================= ========================== =================================
``V_m``            -70 mV            :math:`V_{\text{m}}`       Membrane potential
``g_ex``           0 nS              :math:`g_{\text{ex}}`      Excitatory synaptic conductance
``g_in``           0 nS              :math:`g_{\text{in}}`      Inhibitory synaptic conductance
================== ================= ========================== =================================
 
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
References
++++++++++
 
.. footbibliography::
 
See also
++++++++
 
iaf_cond_exp, iaf_cond_alpha, iaf_cond_alpha_mc
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: iaf_cond_beta
 
EndUserDocs */
// clang-format on
 
void register_iaf_cond_beta( const std::string& name );
 
class iaf_cond_beta : public ArchivingNode
{
 
  // Boilerplate function declarations --------------------------------
 
public:
  iaf_cond_beta();
  iaf_cond_beta( const iaf_cond_beta& );
  ~iaf_cond_beta() override;
 
  /*
   * Import all overloaded virtual functions that we
   * override in this class.  For background information,
   * see http://www.gotw.ca/gotw/005.htm.
   */
 
  using Node::handle;
  using Node::handles_test_event;
 
  size_t send_test_event( Node& tagret, size_t receptor_type, synindex, bool ) 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 handle( SpikeEvent& ) override;
  void handle( CurrentEvent& ) override;
  void handle( DataLoggingRequest& ) override;
 
  void get_status( Dictionary& ) const override;
  void set_status( const Dictionary& ) override;
 
private:
  void init_buffers_() override;
  double get_normalisation_factor( double, double );
  void pre_run_hook() override;
  void update( Time const&, const long, const long ) override;
 
  // END Boilerplate function declarations ----------------------------
 
  // Friends --------------------------------------------------------
 
  // make dynamics function quasi-member
  friend int iaf_cond_beta_dynamics( double, const double*, double*, void* );
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< iaf_cond_beta >;
  friend class UniversalDataLogger< iaf_cond_beta >;
 
private:
  // Parameters class -------------------------------------------------
 
  struct Parameters_
  {
    double V_th;          
    double V_reset;       
    double t_ref;         
    double g_L;           
    double C_m;           
    double E_ex;          
    double E_in;          
    double E_L;           
    double tau_rise_ex;   
    double tau_decay_ex;  
    double tau_rise_in;   
    double tau_decay_in;  
    double I_e;           
 
    Parameters_();  
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
  };
 
  // State variables class --------------------------------------------
 
public:
  struct State_
  {
    enum StateVecElems
    {
      V_M = 0,
      DG_EXC,
      G_EXC,
      DG_INH,
      G_INH,
      STATE_VEC_SIZE
    };
 
    double y[ STATE_VEC_SIZE ];
 
    int r;
 
    State_( const Parameters_& );  
    State_( const State_& );
 
    State_& operator=( const State_& );
 
    void get( Dictionary& ) const;  
 
    void set( const Dictionary&, const Parameters_&, Node* );
  };
 
private:
  // Buffers class --------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( iaf_cond_beta& );                   
    Buffers_( const Buffers_&, iaf_cond_beta& );  
 
    UniversalDataLogger< iaf_cond_beta > logger_;
 
    RingBuffer spike_exc_;
    RingBuffer spike_inh_;
    RingBuffer currents_;
 
    /* GSL ODE stuff */
    gsl_odeiv_step* s_;     
    gsl_odeiv_control* c_;  
    gsl_odeiv_evolve* e_;   
    gsl_odeiv_system sys_;  
 
    // Since IntegrationStep_ is initialized with step_, and the resolution
    // cannot change after nodes have been created, it is safe to place both
    // here.
    double step_;             
    double IntegrationStep_;  
 
    double I_stim_;
  };
 
  // Variables class -------------------------------------------------------
 
  struct Variables_
  {
    double PSConInit_E;
 
    double PSConInit_I;
 
    int RefractoryCounts;
  };
 
  // Access functions for UniversalDataLogger -------------------------------
 
  template < State_::StateVecElems elem >
  double
  get_y_elem_() const
  {
    return S_.y[ elem ];
  }
 
  double
  get_r_() const
  {
    return Time::get_resolution().get_ms() * S_.r;
  }
 
  // Data members -----------------------------------------------------------
 
  // keep the order of these lines, seems to give best performance
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
 
  static RecordablesMap< iaf_cond_beta > recordablesMap_;
};
 
 
// Boilerplate inline function definitions ----------------------------------
 
inline size_t
iaf_cond_beta::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
iaf_cond_beta::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
iaf_cond_beta::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
iaf_cond_beta::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
iaf_cond_beta::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d );
  ArchivingNode::get_status( d );
 
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
iaf_cond_beta::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.
  ArchivingNode::set_status( d );
 
  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace
 
#endif  // HAVE_GSL
 
#endif  // IAF_COND_BETA_H