hh_psc_alpha.h

Go to the documentation of this file.

/*
 *  hh_psc_alpha.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 HH_PSC_ALPHA_H
#define HH_PSC_ALPHA_H
 
// Generated includes:
#include "config.h"
 
#ifdef HAVE_GSL
 
// External includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_sf_exp.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 hh_psc_alpha_dynamics( double, const double*, double*, void* );
 
// Disable clang-formatting for documentation due to over-wide tables.
// clang-format off
/* BeginUserDocs: neuron, Hodgkin-Huxley, current-based, soft threshold
 
Short description
+++++++++++++++++
 
Hodgkin-Huxley neuron model
 
Description
+++++++++++
 
``hh_psc_alpha`` is an implementation of a spiking neuron using the Hodgkin-Huxley
formalism :footcite:p:`Hodgkin1952`, with
 
* membrane dynamics governed by sodium, potassium, and leak currents,
* :math:`\alpha`-shaped synaptic input currents,
* spike detection by a combined threshold-and-local-maximum search.
 
The Hodgkin-Huxley formalism is also described in standard textbooks
:footcite:p:`Gerstner2002`, :footcite:p:`Dayan2001`.
 
Membrane potential and channel dynamics
.......................................
 
The membrane potential evolves according to
 
.. math::
 
   C_{\text{m}} \frac{dV_\text{m}}{dt} = -(I_{\text{Na}} + I_{\text{K}} + I_{\text{L}}) + I_{\text{syn}} + I_\text{e}
 
with the voltage-gated sodium current, the delayed-rectifier potassium current, and the
leak current
 
.. math::
 
   I_{\text{Na}} &= g_{\text{Na}} \, m^3 h \, (V_\text{m} - E_{\text{Na}}) \\
   I_{\text{K}}  &= g_{\text{K}} \, n^4 \, (V_\text{m} - E_{\text{K}}) \\
   I_{\text{L}}  &= g_{\text{L}} \, (V_\text{m} - E_{\text{L}}) \;.
 
The gating variables :math:`m`, :math:`h`, and :math:`n` follow
 
.. math::
 
   \frac{dx}{dt} = \alpha_x(V_\text{m}) (1 - x) - \beta_x(V_\text{m}) \, x, \quad x \in \{m, h, n\} \;,
 
with voltage-dependent rates (:math:`V_\text{m}` in mV)
 
.. math::
 
   \alpha_m &= \frac{0.1 (V_\text{m} + 40)}{1 - e^{-(V_\text{m} + 40)/10}} \\
   \beta_m  &= 4 \, e^{-(V_\text{m} + 65)/18} \\
   \alpha_h &= 0.07 \, e^{-(V_\text{m} + 65)/20} \\
   \beta_h  &= \frac{1}{1 + e^{-(V_\text{m} + 35)/10}} \\
   \alpha_n &= \frac{0.01 (V_\text{m} + 55)}{1 - e^{-(V_\text{m} + 55)/10}} \\
   \beta_n  &= 0.125 \, e^{-(V_\text{m} + 65)/80} \;.
 
Spike detection is done by a combined threshold-and-local-maximum search: if there is a
local maximum of the membrane potential above a certain threshold, it is considered a spike.
 
Synaptic input
..............
 
The total synaptic input current :math:`I_{\text{syn}}` entering the voltage equation above has an
excitatory and an inhibitory component,
 
.. math::
 
   I_{\text{syn}}(t) = I_{\text{syn, ex}}(t) + I_{\text{syn, in}}(t) \;.
 
Each incoming spike evokes an :math:`\alpha`-shaped post-synaptic current
 
.. math::
 
   i_{\text{syn, X}}(t) = \frac{e}{\tau_{\text{syn, X}}} t e^{-\frac{t}{\tau_{\text{syn, X}}}} \Theta(t) \;,
 
for :math:`\text{X} \in \{\text{ex}, \text{in}\}`, with :math:`\Theta(x)` the Heaviside step function. The
currents are normalized such that an event of weight 1.0 evokes a peak current of 1 pA at
:math:`t = \tau_{\text{syn, X}}`.
 
Each channel current is the sum of these single-spike responses over
all presynaptic neurons :math:`j` and their spike times :math:`t_j^k`, weighted by the synaptic weights
:math:`w_j` and :math:`d_j` is the delay from neuron :math:`j`:
 
.. math::
 
   I_{\text{syn, X}}(t) = \sum_j w_j \sum_k i_{\text{syn, X}}(t - t_j^k-d_j) \;.
 
Excitatory and inhibitory inputs are distinguished by the sign of the synaptic weight.
 
For details on asynchronicity in spike and firing events with Hodgkin Huxley models
see :ref:`here <hh_details>`.
 
Parameters
++++++++++
 
The following parameters can be set in the status dictionary.
 
============== =========== ============================= ===================================================
**Parameter**  **Default** **Math equivalent**           **Description**
============== =========== ============================= ===================================================
``E_L``        -54.402 mV  :math:`E_\text{L}`            Leak reversal potential
``C_m``        100 pF      :math:`C_{\text{m}}`          Capacity of the membrane
``t_ref``      2 ms        :math:`t_{\text{ref}}`        Duration of refractory period
``g_L``        30 nS       :math:`g_\text{L}`            Leak conductance
``E_Na``       50 mV       :math:`E_{\text{Na}}`         Sodium reversal potential
``g_Na``       12000 nS    :math:`g_{\text{Na}}`         Sodium peak conductance
``E_K``        -77 mV      :math:`E_{\text{K}}`          Potassium reversal potential
``g_K``        3600 nS     :math:`g_{\text{K}}`          Potassium peak conductance
``tau_syn_ex`` 0.2 ms      :math:`\tau_{\text{syn, ex}}` Rise time of the excitatory input current
``tau_syn_in`` 2.0 ms      :math:`\tau_{\text{syn, in}}` Rise time of the inhibitory input current
``I_e``        0 pA        :math:`I_\text{e}`            Constant input current
============== =========== ============================= ===================================================
 
The following state variables evolve during simulation and are available either as neuron properties or as recordables.
 
================== =========================== ========================== =================================
**State variable** **Initial value**           **Math equivalent**        **Description**
================== =========================== ========================== =================================
``V_m``            -65 mV                      :math:`V_{\text{m}}`       Membrane potential
``Act_m``          :math:`m_\infty` [#steady]_ :math:`m`                  Sodium activation variable
``Inact_h``        :math:`h_\infty` [#steady]_ :math:`h`                  Sodium inactivation variable
``Act_n``          :math:`n_\infty` [#steady]_ :math:`n`                  Potassium activation variable
``I_syn_ex``       0 pA                        :math:`I_{\text{syn, ex}}` Excitatory synaptic input current
``I_syn_in``       0 pA                        :math:`I_{\text{syn, in}}` Inhibitory synaptic input current
================== =========================== ========================== =================================
 
.. [#steady] The gating variables are initialized to their steady-state values
   :math:`x_\infty = \alpha_x / (\alpha_x + \beta_x)` at :math:`V_\text{m} = -65` mV.
 
Problems/Todo
+++++++++++++
 
- better spike detection
- initial wavelet/spike at simulation onset
 
References
++++++++++
 
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
See also
++++++++
 
hh_cond_exp_traub
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: hh_psc_alpha
 
EndUserDocs */
// clang-format on
 
