glif_cond.h

Go to the documentation of this file.

/*
 *  glif_cond.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 GLIF_COND_H
#define GLIF_COND_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>
 
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
/* BeginUserDocs: neuron, integrate-and-fire, conductance-based, adaptation, hard threshold
 
Short description
+++++++++++++++++
 
Conductance-based generalized leaky integrate and fire (GLIF) model (from the Allen Institute)
 
Description
+++++++++++
 
``glif_cond`` provides five generalized leaky integrate
and fire (GLIF) models :footcite:p:`Teeter2018` with conductance-based synapses.
Incoming spike events induce a postsynaptic change of conductance modeled
by an alpha function :footcite:p:`Meffin2004`. The alpha function is normalized such that an event
of weight 1.0 results in a peak conductance change of 1 nS at :math:`t = \tau_\mathrm{syn}`.
On the postsynaptic side, there can be arbitrarily many synaptic time constants.
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 five GLIF models are:
 
* **GLIF Model 1** - Traditional leaky integrate and fire (LIF)
* **GLIF Model 2** - Leaky integrate and fire with biologically defined reset rules
  (LIF_R)
* **GLIF Model 3** - Leaky integrate and fire with after-spike currents (LIF_ASC)
* **GLIF Model 4** - Leaky integrate and fire with biologically defined reset rules
  and after-spike currents (LIF_R_ASC)
* **GLIF Model 5** - Leaky integrate and fire with biologically defined reset rules,
  after-spike currents and a voltage dependent threshold (LIF_R_ASC_A)
 
Remarks:
 
GLIF model mechanism setting is based on three parameters
(``spike_dependent_threshold``, ``after_spike_currents``, ``adapting_threshold``).
The settings of these three parameters for the five GLIF models are listed
below. Other combinations of these parameters will not be supported.
 
+--------+---------------------------+----------------------+--------------------+
| Model  | spike_dependent_threshold | after_spike_currents | adapting_threshold |
+========+===========================+======================+====================+
| GLIF1  | False                     | False                | False              |
+--------+---------------------------+----------------------+--------------------+
| GLIF2  | True                      | False                | False              |
+--------+---------------------------+----------------------+--------------------+
| GLIF3  | False                     | True                 | False              |
+--------+---------------------------+----------------------+--------------------+
| GLIF4  | True                      | True                 | False              |
+--------+---------------------------+----------------------+--------------------+
| GLIF5  | True                      | True                 | True               |
+--------+---------------------------+----------------------+--------------------+
 
Typical parameter setting of different levels of GLIF models for different cells
can be found and downloaded in the `Allen Cell Type Database
<https://celltypes.brain-map.org>`_. For example, the default parameter setting of this
``glif_cond`` neuron model was from the parameter values of GLIF Model 5 of Cell
490626718, which can be retrieved from the `Allen Brain Atlas
<https://celltypes.brain-map.org/mouse/experiment/electrophysiology/
490626718>`_, with units being converted from SI units (i.e., V, S (1/Ohm),
F, s, A) to NEST used units (i.e., mV, nS (1/GOhm), pF, ms, pA) and values
being rounded to appropriate digits for simplification.
 
For models with spike dependent threshold (i.e., GLIF2, GLIF4 and GLIF5),
parameter setting of ``voltage_reset_fraction`` and ``voltage_reset_add`` may lead
to the situation that voltage is bigger than threshold after reset. In this case,
the neuron will continue to spike until the end of the simulation regardless the
stimulated inputs. We recommend the setting of the parameters of these three models
to follow the condition of
 
.. math::
 
    E_L + \mathrm{voltage\_reset\_fraction} \cdot \left( V_\mathrm{th} - E_L \right)
    + \mathrm{voltage\_reset\_add} < V_\mathrm{th} + \mathrm{th\_spike\_add}
 
Parameters
++++++++++
 
The following parameters can be set in the status Dictionary.
 
========= ======== ============================================================
**Membrane parameters**
-------------------------------------------------------------------------------
V_m        double   Membrane potential in mV (absolute value)
V_th       double   Instantaneous threshold in mV
g          double   Membrane conductance in nS
E_L        double   Resting membrane potential in mV
C_m        double   Capacitance of the membrane in pF
t_ref      double   Duration of refractory time in ms
V_reset    double   Reset potential of the membrane in mV (GLIF 1 or GLIF 3)
========= ======== ============================================================
 
========================= =============== =====================================
**Spike adaptation and firing intensity parameters**
-------------------------------------------------------------------------------
th_spike_add               double         Threshold addition following spike
                                          in mV (delta_theta_s in Equation (6)
                                          in :footcite:p:`Teeter2018`)
th_spike_decay             double         Spike-induced threshold time
                                          constant in 1/ms (bs in Equation (2)
                                          in :footcite:p:`Teeter2018`)
voltage_reset_fraction     double         Voltage fraction coefficient
                                          following spike (fv in Equation (5)
                                          in :footcite:p:`Teeter2018`)
voltage_reset_add          double         Voltage addition following spike in
                                          mV (-delta_V (sign flipped) in
                                          Equation (5) in :footcite:p:`Teeter2018`)
asc_init                   double vector  Initial values of after-spike
                                          currents in pA
asc_decay                  double vector  After-spike current time constants
                                          in 1/ms (kj in Equation (3) in :footcite:p:`Teeter2018`)
asc_amps                   double vector  After-spike current amplitudes in
                                          pA (deltaIj in Equation (7) in :footcite:p:`Teeter2018`)
asc_r                      double vector  Current fraction following spike
                                          coefficients for fj in Equation (7)
                                          in :footcite:p:`Teeter2018`
th_voltage_index           double         Adaptation index of threshold - A
                                          'leak-conductance' for the
                                          voltage-dependent component of the
                                          threshold in 1/ms (av in Equation
                                          (4) in :footcite:p:`Teeter2018`)
th_voltage_decay           double         Voltage-induced threshold time
                                          constant - Inverse of which is the
                                          time constant of the
                                          voltage-dependent component of the
                                          threshold in 1/ms (bv in Equation
                                          (4) in :footcite:p:`Teeter2018`)
tau_syn                    double vector  Rise time constants of the synaptic
                                          alpha function in ms
E_rev                      double vector  Reversal potential in mV
spike_dependent_threshold  bool           flag whether the neuron has
                                          biologically defined reset rules
                                          with a spike dependent threshold
                                          component
after_spike_currents       bool           flag whether the neuron has after
                                          spike currents
adapting_threshold         bool           flag whether the neuron has a
                                          voltage dependent threshold component
========================= =============== =====================================
 
References
++++++++++
 
.. footbibliography::
 
See also
++++++++
 
gif_psc_exp_multisynapse, gif_cond_exp, gif_cond_exp_multisynapse, gif_pop_psc_exp
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: glif_cond
 
EndUserDocs */
 
