pp_cond_exp_mc_urbanczik.h

Go to the documentation of this file.

/*
 *  pp_cond_exp_mc_urbanczik.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_COND_EXP_MC_URBANCZIK_H
#define PP_COND_EXP_MC_URBANCZIK_H
 
// Generated includes:
#include "config.h"
 
#ifdef HAVE_GSL
 
// C++ includes:
#include <vector>
 
// C includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
 
// Includes from nestkernel:
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "random_generators.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
#include "urbanczik_archiving_node.h"
#include "urbanczik_archiving_node_impl.h"
 
 
namespace nest
{
extern "C" int pp_cond_exp_mc_urbanczik_dynamics( double, const double*, double*, void* );
 
class pp_cond_exp_mc_urbanczik_parameters
{
  friend class pp_cond_exp_mc_urbanczik;
  friend class UrbanczikArchivingNode< pp_cond_exp_mc_urbanczik_parameters >;
 
private:
  enum Compartments_
  {
    SOMA = 0,
    DEND,
    NCOMP
  };
 
  double phi_max;     
  double rate_slope;  
  double beta;        
  double theta;       
  double phi( double u ) const;
  double h( double u ) const;
 
public:
  // The Urbanczik parameters need to be public within this class as they are passed to the GSL solver
  double g_conn[ NCOMP ];      
  double g_L[ NCOMP ];         
  double C_m[ NCOMP ];         
  double E_L[ NCOMP ];         
  double tau_syn_ex[ NCOMP ];  
  double tau_syn_in[ NCOMP ];  
};
 
/* BeginUserDocs: neuron, point process, conductance-based, stochastic
 
Short description
+++++++++++++++++
 
Two-compartment point process neuron with conductance-based synapses
 
Description
+++++++++++
 
pp_cond_exp_mc_urbanczik is an implementation of a two-compartment spiking
point process neuron with conductance-based synapses as it is used
in :footcite:p:`Urbanczik2014`. It is capable of connecting to an :doc:`Urbanczik synapse
<urbanczik_synapse>`.
 
The model has two compartments: soma and dendrite, labeled as s and p,
respectively. Each compartment can receive spike events and current input
from a current generator. Additionally, an external (rheobase) current can be
set for each compartment.
 
Synapses, including those for injection external currents, are addressed through
the receptor types given in the receptor_types entry of the state Dictionary.
Note that in contrast to the single-compartment models, all
synaptic weights must be positive numbers! The distinction between excitatory
and inhibitory synapses is made explicitly by specifying the receptor type of
the synapse. For example, receptor_type=dendritic_exc results in an excitatory
input and receptor_type=dendritic_inh results in an inhibitory input to the
dendritic compartment.
 
.. _multicompartment-models:
 
Multicompartment models and synaptic delays
+++++++++++++++++++++++++++++++++++++++++++
 
Note that in case of multicompartment models that represent the dendrite
explicitly, the interpretation of the synaptic delay in NEST requires careful
consideration. In NEST, the delay is at least one simulation time step and is
assumed to be located entirely at the postsynaptic side. For point neurons, it
represents the time it takes for an incoming spike to travel along the
postsynaptic dendrite before it reaches the soma, see :ref:`panel a)
<fig-multicompartment>`. Conversely, if the synaptic weight depends on the
state of the postsynaptic neuron, the delay also represents the time it takes
for the information on the state to propagate back through the dendrite to the
synapse.
 
For multicompartment models in NEST, this means the delay is positioned directly
behind the incoming synapse, that is, before the first dendritic compartment on the
postsynaptic side, see :ref:`panel b) <fig-multicompartment>`. Therefore, the
delay specified in the synapse model does *not* account for any delay that might
be associated with information traveling through the explicitly modeled
dendritic compartments.
 
In the :doc:`Urbanczik synapse <urbanczik_synapse>`, the change of the synaptic
weight is driven by an error signal, which is the difference between the firing
rate of the soma (derived from the somatic spike train :math:`S_{post}`) and the
dendritic prediction of the firing rate of the soma (derived from the dendritic
membrane potential :math:`V`). The original publication :footcite:p:`Urbanczik2014` does not assume any
delay in the interaction between the soma and the dendritic compartment.
Therefore, we evaluate the firing rate and the dendritic prediction at equal
time points to calculate the error signal at that time point. Due to the
synaptic delay :math:`d`, the synapse combines a delayed version of the error
signal with the presynaptic spike train (:math:`S_{pre}`), see :ref:`panel c)
<fig-multicompartment>`.
 
.. _fig-multicompartment:
 
.. figure:: ../static/img/multicompartment.png
   :width: 75 %
 
   a) Two point neurons (red circles *pre* and *post*) connected via a synapse.
   In NEST, the delay is entirely on the postsynaptic side, and in the case of point
   neurons, it is interpreted as the dendritic delay. b) Two two-compartment
   neuron models composed of a somatic (green) and a dendritic (blue)
   compartment. The soma of the presynaptic neuron is connected to the dendrite
   of the postsynaptic neuron. The synaptic delay is located behind the synapse
   and before the dendrite. c) Time trace of the State variables that enter the
   Urbanczik-Senn rule. Due to the synaptic delay :math:`d`, the presynaptic
   spike train (top) is combined with a delayed version of the postsynaptic
   quantities; the dendritic membrane potential (middle) and the somatic spike
   train (bottom).
 
See :doc:`../auto_examples/urbanczik_synapse_example` to learn more.
 
Parameters
++++++++++
 
The following parameters can be set in the status Dictionary. Parameters
for each compartment are collected in a sub-Dictionary; these sub-dictionaries
are called "soma" and "dendritic", respectively. In the list below,
these parameters are marked with an asterisk.
 
============   =====   =====================================================
 V_m*           mV      Membrane potential
 E_L*           mV      Leak reversal potential
 C_m*           pF      Capacity of the membrane
 E_ex*          mV      Excitatory reversal potential
 E_in*          mV      Inhibitory reversal potential
 g_L*           nS      Leak conductance
 tau_syn_ex*    ms      Rise time of the excitatory synaptic alpha function
 tau_syn_in*    ms      Rise time of the inhibitory synaptic alpha function
 I_e*           pA      Constant input current
 g_sp           nS      Coupling between soma and dendrite
 g_ps           nS      Coupling between dendrite and soma
 t_ref          ms      Duration of refractory period
============   =====   =====================================================
 
.. note::
 
   The neuron model uses standard units of NEST instead of the unitless quantities
   used in :footcite:p:`Urbanczik2014`.
 
.. note::
 
   All parameters that occur for both compartments are stored as C arrays, with
   index 0 being soma.
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
References
++++++++++
 
.. footbibliography::
 
See also
++++++++
 
urbanczik_synapse
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: pp_cond_exp_mc_urbanczik
 
EndUserDocs */
 
