gif_cond_exp_multisynapse.h

Go to the documentation of this file.

/*
 *  gif_cond_exp_multisynapse.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 GIF_COND_EXP_MULTISYNAPSE_H
#define GIF_COND_EXP_MULTISYNAPSE_H
 
#include "config.h"
 
#ifdef HAVE_GSL
 
// Includes from gnu gsl:
#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 "ring_buffer.h"
#include "universal_data_logger.h"
 
#include "nest.h"
 
namespace nest
{
 
extern "C" int gif_cond_exp_multisynapse_dynamics( double, const double*, double*, void* );
 
/* BeginUserDocs: neuron, integrate-and-fire, conductance-based, adaptation, stochastic
 
Short description
+++++++++++++++++
 
Conductance-based generalized integrate-and-fire neuron (GIF) with multiple synaptic time constants (from the Gerstner
lab)
 
Description
+++++++++++
 
``gif_cond_exp_multisynapse`` is the generalized integrate-and-fire neuron
according to Mensi et al. (2012) :footcite:p:`Mensi2012` and Pozzorini et al. (2015) :footcite:p:`Pozzorini2015`, with
postsynaptic conductances in the form of truncated exponentials.
 
This model features both an adaptation current and a dynamic threshold for
spike-frequency adaptation. The membrane potential (V) is described by the
differential equation:
 
.. math::
 
 C \cdot dV(t)/dt = -g_L \cdot (V(t)-E_L) - \eta_1(t) - \eta_2(t) - \ldots - \eta_n(t) + I(t)
 
where each :math:`\eta_i` is a spike-triggered current (stc), and the neuron
model can have arbitrary number of them.
Dynamic of each :math:`\eta_i` is described by:
 
.. math::
 
 \tau_{\eta_i} \cdot d{\eta_i}/dt = -\eta_i
 
and in case of spike emission, its value increased by a constant (which can be
positive or negative):
 
.. math::
 
 \eta_i = \eta_i + q_{\eta_i} \text{ (in case of spike emission).}
 
Neuron produces spikes stochastically according to a point process with the
firing intensity:
 
.. math::
 
 \lambda(t) = \lambda_0 \cdot \exp(V(t)-V_T(t)) / \Delta_V
 
where :math:`V_T(t)` is a time-dependent firing threshold:
 
.. math::
 
   V_T(t) = V_{T_{star}} + \gamma_1(t) + \gamma_2(t) + \ldots + \gamma_m(t)
 
where :math:`\gamma_i` is a kernel of spike-frequency adaptation (sfa), and the
neuron model can have arbitrary number of them.
Dynamic of each :math:`\gamma_i` is described by:
 
.. math::
 
 \tau_{\gamma_i} \cdot d\gamma_i/dt = -\gamma_i
 
and in case of spike emission, its value increased by a constant (which can be
positive or negative):
 
.. math::
 
 \gamma_i = \gamma_i + q_{\gamma_i} \text{ (in case of spike emission).}
 
 
.. note::
 
   In the current implementation of the model,
   the values of :math:`\eta_i` and :math:`\gamma_i` are affected immediately
   after spike emission. However, `GIF toolbox <http://wiki.epfl.ch/giftoolbox>`_,
   which fits the model using experimental data, requires a different set of
   :math:`\eta_i` and :math:`\gamma_i`. It applies the jump of :math:`\eta_i` and
   :math:`\gamma_i` after the refractory period. One can easily convert between
   :math:`q_{\eta/\gamma}` with these two approaches:
 
   .. math::
 
      q_{\eta,giftoolbox} = q_{\eta,NEST} \cdot (1 - \exp( -\tau_{ref} / \tau_\eta ))
 
   The same formula applies for :math:`q_\gamma`.
 
On the postsynaptic side, there can be arbitrarily many synaptic time
constants (gif_psc_exp has exactly two: ``tau_syn_ex`` and
``tau_syn_in``). This can be reached by specifying separate receptor
ports, each for a different time constant. The port number has to
match the respective ``receptor_type`` in the connectors. The shape of
synaptic conductance is exponential.
 
When connecting to conductance-based multisynapse models, all synaptic weights
must be non-negative.
 
Parameters
++++++++++
 
The following parameters can be set in the status Dictionary.
 
=========   ======   ======================================================
**Membrane Parameters**
---------------------------------------------------------------------------
 C_m         pF      Capacity of the membrane
 t_ref       ms      Duration of refractory period
 V_reset     mV      Reset value for V_m after a spike
 E_L         mV      Leak reversal potential
 g_L         nS      Leak conductance
 I_e         pA      Constant external input current
=========   ======   ======================================================
 
========= ==============  =====================================================
**Spike adaptation and firing intensity parameters**
-------------------------------------------------------------------------------
 q_stc      list of nA        Values added to spike-triggered currents (stc)
                              after each spike emission
 tau_stc    list of ms        Time constants of stc variables
 q_sfa      list of mV        Values added to spike-frequency adaptation
                              (sfa) after each spike emission
 tau_sfa    list of ms        Time constants of sfa variables
 Delta_V    mV                Stochasticity level
 lambda_0   real              Stochastic intensity at firing threshold V_T i
                              n 1/s.
 V_T_star   mV                Base threshold
========= ==============  =====================================================
 
=========   =============  ===================================================
**Synaptic parameters**
------------------------------------------------------------------------------
 tau_syn    list of ms     Time constants of the synaptic conductance
                           (same size as E_rev)
 E_rev      list of mV     Reversal potentials (same size as tau_syn)
=========   =============  ===================================================
 
==============  ======  ======================================================
**Integration parameters**
------------------------------------------------------------------------------
 gsl_error_tol  real    This parameter controls the admissible error of the
                        GSL integrator. Reduce it if NEST complains about
                        numerical instabilities
==============  ======  ======================================================
 
References
++++++++++
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
See also
++++++++
 
pp_psc_delta, gif_cond_exp, iaf_psc_exp_multisynapse, gif_psc_exp_multisynapse
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: gif_cond_exp_multisynapse
 
EndUserDocs */
 
