iaf_cond_alpha_mc.h

Go to the documentation of this file.

/*
 *  iaf_cond_alpha_mc.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_ALPHA_MC_H
#define IAF_COND_ALPHA_MC_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 "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_alpha_mc_dynamics( double, const double*, double*, void* );
 
 
/* BeginUserDocs: neuron, integrate-and-fire, conductance-based, hard threshold
 
Short description
+++++++++++++++++
 
Multi-compartment conductance-based leaky integrate-and-fire neuron model
 
Description
+++++++++++
 
.. admonition:: Deprecated model
 
   ``iaf_cond_alpha_mc`` is deprecated because ``cm_default`` is an improved
   implementation of compartmental models in NEST (see :doc:`cm_default`).
 
THIS MODEL IS A PROTOTYPE FOR ILLUSTRATION PURPOSES. IT IS NOT YET
FULLY TESTED. USE AT YOUR OWN PERIL!
 
``iaf_cond_alpha_mc`` is an implementation of a multi-compartment spiking
neuron using IAF dynamics with conductance-based synapses. It serves
mainly to illustrate the implementation of :ref:`multicompartment models
<multicompartment-models>` in NEST.
 
The model has three compartments: soma, proximal and distal dendrite,
labeled as s, p, and d, respectively. Compartments are connected through
passive conductances as follows
 
.. math::
 
 C_{m.s} d/dt V_{m.s} = \ldots - g_{sp} ( V_{m.s} - V_{m.p} ) \\
 
 C_{m.p} d/dt V_{m.p} = \ldots - g_{sp} ( V_{m.p} - V_{m.s} )
     - g_{pd} ( V_{m.p} - V_{m.d} ) \\
 
 C_{m.d} d/dt V_{m.d} = \ldots \qquad - g_{pd} ( V_{m.d} - V_{m.p} )
 
A spike is fired when the somatic membrane potential exceeds threshold,
:math:`V_{m.s} >= V_{th}`. After a spike, somatic membrane potential is
clamped to a reset potential, :math:` V_{m.s} == V_{reset}`, for the refractory
period. Dendritic membrane potentials are not manipulated after a spike.
 
There is one excitatory and one inhibitory conductance-based synapse
onto each compartment, with alpha-function time course. The alpha
function is normalized such that an event of weight 1.0 results in a
peak current of 1 nS at :math:`t = \tau_{syn}`. Each compartment can also
receive current input from a current generator, and 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 ``iaf_cond_alpha`` model, all
synaptic weights must be positive numbers!
 
See also :footcite:p:`Meffin2004`, :footcite:p:`Bernander1991`.
 
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", "proximal", and "distal", 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      Conductance connecting soma and proximal dendrite
 g_pd        nS      Conductance connecting proximal and distal dendrite
 t_ref       ms      Duration of refractory period
 V_th        mV      Spike threshold in mV
 V_reset     mV      Reset potential of the membrane
============ ======= ==========================================================
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
References
++++++++++
 
 
.. footbibliography::
 
See also
++++++++
 
iaf_cond_alpha
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: iaf_cond_alpha_mc
 
EndUserDocs */
 
void register_iaf_cond_alpha_mc( const std::string& name );
 
class iaf_cond_alpha_mc : public ArchivingNode
{
 
  // Boilerplate function declarations --------------------------------
 
public:
  iaf_cond_alpha_mc();
  iaf_cond_alpha_mc( const iaf_cond_alpha_mc& );
  ~iaf_cond_alpha_mc() 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,
    PROX,
    DIST,
    NCOMP
  };
 
  static const size_t MIN_SPIKE_RECEPTOR = 1;
 
  enum SpikeSynapseTypes
  {
    SOMA_EXC = MIN_SPIKE_RECEPTOR,
    SOMA_INH,
    PROX_EXC,
    PROX_INH,
    DIST_EXC,
    DIST_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_PROX,
    I_DIST,
    SUP_CURR_RECEPTOR
  };
 
  static const size_t NUM_CURR_RECEPTORS = SUP_CURR_RECEPTOR - MIN_CURR_RECEPTOR;
 
  // Friends --------------------------------------------------------
 
  friend int iaf_cond_alpha_mc_dynamics( double, const double*, double*, void* );
 
  friend class RecordablesMap< iaf_cond_alpha_mc >;
  friend class UniversalDataLogger< iaf_cond_alpha_mc >;
 
 
  // Parameters ------------------------------------------------------
 
  struct Parameters_
  {
    double V_th;                 
    double V_reset;              
    double t_ref;                
    double g_conn[ NCOMP - 1 ];  
    double g_L[ NCOMP ];         
    double C_m[ NCOMP ];         
    double E_ex[ NCOMP ];        
    double E_in[ NCOMP ];        
    double E_L[ NCOMP ];         
    double tau_synE[ NCOMP ];    
    double tau_synI[ NCOMP ];    
    double I_e[ NCOMP ];         
 
    Parameters_();                                 
    Parameters_( const Parameters_& );             
    Parameters_& operator=( const Parameters_& );  
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
  };
 
 
  // State variables  ------------------------------------------------------
 
public:
  struct State_
  {
    enum StateVecElems_
    {
      V_M = 0,
      DG_EXC,
      G_EXC,
      DG_INH,
      G_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_&, Node* );
 
    static size_t
    idx( size_t comp, StateVecElems_ elem )
    {
      return comp * STATE_VEC_COMPS + elem;
    }
  };
 
private:
  // Internal buffers --------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( iaf_cond_alpha_mc& );  
    Buffers_( const Buffers_&, iaf_cond_alpha_mc& );
 
    UniversalDataLogger< iaf_cond_alpha_mc > 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_;  
 
    // 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_[ NCOMP ];  
  };
 
  // Internal variables ---------------------------------------------
 
  struct Variables_
  {
    double PSConInit_E_[ NCOMP ];
 
    double PSConInit_I_[ NCOMP ];
 
    int RefractoryCounts_;
  };
 
  // 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< iaf_cond_alpha_mc > recordablesMap_;
};
 
inline size_t
iaf_cond_alpha_mc::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_alpha_mc::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
iaf_cond_alpha_mc::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
iaf_cond_alpha_mc::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
iaf_cond_alpha_mc::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d );
  ArchivingNode::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::proximal_exc ] = static_cast< long >( PROX_EXC );
  receptor_dict_[ names::proximal_inh ] = static_cast< long >( PROX_INH );
  receptor_dict_[ names::proximal_curr ] = static_cast< long >( I_PROX );
 
  receptor_dict_[ names::distal_exc ] = static_cast< long >( DIST_EXC );
  receptor_dict_[ names::distal_inh ] = static_cast< long >( DIST_INH );
  receptor_dict_[ names::distal_curr ] = static_cast< long >( I_DIST );
 
  d[ names::receptor_types ] = receptor_dict_;
}
 
inline void
iaf_cond_alpha_mc::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_ALPHA_MC_H