void register_pp_cond_exp_mc_urbanczik( const std::string& name );
 
class pp_cond_exp_mc_urbanczik : public UrbanczikArchivingNode< pp_cond_exp_mc_urbanczik_parameters >
{
 
  // Boilerplate function declarations --------------------------------
 
public:
  pp_cond_exp_mc_urbanczik();
  pp_cond_exp_mc_urbanczik( const pp_cond_exp_mc_urbanczik& );
  ~pp_cond_exp_mc_urbanczik() 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;
 
  // Enumerations and constants specifying structure and properties ----
 
  enum Compartments_
  {
    SOMA = 0,
    DEND,
    NCOMP
  };
 
  static const size_t MIN_SPIKE_RECEPTOR = 1;
 
  enum SpikeSynapseTypes
  {
    SOMA_EXC = MIN_SPIKE_RECEPTOR,
    SOMA_INH,
    DEND_EXC,
    DEND_INH,
    SUP_SPIKE_RECEPTOR
  };
 
  static const size_t NUM_SPIKE_RECEPTORS = SUP_SPIKE_RECEPTOR - MIN_SPIKE_RECEPTOR;
 
  static const size_t MIN_CURR_RECEPTOR = SUP_SPIKE_RECEPTOR;
 
  enum CurrentSynapseTypes
  {
    I_SOMA = MIN_CURR_RECEPTOR,
    I_DEND,
    SUP_CURR_RECEPTOR
  };
 
  static const size_t NUM_CURR_RECEPTORS = SUP_CURR_RECEPTOR - MIN_CURR_RECEPTOR;
 
  // Friends --------------------------------------------------------
 
  friend int pp_cond_exp_mc_urbanczik_dynamics( double, const double*, double*, void* );
 
  friend class RecordablesMap< pp_cond_exp_mc_urbanczik >;
  friend class UniversalDataLogger< pp_cond_exp_mc_urbanczik >;
 
 
  // Parameters ------------------------------------------------------
 
  struct Parameters_
  {
    double t_ref;          
    double E_ex[ NCOMP ];  
    double E_in[ NCOMP ];  
    double I_e[ NCOMP ];   
 
    pp_cond_exp_mc_urbanczik_parameters urbanczik_params;
 
    double dead_time_;
 
    Parameters_();                                 
    Parameters_( const Parameters_& );             
    Parameters_& operator=( const Parameters_& );  
 
    void get( Dictionary& ) const;  
    void set( const Dictionary& );  
  };
 
 
  // State variables  ------------------------------------------------------
 
public:
  struct State_
  {
 
    enum StateVecElems_
    {
      V_M = 0,
      G_EXC,
      G_INH,
      I_EXC,  // in the paper it is I_dnd which accounts for both excitation and inhibition
      I_INH,
      STATE_VEC_COMPS
    };
 
    static const size_t STATE_VEC_SIZE = to_underlying( STATE_VEC_COMPS ) * to_underlying( NCOMP );
 
    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_& );
 
    static size_t
    idx( size_t comp, StateVecElems_ elem )
    {
      assert( comp * STATE_VEC_COMPS + elem < STATE_VEC_SIZE );
      return comp * STATE_VEC_COMPS + elem;
    }
  };
 
