glif_psc.h

Go to the documentation of this file.

/*
 *  glif_psc.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_PSC_H
#define GLIF_PSC_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, current-based, adaptation, hard threshold
 
Short description
+++++++++++++++++
 
Current-based generalized leaky integrate-and-fire (GLIF) models (from the Allen Institute)
 
Description
+++++++++++
 
``glif_psc`` provides five generalized leaky integrate-and-fire
(GLIF) models :footcite:p:`Teeter2018` with alpha-function shaped synaptic currents.
Incoming spike events induce a postsynaptic change of current 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 current of 1 pA at :math:`t = \tau_\mathrm{syn}`.
By default, ``glif_psc`` has a single synapse that is accessible through
``receptor_port`` 1. An arbitrary number of synapses with different time constants
can be configured by setting the desired time constants as ``tau_syn`` array.
The resulting synapses are addressed through ``receptor_port`` 1, 2, 3, ....
 
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_psc`` 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}
 
.. note::
 
  If ``tau_m`` is very close to ``tau_syn_ex`` or ``tau_syn_in``, the model
  will numerically behave as if ``tau_m`` is equal to ``tau_syn_ex`` or
  ``tau_syn_in``, respectively, to avoid numerical instabilities.
 
  For implementation details see the
  `IAF Integration Singularity notebook <../model_details/IAF_Integration_Singularity.ipynb>`_.
 
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  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,
glif_psc_double_alpha
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: glif_psc
 
EndUserDocs */
 
namespace nest
{
 
void register_glif_psc( const std::string& name );
 
class glif_psc : public ArchivingNode
{
public:
  glif_psc();
 
  glif_psc( const glif_psc& );
 
  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;
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< glif_psc >;
  friend class UniversalDataLogger< glif_psc >;
 
  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_;    
 
    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 U_;                          
    double threshold_;                  
    double threshold_spike_;            
    double threshold_voltage_;          
    double I_;                          
    double I_syn_;                      
    std::vector< double > ASCurrents_;  
    double ASCurrents_sum_;             
    int refractory_steps_;              
    std::vector< double > y1_;          
    std::vector< double > y2_;          
 
    State_( const Parameters_& );
 
    void get( Dictionary&, const Parameters_& ) const;
    void set( const Dictionary&, const Parameters_&, double, Node* );
  };
 
 
  struct Buffers_
  {
    Buffers_( glif_psc& );
    Buffers_( const Buffers_&, glif_psc& );
 
    std::vector< RingBuffer > spikes_;  
    RingBuffer currents_;               
 
    UniversalDataLogger< glif_psc > logger_;
  };
 
  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 > P11_;  
    std::vector< double > P21_;  
    std::vector< double > P22_;  
    double P30_;                 
    double P33_;                 
    std::vector< double > P31_;  
    std::vector< double > P32_;  
 
    std::vector< double > PSCInitialValues_;
  };
 
  double
  get_V_m_() const
  {
    return S_.U_ + P_.E_L_;
  }
 
  double
  get_ASCurrents_sum_() const
  {
    return S_.ASCurrents_sum_;
  }
 
  double
  get_I_() const
  {
    return S_.I_;
  }
 
  double
  get_I_syn_() const
  {
    return S_.I_syn_;
  }
 
  double
  get_threshold_() const
  {
    return S_.threshold_ + P_.E_L_;
  }
 
  double
  get_threshold_spike_() const
  {
    return S_.threshold_spike_;
  }
 
  double
  get_threshold_voltage_() const
  {
    return S_.threshold_voltage_;
  }
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
 
  static RecordablesMap< glif_psc > recordablesMap_;
};
 
 
inline size_t
glif_psc::Parameters_::n_receptors_() const
{
  return tau_syn_.size();
}
 
inline size_t
glif_psc::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_psc::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
glif_psc::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_psc::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_psc::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 );
 
  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
 
}  // namespace nest
 
#endif