hh_cond_beta_gap_traub.h

Go to the documentation of this file.

/*
 *  hh_cond_beta_gap_traub.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_COND_BETA_GAP_TRAUB_H
#define HH_COND_BETA_GAP_TRAUB_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 <gsl/gsl_sf_exp.h>
 
// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "node.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
namespace nest
{
 
extern "C" int hh_cond_beta_gap_traub_dynamics( double, const double*, double*, void* );
 
/* BeginUserDocs: neuron, Hodgkin-Huxley, conductance-based, soft threshold
 
Short description
+++++++++++++++++
 
Hodgkin-Huxley neuron with gap junction support and beta function synaptic conductances
 
Description
+++++++++++
 
``hh_cond_beta_gap_traub`` is an implementation of a modified Hodgkin-Huxley model
that also supports gap junctions.
 
This model is derived from the ``hh_conda_exp`` model, but supports double-exponential-shaped
(beta-shaped) synaptic conductances and also supports gap junctions. The model is originally
based on a model of hippocampal pyramidal cells by Traub and Miles :footcite:p:`Traub1991`.
The key differences between the current model and the model in :footcite:p:`Traub1991` are:
 
- This model is a point neuron, not a compartmental model.
- Following :footcite:p:`Brette2007`, this model includes only ``I_Na`` and ``I_K``, with simpler ``I_K`` dynamics than
  in :footcite:p:`Traub1991`, so it has only three instead of eight gating variables;
  in particular, all Ca dynamics have been removed.
- Incoming spikes induce an instantaneous conductance change followed by
  exponential decay instead of activation over time.
- The model incorporates gap junctions :footcite:p:`Hahne2015`.
 
For details on asynchronicity in spike and firing events with Hodgkin Huxley models
see :ref:`here <hh_details>`.
 
 
Postsynaptic currents
---------------------
 
Incoming spike events induce a postsynaptic change of conductance modelled by a
beta function as outlined in :footcite:p:`Rotter1999` :footcite:p:`Roth2010`. The beta function is normalized such that
an event of weight 1.0 results in a peak conductance of 1 nS at :math:`t = \tau_{rise,xx}` where xx is `ex` or `in`.
 
Spike Detection
---------------
 
Spike detection is done by a combined threshold-and-local-maximum search: if
there is a local maximum above a certain threshold of the membrane potential,
it is considered a spike.
 
Gap Junctions
-------------
 
Gap Junctions are implemented by a gap current of the form
:math:`g_{ij}( V_i - V_j)`.
 
.. note::
   In this model, a spike is emitted if :math:`V_m \geq V_T + 30` mV and
   :math:`V_m` has fallen during the current time step.
 
   To avoid multiple spikes from occurring during the falling flank of a
   spike, it is essential to choose a sufficiently long refractory period.
   Traub and Miles used :math:`t_{ref} = 3` ms (:footcite:p:`Traub1991`, p 118), while we used
   :math:`t_{ref} = 2` ms in :footcite:p:`Traub1991`.
 
Parameters
++++++++++
 
The following parameters can be set in the status Dictionary.
 
============ ======  =======================================================
V_m          mV      Membrane potential
V_T          mV      Voltage offset that controls dynamics. For default
                     parameters, V_T = -63mV results in a threshold around
                     -50mV
E_L          mV      Leak reversal potential
C_m          pF      Capacity of the membrane
g_L          nS      Leak conductance
tau_rise_ex  ms      Excitatory synaptic beta function rise time
tau_decay_ex ms      Excitatory synaptic beta function decay time
tau_rise_in  ms      Inhibitory synaptic beta function rise time
tau_decay_in ms      Inhibitory synaptic beta function decay time
t_ref        ms      Duration of refractory period (see Note)
E_ex         mV      Excitatory synaptic reversal potential
E_in         mV      Inhibitory synaptic reversal potential
E_Na         mV      Sodium reversal potential
g_Na         nS      Sodium peak conductance
E_K          mV      Potassium reversal potential
g_K          nS      Potassium peak conductance
I_e          pA      External input current
============ ======  =======================================================
 
References
++++++++++
 
.. footbibliography::
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
See also
++++++++
 
hh_psc_alpha_gap, hh_cond_exp_traub, gap_junction, iaf_cond_beta
 
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: hh_cond_beta_gap_traub
 
EndUserDocs */
 
void register_hh_cond_beta_gap_traub( const std::string& name );
 
class hh_cond_beta_gap_traub : public ArchivingNode
{
 
public:
  typedef Node base;
 
  hh_cond_beta_gap_traub();
  hh_cond_beta_gap_traub( const hh_cond_beta_gap_traub& );
  ~hh_cond_beta_gap_traub() override;
 