void register_gif_cond_exp_multisynapse( const std::string& name );
 
class gif_cond_exp_multisynapse : public ArchivingNode
{
 
public:
  gif_cond_exp_multisynapse();
  gif_cond_exp_multisynapse( const gif_cond_exp_multisynapse& );
  ~gif_cond_exp_multisynapse() override;
 
  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;
 
  // make dynamics function quasi-member
  friend int gif_cond_exp_multisynapse_dynamics( double, const double*, double*, void* );
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< gif_cond_exp_multisynapse >;
  friend class UniversalDataLogger< gif_cond_exp_multisynapse >;
 
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    double g_L_;
    double E_L_;
    double V_reset_;
    double Delta_V_;
    double V_T_star_;
    double lambda_0_; 
    double t_ref_;
 
    double c_m_;
 
    std::vector< double > tau_stc_;
 
    std::vector< double > q_stc_;
 
    std::vector< double > tau_sfa_;
 
    std::vector< double > q_sfa_;
 
    std::vector< double > tau_syn_;
 
    std::vector< double > E_rev_;  
 
    double I_e_;
 
    bool has_connections_;
 
    double gsl_error_tol;  
 
    Parameters_();  
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
 
    inline size_t
    n_receptors() const
    {
      return E_rev_.size();
    }
  };
 
  // ----------------------------------------------------------------
 
  struct State_
  {
    enum StateVecElems
    {
      V_M = 0,
      G,
      STATE_VEC_SIZE
    };
 
    static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 1;  
    static const size_t NUM_STATE_ELEMENTS_PER_RECEPTOR = 1;  
 
    std::vector< double > y_;  
 
    double I_stim_;  
    double sfa_;     
    double stc_;     
 
    std::vector< double > sfa_elems_;  
    std::vector< double > stc_elems_;  
 
    unsigned int r_ref_;
 
    State_( const Parameters_& );  
 
    void get( Dictionary&, const Parameters_& ) const;
    void set( const Dictionary&, const Parameters_&, Node* );
 
  };  // State_
 
  // ----------------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( gif_cond_exp_multisynapse& );
    Buffers_( const Buffers_&, gif_cond_exp_multisynapse& );
 
    std::vector< RingBuffer > spikes_;
    RingBuffer currents_;
 
    UniversalDataLogger< gif_cond_exp_multisynapse > logger_;
 
    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_;  
  };
 
  // ----------------------------------------------------------------
 
  struct Variables_
  {
    std::vector< double > P_sfa_;  
    std::vector< double > P_stc_;  
    RngPtr rng_;                   
 
    unsigned int RefractoryCounts_;
  };
 
  // Access functions for UniversalDataLogger -----------------------
 
  template < State_::StateVecElems elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }
 
  double
  get_E_sfa_() const
  {
    return S_.sfa_;
  }
 
  double
  get_I_stc_() const
  {
    return S_.stc_;
  }
 
  // ----------------------------------------------------------------
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
  static RecordablesMap< gif_cond_exp_multisynapse > recordablesMap_;
};
 
inline size_t
gif_cond_exp_multisynapse::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_cond_exp_multisynapse::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;
}
 
inline size_t
gif_cond_exp_multisynapse::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
gif_cond_exp_multisynapse::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_cond_exp_multisynapse::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d, P_ );
  ArchivingNode::get_status( d );
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
gif_cond_exp_multisynapse::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  /* #ifndef GIF_COND_EXP_MULTISYNAPSE_H */