namespace nest
{
 
extern "C" int glif_cond_dynamics( double, const double*, double*, void* );
 
void register_glif_cond( const std::string& name );
 
class glif_cond : public ArchivingNode
{
public:
  glif_cond();
 
  glif_cond( const glif_cond& );
 
  ~glif_cond() 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 glif_cond_dynamics( double, const double*, double*, void* );
 
  // The next three classes need to be friends to access the State_ class/member
  friend class DynamicRecordablesMap< glif_cond >;
  friend class DynamicUniversalDataLogger< glif_cond >;
  friend class DataAccessFunctor< glif_cond >;
 
  struct Parameters_
  {
    double G_;                         
    double E_L_;                       
    double th_inf_;                    
    double C_m_;                       
    double t_ref_;                     
    double V_reset_;                   
    double th_spike_add_;              
    double th_spike_decay_;            
    double voltage_reset_fraction_;    
    double voltage_reset_add_;         
    double th_voltage_index_;          
    double th_voltage_decay_;          
    std::vector< double > asc_init_;   
    std::vector< double > asc_decay_;  
    std::vector< double > asc_amps_;   
    std::vector< double > asc_r_;      
    std::vector< double > tau_syn_;    
    std::vector< double > E_rev_;      
 
    bool has_connections_;
 
    bool has_theta_spike_;
 
    bool has_asc_;
 
    bool has_theta_voltage_;
 
    size_t n_receptors_() const;  
 
    Parameters_();
 
    void get( Dictionary& ) const;
    double set( const Dictionary&, Node* );
  };
 
 
  struct State_
  {
 
    double threshold_;                  
    double threshold_spike_;            
    double threshold_voltage_;          
    std::vector< double > ASCurrents_;  
    double ASCurrents_sum_;             
    int refractory_steps_;              
 
    enum StateVecElems
    {
      V_M = 0,
      I,
      ASC_SUM,
      TH,
      TH_SPK,
      TH_VLT,
      DG_SYN,
      G_SYN,
      STATE_VECTOR_MIN_SIZE
    };
 