void register_hh_psc_alpha( const std::string& name );
 
class hh_psc_alpha : public ArchivingNode
{
 
public:
  hh_psc_alpha();
  hh_psc_alpha( const hh_psc_alpha& );
  ~hh_psc_alpha() 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;
 
  // END Boilerplate function declarations ----------------------------
 
  // Friends --------------------------------------------------------
 
  // make dynamics function quasi-member
  friend int hh_psc_alpha_dynamics( double, const double*, double*, void* );
 
  // The next two classes need to be friend to access the State_ class/member
  friend class RecordablesMap< hh_psc_alpha >;
  friend class UniversalDataLogger< hh_psc_alpha >;
 
private:
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    double t_ref_;    
    double g_Na;      
    double g_K;       
    double g_L;       
    double C_m;       
    double E_Na;      
    double E_K;       
    double E_L;       
    double tau_synE;  
    double tau_synI;  
    double I_e;       
 
    Parameters_();  
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
  };
 
public:
  // ----------------------------------------------------------------
 
  struct State_
  {
    enum StateVecElems
    {
      V_M = 0,
      HH_M,    // 1
      HH_H,    // 2
      HH_N,    // 3
      DI_EXC,  // 4
      I_EXC,   // 5
      DI_INH,  // 6
      I_INH,   // 7
      STATE_VEC_SIZE
    };
 
 
    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&, Node* node );
  };
 
  // ----------------------------------------------------------------
 
private:
  struct Buffers_
  {
    Buffers_( hh_psc_alpha& );                   
    Buffers_( const Buffers_&, hh_psc_alpha& );  
 
    UniversalDataLogger< hh_psc_alpha > logger_;
 
    RingBuffer spike_exc_;
    RingBuffer spike_inh_;
    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_;
  };
 
  // ----------------------------------------------------------------
 
  struct Variables_
  {
    double PSCurrInit_E_;
 
    double PSCurrInit_I_;
 
    int RefractoryCounts_;
  };
 
  // Access functions for UniversalDataLogger -------------------------------
 
  template < State_::StateVecElems elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }
 
  // ----------------------------------------------------------------
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
 
  static RecordablesMap< hh_psc_alpha > recordablesMap_;
};
 
 
inline size_t
hh_psc_alpha::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
hh_psc_alpha::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
hh_psc_alpha::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
hh_psc_alpha::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
hh_psc_alpha::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d );
  ArchivingNode::get_status( d );
 
  d[ names::recordables ] = recordablesMap_.get_list();
}
 
inline void
hh_psc_alpha::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, 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  // HH_PSC_ALPHA_H