ht_neuron.h

Go to the documentation of this file.

/*
 *  ht_neuron.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 HT_NEURON_H
#define HT_NEURON_H
 
// Generated includes:
#include "config.h"
 
#ifdef HAVE_GSL
 
// C++ includes:
#include <string>
#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 "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"
 
 
namespace nest
{
extern "C" int ht_neuron_dynamics( double, const double*, double*, void* );
 
/* BeginUserDocs: neuron, adaptation, integrate-and-fire, soft threshold, Hill-Tononi
 
Short description
+++++++++++++++++
 
Neuron model after Hill & Tononi (2005)
 
Description
+++++++++++
 
This model neuron implements a slightly modified version of the
neuron model described in :footcite:p:`Hill2005`. The most important properties are:
 
- Integrate-and-fire with adaptive threshold.
- Repolarizing potassium current instead of hard reset.
- AMPA, NMDA, GABA_A, and GABA_B conductance-based synapses with
  beta-function (difference of exponentials) time course.
- Voltage-dependent NMDA with instantaneous or two-stage unblocking :footcite:p:`Hill2005`, :footcite:p:`Vargas2003`.
- Intrinsic currents I_h, I_T, I_Na(p), and I_KNa.
- Synaptic "minis" are not implemented.
 
For implementation details see:
 
- `HillTononi_model <../model_details/HillTononiModels.ipynb>`_
 
For examples, see:
 
- :doc:`../auto_examples/intrinsic_currents_spiking`
- :doc:`../auto_examples/intrinsic_currents_subthreshold`
 
For an example network model using ``ht_neuron`` (based on :footcite:p:`Hill2005`), see:
 
- `Multiarea Hill-Tononi thalamocortical network model
  <https://github.com/ricardomurphy/Multiarea-Hill-Tononi-thalamocortical-network-model>`_
 
Parameters
++++++++++
 
=============== ======= =========================================================
 V_m            mV      Membrane potential
 tau_m          ms      Membrane time constant applying to all currents except
                        repolarizing K-current (see :footcite:p:`Hill2005`, p 1677)
 t_ref          ms      Refractory time and duration of post-spike repolarizing
                        potassium current (t_spike in :footcite:p:`Hill2005`)
 tau_spike      ms      Membrane time constant for post-spike repolarizing
                        potassium current
 voltage_clamp  boolean If true, clamp voltage to value at beginning of
 simulation
                        (default: false, mainly for testing)
 theta          mV      Threshold
 theta_eq       mV      Equilibrium value
 tau_theta      ms      Time constant
 g_KL           nS      Conductance for potassium leak current
 E_K            mV      Reversal potential for potassium leak currents
 g_NaL          nS      Conductance for sodium leak currents
 E_Na           mV      Reversal potential for Na leak currents
 tau_D_KNa      ms      Relaxation time constant for I_KNa
 receptor_types         Dictionary mapping synapse names to ports on neuron model
 recordables            List of recordable quantities
=============== ======= =========================================================
 
+------------------------------------------------------------+
|{E_rev,g_peak,tau_rise,tau_decay}_{AMPA,NMDA,GABA_A,GABA_B} |
+------------------------------------------------------------+
| Reversal potentials, peak conductances and time constants  |
| for synapses (tau_rise/tau_decay correspond to tau_1/tau_2 |
| in the paper)                                              |
+------------------------------------------------------------+
 
+------------------------+------------------------------------------------+
|V_act_NMDA, S_act_NMDA, |  Parameters for voltage dependence of NMDA-    |
|tau_Mg_{fast, slow}_NMDA|  conductance, see above                        |
+------------------------+------------------------------------------------+
 
============================ =================================================
instant_unblock_NMDA         Instantaneous NMDA unblocking (default: false)
{E_rev,g_peak}_{h,T,NaP,KNa} Reversal potential and peak conductance for
                             intrinsic currents
{N}_{T,NaP}                  Exponent activation term m_inf, corresponding to
                             number of activation particles
equilibrate                  If given and true, time-dependent activation
                             and inactivation state variables (h, m) of
                             intrinsic currents and NMDA channels are set
                             to their equilibrium values during this
                             SetStatus call; otherwise they retain their
                             present values.
============================ =================================================
 
.. note::
   Conductances are unitless in this model and currents are in mV.
 
Sends
+++++
 
SpikeEvent
 
Receives
++++++++
 
SpikeEvent, CurrentEvent, DataLoggingRequest
 
References
++++++++++
 
.. footbibliography::
 
See also
++++++++
 
ht_synapse
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: ht_neuron
 
EndUserDocs */
 
void register_ht_neuron( const std::string& name );
 
class ht_neuron : public ArchivingNode
{
public:
  ht_neuron();
  ht_neuron( const ht_neuron& );
  ~ht_neuron() override;
 
  using Node::handle;
  using Node::handles_test_event;
 
  size_t send_test_event( Node&, size_t, synindex, bool ) override;
 
  void handle( SpikeEvent& e ) override;
  void handle( CurrentEvent& e ) 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:
  enum SynapseTypes
  {
    INF_SPIKE_RECEPTOR = 0,
    AMPA,
    NMDA,
    GABA_A,
    GABA_B,
    SUP_SPIKE_RECEPTOR
  };
 
  void init_buffers_() override;
  void pre_run_hook() override;
 
  void update( Time const&, const long, const long ) override;
 
  double get_synapse_constant( double, double, double );
 
  // END Boilerplate function declarations ----------------------------
 
  // Friends --------------------------------------------------------
 
  // make dynamics function quasi-member
  friend int ht_neuron_dynamics( double, const double*, double*, void* );
 