    static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 1;          // V_M
    static const size_t NUMBER_OF_RECORDABLES_ELEMENTS = DG_SYN - 1;  // I, ASC, TH, Th_SPK, TH_VLT
    static const size_t NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR = 2;   // DG_SYN, G_SYN
 
    std::vector< double > y_;  
 
    State_( const Parameters_& );
 
    void get( Dictionary&, const Parameters_& ) const;
    void set( const Dictionary&, const Parameters_&, double, Node* );
  };
 
 
  struct Buffers_
  {
    Buffers_( glif_cond& );
    Buffers_( const Buffers_&, glif_cond& );
 
    std::vector< RingBuffer > spikes_;  
    RingBuffer currents_;               
 
    DynamicUniversalDataLogger< glif_cond > logger_;
 
    /* GSL ODE stuff */
    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_;
  };
 
  struct Variables_
  {
    int RefractoryCounts_;                              
    double theta_spike_decay_rate_;                     
    double theta_spike_refractory_decay_rate_;          
    double theta_voltage_decay_rate_inverse_;           
    double potential_decay_rate_;                       
    double abpara_ratio_voltage_;                       
    std::vector< double > asc_decay_rates_;             
    std::vector< double > asc_stable_coeff_;            
    std::vector< double > asc_refractory_decay_rates_;  
    double phi;                                         
 
    std::vector< double > CondInitialValues_;
  };
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
 
  // Mapping of recordables names to access functions
  DynamicRecordablesMap< glif_cond > recordablesMap_;
 
  // Data Access Functor getter
  DataAccessFunctor< glif_cond > get_data_access_functor( size_t elem );
  inline double
  get_state_element( size_t elem )
  {
    if ( elem == glif_cond::State_::V_M )
    {
      return S_.y_[ elem ] + P_.E_L_;
    }
    else if ( elem == glif_cond::State_::I )
    {
      return B_.I_;
    }
    else if ( elem == glif_cond::State_::ASC_SUM )
    {
      return S_.ASCurrents_sum_;
    }
    else if ( elem == glif_cond::State_::TH )
    {
      return S_.threshold_ + P_.E_L_;
    }
    else if ( elem == glif_cond::State_::TH_SPK )
    {
      return S_.threshold_spike_;
    }
    else if ( elem == glif_cond::State_::TH_VLT )
    {
      return S_.threshold_voltage_;
    }
    else
    {
      return S_.y_[ elem - glif_cond::State_::NUMBER_OF_RECORDABLES_ELEMENTS ];
    }
  };
 
  // Utility function that inserts the synaptic conductances to the
  // recordables map
 
  std::string get_g_receptor_name( size_t receptor );
  void insert_conductance_recordables( size_t first = 0 );
};
 
 
inline size_t
glif_cond::Parameters_::n_receptors_() const
{
  return tau_syn_.size();
}
 
 
inline size_t
glif_cond::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
glif_cond::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
glif_cond::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
glif_cond::get_status( Dictionary& d ) const
{
  // get our own parameter and state data
  P_.get( d );
  S_.get( d, P_ );
 
  // get information managed by parent class
  ArchivingNode::get_status( d );
 
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
glif_cond::set_status( const Dictionary& d )
{
  Parameters_ ptmp = P_;                        // temporary copy in case of errors
  const double delta_EL = ptmp.set( d, this );  // throws if BadProperty
  State_ stmp = S_;                             // temporary copy in case of errors
  stmp.set( d, ptmp, delta_EL, this );          // throws if BadProperty
 
  ArchivingNode::set_status( d );
 
  /*
   * Here is where we must update the recordablesMap_ if new receptors
   * are added!
   */
  if ( ptmp.n_receptors_() > P_.n_receptors_() )  // Number of receptors increased
  {
    for ( size_t receptor = P_.n_receptors_(); receptor < ptmp.n_receptors_(); ++receptor )
    {
      size_t elem = glif_cond::State_::G_SYN + receptor * glif_cond::State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
      recordablesMap_.insert( get_g_receptor_name( receptor ), get_data_access_functor( elem ) );
    }
  }
  else if ( ptmp.n_receptors_() < P_.n_receptors_() )
  {  // Number of receptors decreased
    for ( size_t receptor = ptmp.n_receptors_(); receptor < P_.n_receptors_(); ++receptor )
    {
      recordablesMap_.erase( get_g_receptor_name( receptor ) );
    }
  }
 
  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace nest
 
#endif  // HAVE_GSL
#endif