pp_psc_delta.h

Go to the documentation of this file.

/*
 *  pp_psc_delta.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_PSC_DELTA_H
#define PP_PSC_DELTA_H
 
// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "random_generators.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
namespace nest
{
 
/* BeginUserDocs: neuron, point process, current-based, adaptation, stochastic
 
Short description
+++++++++++++++++
 
Point process neuron with leaky integration of delta-shaped PSCs
 
Description
+++++++++++
 
``pp_psc_delta`` is an implementation of a leaky integrator, where the potential
jumps on each spike arrival. It produces spike stochastically, and supports
spike-frequency adaptation, and other optional features.
 
Spikes are generated randomly according to the current value of the
transfer function which operates on the membrane potential. Spike
generation is followed by an optional dead time. Setting with_reset to
true will reset the membrane potential after each spike.
 
The transfer function can be chosen to be linear, exponential or a sum of
both by adjusting three parameters:
 
.. math::
 
  rate = Rect[ c_1 \cdot V' + c_2 \cdot \exp(c_3 * V') ],
 
where the effective potential :math:`V' = V_m - E_{sfa}` and :math:`E_{sfa}`
is called the adaptive threshold. Here Rect means rectifier:
:math:`Rect(x) = {x \text{ if } x>=0, 0 \text{ else}}` (this is necessary
because
negative rates are not possible).
 
By setting c_3 = 0, c_2 can be used as an offset spike rate for an otherwise
linear rate model.
 
The dead time enables to include refractoriness. If dead time is 0, the
number of spikes in one time step might exceed one and is drawn from the
Poisson distribution accordingly. Otherwise, the probability for a spike
is given by :math:`1 - \exp(-rate \cdot h)`, where h is the simulation time step. If
dead_time is smaller than the simulation resolution (time step), it is
internally set to the resolution.
 
Note that, even if non-refractory neurons are to be modeled, a small value
of dead_time, like ``dead_time=1e-8``, might be the value of choice since it
uses faster uniform random numbers than ``dead_time=0``, which draws Poisson
numbers. Only for very large spike rates (> 1 spike/time_step) this will
cause errors.
 
The model can optionally include an adaptive firing threshold.
If the neuron spikes, the threshold increases and the membrane potential
will take longer to reach it.
Here this is implemented by subtracting the value of the adaptive threshold
E_sfa from the membrane potential ``V_m`` before passing the potential to the
transfer function, see also above. ``E_sfa`` jumps by ``q_sfa`` when the neuron
fires a spike, and decays exponentially with the time constant tau_sfa
after (see :footcite:p:`Jolivet2006` or :footcite:p:`Pozzorini2013`). Thus, the ``E_sfa`` corresponds to the convolution
of the neuron's spike train with an exponential kernel. This adaptation kernel may also be chosen as the sum of n
exponential kernels. To use this feature,`` ``q_sfa and ``tau_sfa`` have to be given as a list of n values each.
 
The firing of ``pp_psc_delta`` is usually not a renewal process. For example,
its firing may depend on its past spikes if it has non-zero adaptation terms
(q_sfa). But if so, it will depend on all its previous spikes, not just the
last one -- so it is not a renewal process model. However, if ``with_reset``
is True, and all adaptation terms (``q_sfa``) are 0, then it will reset
("forget") its membrane potential each time a spike is emitted, which makes
it a renewal process model (where ``rate`` above is its hazard function,
also known as conditional intensity).
 
``pp_psc_delta`` may also be called a spike-response model with escape-noise :footcite:p:`Gerstner2014`
(for vanishing, non-random dead_time). If ``c_1>0`` and ``c_2==0``, the rate is a
convolution of the inputs with exponential filters -- which is a model known
as a Hawkes point process (see :footcite:p:`Grytskyy2013`). If instead ``c_1==0``, then ``pp_psc_delta`` is
a point process generalized linear model (with the canonical link function,
and exponential input filters) (see [5,6]_).
 
This model has been adapted from ``iaf_psc_delta``. The default parameters are
set to the mean values given in :footcite:p:`Jolivet2006`, which have been matched to spike-train
recordings. Due to the many features of ``pp_psc_delta`` and its versatility,
parameters should be set carefully and consciously.
 
See also :footcite:p:`Cardanobile2010`, :footcite:p:`Deger2014`.
 
Parameters
++++++++++
 
The following parameters can be set in the status dictionary.
 
 
=================  ======= ===================================================
 V_m               mV      Membrane potential
 C_m               pF      Capacitance of the membrane
 tau_m             ms      Membrane time constant
 q_sfa             mV      Adaptive threshold jump
 tau_sfa           ms      Adaptive threshold time constant
 dead_time         ms      Duration of the dead time
 dead_time_random  boolean Should a random dead time be drawn after each
                           spike?
 dead_time_shape   integer Shape parameter of dead time gamma distribution
 t_ref_remaining   ms      Remaining dead time at simulation start
 with_reset        boolean Should the membrane potential be reset after a
                           spike?
 I_e               pA      Constant input current
 c_1               Hz/mV   Slope of linear part of transfer function in
                           Hz/mV
 c_2               Hz      Prefactor of exponential part of transfer function
 c_3               1/mV    Coefficient of exponential non-linearity of
                           transfer function
=================  ======= ===================================================
 
 
References
++++++++++
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
See also
++++++++
 
iaf_psc_delta, iaf_psc_alpha, iaf_psc_exp, iaf_psc_delta_ps
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: pp_psc_delta
 
EndUserDocs */
 