  using Node::handle;
  using Node::handles_test_event;
  using Node::sends_secondary_event;
 
  size_t send_test_event( Node& target, size_t receptor_type, synindex, bool ) override;
 
  void handle( SpikeEvent& ) override;
  void handle( CurrentEvent& ) override;
  void handle( DataLoggingRequest& ) override;
  void handle( GapJunctionEvent& ) 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;
  size_t handles_test_event( GapJunctionEvent&, size_t ) override;
 
  void
  sends_secondary_event( GapJunctionEvent& ) override
  {
  }
 
  void get_status( Dictionary& ) const override;
  void set_status( const Dictionary& ) override;
 
private:
  void init_buffers_() override;
  double get_normalisation_factor( double, double );
  void pre_run_hook() override;
 
  bool update_( Time const&, const long, const long, const bool );
 
  void update( Time const&, const long, const long ) override;
  bool wfr_update( Time const&, const long, const long ) override;
 
  // END Boilerplate function declarations ----------------------------
 
  // Friends --------------------------------------------------------
 
  // make dynamics function quasi-member
  friend int hh_cond_beta_gap_traub_dynamics( double, const double*, double*, void* );
 
  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< hh_cond_beta_gap_traub >;
  friend class UniversalDataLogger< hh_cond_beta_gap_traub >;
 
private:
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    double g_Na;          
    double g_K;           
    double g_L;           
    double C_m;           
    double E_Na;          
    double E_K;           
    double E_L;           
    double V_T;           
    double E_ex;          
    double E_in;          
    double tau_rise_ex;   
    double tau_decay_ex;  
    double tau_rise_in;   
    double tau_decay_in;  
    double t_ref_;        
    double I_e;           
 
    Parameters_();
 
    void get( Dictionary& ) const;         
    void set( const Dictionary&, Node* );  
  };
 
public:
  // ----------------------------------------------------------------
 
  struct State_
  {
 
    enum StateVecElems
    {
      V_M = 0,
      HH_M,    // 1
      HH_H,    // 2
      HH_N,    // 3
      DG_EXC,  // 4
      G_EXC,   // 5
      DG_INH,  // 6
      G_INH,   // 7
      STATE_VEC_SIZE
    };
 
    double y_[ STATE_VEC_SIZE ];
    int r_;  
 
    State_( const Parameters_& p );
    State_( const State_& s );
 
    State_& operator=( const State_& );
 
    void get( Dictionary& ) const;
    void set( const Dictionary&, const Parameters_&, Node* );
  };
 
  // Variables class -------------------------------------------------------
 
  struct Variables_
  {
    double PSConInit_E;
 
    double PSConInit_I;
 
    int refractory_counts_;
    double U_old_;  // for spike-detection
  };
 
  // ----------------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( hh_cond_beta_gap_traub& );  
    Buffers_( const Buffers_&, hh_cond_beta_gap_traub& );
 
    UniversalDataLogger< hh_cond_beta_gap_traub > 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_;  
 
    // remembers current lag for piecewise interpolation
    long lag_;
 
    // remembers y_values from last wfr_update
    std::vector< double > last_y_values;
 
    // summarized gap weight
    double sumj_g_ij_;
 
    // summarized coefficients of the interpolation polynomial
    std::vector< double > interpolation_coefficients;
 
    double I_stim_;
  };
 
  // 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_cond_beta_gap_traub > recordablesMap_;
};
 
inline void
hh_cond_beta_gap_traub::update( Time const& origin, const long from, const long to )
{
  update_( origin, from, to, false );
}
 
inline bool
hh_cond_beta_gap_traub::wfr_update( Time const& origin, const long from, const long to )
{
  State_ old_state = S_;  // save state before wfr_update
  const bool wfr_tol_exceeded = update_( origin, from, to, true );
  S_ = old_state;  // restore old state
 
  return not wfr_tol_exceeded;
}
 
inline size_t
hh_cond_beta_gap_traub::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_cond_beta_gap_traub::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
hh_cond_beta_gap_traub::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
hh_cond_beta_gap_traub::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 size_t
hh_cond_beta_gap_traub::handles_test_event( GapJunctionEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline void
hh_cond_beta_gap_traub::get_status( Dictionary& d ) const
{
  P_.get( d );
  S_.get( d );
  ArchivingNode::get_status( d );
 
  d[ names::recordables ] = recordablesMap_.get_list();
 
  d[ names::t_spike ] = get_spiketime_ms();
}
 
inline void
hh_cond_beta_gap_traub::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;
 
  pre_run_hook();
}
 
}  // namespace
 
 
#endif  // HAVE_GSL
#endif  // HH_COND_BETA_GAP_TRAUB_H