  // ----------------------------------------------------------------
 
  struct Parameters_
  {
    Parameters_();
 
    void get( Dictionary& ) const;              
    void set( const Dictionary&, Node* node );  
 
    // Note: Conductances are unitless
    // Leaks
    double E_Na;  // mV
    double E_K;   // mV
    double g_NaL;
    double g_KL;
    double tau_m;  // ms
 
    // Dynamic threshold
    double theta_eq;   // mV
    double tau_theta;  // ms
 
    // Post-spike potassium current
    double tau_spike;  // ms, membrane time constant for this current
    double t_ref;      // ms, refractory time
 
    // Parameters for synapse of type AMPA, GABA_A, GABA_B and NMDA
    double g_peak_AMPA;
    double tau_rise_AMPA;   // ms
    double tau_decay_AMPA;  // ms
    double E_rev_AMPA;      // mV
 
    double g_peak_NMDA;
    double tau_rise_NMDA;     // ms
    double tau_decay_NMDA;    // ms
    double E_rev_NMDA;        // mV
    double V_act_NMDA;        // mV, inactive for V << Vact, inflection of sigmoid
    double S_act_NMDA;        // mV, scale of inactivation
    double tau_Mg_slow_NMDA;  // ms
    double tau_Mg_fast_NMDA;  // ms
    bool instant_unblock_NMDA;
 
    double g_peak_GABA_A;
    double tau_rise_GABA_A;   // ms
    double tau_decay_GABA_A;  // ms
    double E_rev_GABA_A;      // mV
 
    double g_peak_GABA_B;
    double tau_rise_GABA_B;   // ms
    double tau_decay_GABA_B;  // ms
    double E_rev_GABA_B;      // mV
 
    // parameters for intrinsic currents
    double g_peak_NaP;
    double E_rev_NaP;  // mV
    double N_NaP;
 
    double g_peak_KNa;
    double E_rev_KNa;  // mV
    double tau_D_KNa;  // ms
 
    double g_peak_T;
    double E_rev_T;  // mV
    double N_T;
 
    double g_peak_h;
    double E_rev_h;  // mV
 
    bool voltage_clamp;
  };
 
  // ----------------------------------------------------------------
 
public:
  struct State_
  {
    // y_ = [V, theta, Synapses]
    enum StateVecElems_
    {
      V_M = 0,
      THETA,
      DG_AMPA,
      G_AMPA,
      DG_NMDA_TIMECOURSE,
      G_NMDA_TIMECOURSE,
      DG_GABA_A,
      G_GABA_A,
      DG_GABA_B,
      G_GABA_B,  // DO NOT INSERT ANYTHING UP TO HERE, WILL MIX UP
                 // SPIKE DELIVERY
      m_fast_NMDA,
      m_slow_NMDA,
      m_Ih,
      D_IKNa,
      m_IT,
      h_IT,
      STATE_VEC_SIZE
    };
 
    double y_[ STATE_VEC_SIZE ];
 
    long ref_steps_;
 
    double I_NaP_;  
    double I_KNa_;  
    double I_T_;    
    double I_h_;    
 
    State_( const ht_neuron&, const Parameters_& p );
    State_( const State_& s );
 
    State_& operator=( const State_& s );
 
    ~State_();
 
    void get( Dictionary& ) const;
    void set( const Dictionary&, const ht_neuron&, Node* node );
  };
 
private:
  // These friend declarations must be precisely here.
  friend class RecordablesMap< ht_neuron >;
  friend class UniversalDataLogger< ht_neuron >;
 
 
  // ----------------------------------------------------------------
 
  struct Buffers_
  {
    Buffers_( ht_neuron& );
    Buffers_( const Buffers_&, ht_neuron& );
 
    UniversalDataLogger< ht_neuron > logger_;
 
    std::vector< RingBuffer > spike_inputs_;
    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 integration_step_;  
 
    double I_stim_;
  };
 
  // ----------------------------------------------------------------
 
  struct Variables_
  {
    std::vector< double > cond_steps_;
 
    int PotassiumRefractoryCounts_;
 
    double V_clamp_;
  };
 
 
  // readout functions, can use template for vector elements
  template < State_::StateVecElems_ elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }
  double
  get_I_NaP_() const
  {
    return S_.I_NaP_;
  }
  double
  get_I_KNa_() const
  {
    return S_.I_KNa_;
  }
  double
  get_I_T_() const
  {
    return S_.I_T_;
  }
  double
  get_I_h_() const
  {
    return S_.I_h_;
  }
 
  double get_g_NMDA_() const;
 
  double m_NMDA_( double V, double m_eq, double m_fast, double m_slow ) const;
 
  double m_eq_h_( double V ) const;
 
  double m_eq_T_( double V ) const;
 
  double h_eq_T_( double V ) const;
 
  double m_eq_NMDA_( double V ) const;
 
  double D_eq_KNa_( double V ) const;
 
  static RecordablesMap< ht_neuron > recordablesMap_;
 
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
};
 
 
inline size_t
ht_neuron::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
ht_neuron::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  assert( B_.spike_inputs_.size() == 4 );
 
  if ( not( INF_SPIKE_RECEPTOR < receptor_type and receptor_type < SUP_SPIKE_RECEPTOR ) )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
    return 0;
  }
  else
  {
    return receptor_type - 1;
  }
}
 
inline size_t
ht_neuron::handles_test_event( CurrentEvent&, size_t receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}
 
inline size_t
ht_neuron::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_ );
}
}
 
#endif  // HAVE_GSL
#endif  // HT_NEURON_H