void register_pp_psc_delta( const std::string& name );
 
class pp_psc_delta : public ArchivingNode
{
 
public:
  pp_psc_delta();
  pp_psc_delta( const pp_psc_delta& );
 
  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_state_() override;
  void init_buffers_() override;
  void pre_run_hook() override;
 
  void update( Time const&, const long, const long ) override;
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< pp_psc_delta >;
  friend class UniversalDataLogger< pp_psc_delta >;
 
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    double tau_m_;
 
    double c_m_;
 
    double dead_time_;
 
    bool dead_time_random_;
 
    long dead_time_shape_;
 
    bool with_reset_;
 
    std::vector< double > tau_sfa_;
 
    std::vector< double > q_sfa_;
 
    bool multi_param_;
 
    double c_1_;
 
    double c_2_;
 
    double c_3_;
 
    double I_e_;
 
    double t_ref_remaining_;
 
    Parameters_();  
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
  };
 
  // ----------------------------------------------------------------
 
  struct State_
  {
    double y0_;  
    double y3_;
    double q_;  
 
    std::vector< double > q_elems_;
 
    int r_;  
 
    bool initialized_;  
 
    State_();  
 
    void get( Dictionary&, const Parameters_& ) const;
    void set( const Dictionary&, const Parameters_&, Node* );
  };
 
  // ----------------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( pp_psc_delta& );
    Buffers_( const Buffers_&, pp_psc_delta& );
 
    RingBuffer spikes_;
    RingBuffer currents_;
 
    UniversalDataLogger< pp_psc_delta > logger_;
  };
 
  // ----------------------------------------------------------------
 
  struct Variables_
  {
 
    double P30_;
    double P33_;
 
    std::vector< double > Q33_;
 
    double h_;        
    double dt_rate_;  
 
    RngPtr rng_;                         
    gamma_distribution gamma_dist_;      
    poisson_distribution poisson_dist_;  
 
    int DeadTimeCounts_;
  };
 
  // Access functions for UniversalDataLogger -----------------------
 
  double
  get_V_m_() const
  {
    return S_.y3_;
  }
 
  double
  get_E_sfa_() const
  {
    return S_.q_;
  }
 
  // ----------------------------------------------------------------
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
  static RecordablesMap< pp_psc_delta > recordablesMap_;
};
 
inline size_t
pp_psc_delta::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
pp_psc_delta::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
pp_psc_delta::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
pp_psc_delta::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
pp_psc_delta::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d, P_ );
  ArchivingNode::get_status( d );
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
pp_psc_delta::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 /* #ifndef PP_PSC_DELTA_H */