private:
  // Internal buffers --------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( pp_cond_exp_mc_urbanczik& );  
    Buffers_( const Buffers_&, pp_cond_exp_mc_urbanczik& );
 
    UniversalDataLogger< pp_cond_exp_mc_urbanczik > logger_;
 
    std::vector< RingBuffer > spikes_;
    std::vector< RingBuffer > currents_;
 
    gsl_odeiv_step* s_;     
    gsl_odeiv_control* c_;  
    gsl_odeiv_evolve* e_;   
    gsl_odeiv_system sys_;  
 
    // IntegrationStep_ should be reset with the neuron on ResetNetwork,
    // but remain unchanged during calibration. Since it 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_[ NCOMP ];  
  };
 
  // Internal variables ---------------------------------------------
 
  struct Variables_
  {
    int RefractoryCounts_;
 
    double h_;                           
    RngPtr rng_;                         
    poisson_distribution poisson_dist_;  
  };
 
  // Access functions for UniversalDataLogger -------------------------------
 
  template < State_::StateVecElems_ elem, Compartments_ comp >
  double
  get_y_elem_() const
  {
    return S_.y_[ S_.idx( comp, elem ) ];
  }
 
  double
  get_r_() const
  {
    return Time::get_resolution().get_ms() * S_.r_;
  }
 
  // Data members ----------------------------------------------------
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
 
  static std::vector< std::string > comp_names_;
 
  // Dictionary of receptor types, leads to seg fault on exit, see #328
 
  static RecordablesMap< pp_cond_exp_mc_urbanczik > recordablesMap_;
};
 
 
// Inline functions of pp_cond_exp_mc_urbanczik_parameters
inline double
pp_cond_exp_mc_urbanczik_parameters::phi( double u ) const
{
  return phi_max / ( 1.0 + rate_slope * exp( beta * ( theta - u ) ) );
}
 
inline double
pp_cond_exp_mc_urbanczik_parameters::h( double u ) const
{
  return 15.0 * beta / ( 1.0 + ( 1.0 / rate_slope ) * exp( -beta * ( theta - u ) ) );
}
 
 
// Inline functions of pp_cond_exp_mc_urbanczik
inline size_t
pp_cond_exp_mc_urbanczik::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_cond_exp_mc_urbanczik::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type < MIN_SPIKE_RECEPTOR or receptor_type >= SUP_SPIKE_RECEPTOR )
  {
    if ( receptor_type >= SUP_CURR_RECEPTOR )
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
    else
    {
      throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" );
    }
  }
  return receptor_type - MIN_SPIKE_RECEPTOR;
}
 
inline size_t
pp_cond_exp_mc_urbanczik::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type < MIN_CURR_RECEPTOR or receptor_type >= SUP_CURR_RECEPTOR )
  {
    if ( receptor_type < MIN_CURR_RECEPTOR )
    {
      throw IncompatibleReceptorType( receptor_type, get_name(), "CurrentEvent" );
    }
    else
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
  }
  return receptor_type - MIN_CURR_RECEPTOR;
}
 
inline size_t
pp_cond_exp_mc_urbanczik::handles_test_event( DataLoggingRequest& dlr, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    if ( receptor_type >= SUP_CURR_RECEPTOR )
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
    else
    {
      throw IncompatibleReceptorType( receptor_type, get_name(), "DataLoggingRequest" );
    }
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}
 
inline void
pp_cond_exp_mc_urbanczik::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d );
  UrbanczikArchivingNode< pp_cond_exp_mc_urbanczik_parameters >::get_status( d );
 
  d[ names::recordables ] = recordablesMap_.get_list();
 
  Dictionary receptor_dict_;
  receptor_dict_[ names::soma_exc ] = static_cast< long >( SOMA_EXC );
  receptor_dict_[ names::soma_inh ] = static_cast< long >( SOMA_INH );
  receptor_dict_[ names::soma_curr ] = static_cast< long >( I_SOMA );
 
  receptor_dict_[ names::dendritic_exc ] = static_cast< long >( DEND_EXC );
  receptor_dict_[ names::dendritic_inh ] = static_cast< long >( DEND_INH );
  receptor_dict_[ names::dendritic_curr ] = static_cast< long >( I_DEND );
 
  d[ names::receptor_types ] = receptor_dict_;
}
 
inline void
pp_cond_exp_mc_urbanczik::set_status( const Dictionary& d )
{
  Parameters_ ptmp = P_;  // temporary copy in case of errors
  ptmp.set( d );          // throws if BadProperty
  State_ stmp = S_;       // temporary copy in case of errors
  stmp.set( d, ptmp );    // 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.
  UrbanczikArchivingNode< pp_cond_exp_mc_urbanczik_parameters >::set_status( d );
 
  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace
 
 
#endif  // HAVE_GSL
#endif  // PP_COND_EXP_MC_URBANCZIK_H