aeif_cond_exp.h

Go to the documentation of this file.

/*
 *  aeif_cond_exp.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 AEIF_COND_EXP_H
#define AEIF_COND_EXP_H

#include "config.h"

#ifdef HAVE_GSL_1_11

#include "nest.h"
#include "event.h"
#include "archiving_node.h"
#include "ring_buffer.h"
#include "connection.h"
#include "universal_data_logger.h"
#include "recordables_map.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

/* BeginDocumentation
Name: aeif_cond_exp - Conductance based exponential integrate-and-fire neuron model according to Brette and Gerstner (2005).

Description: 

aeif_cond_exp is the adaptive exponential integrate and fire neuron
according to Brette and Gerstner (2005), with post-synaptic
conductances in the form of truncated exponentials.

This implementation uses the embedded 4th order Runge-Kutta-Fehlberg
solver with adaptive stepsize to integrate the differential equation.

The membrane potential is given by the following differential equation:
C dV/dt= -g_L(V-E_L)+g_L*Delta_T*exp((V-V_T)/Delta_T)-g_e(t)(V-E_e) -g_i(t)(V-E_i)-w +I_e

and

tau_w * dw/dt= a(V-E_L) -W


Note that the spike detection threshold V_peak is automatically set to
V_th+10 mV to avoid numerical instabilites that may result from
setting V_peak too high.

Parameters: 
The following parameters can be set in the status dictionary.

Dynamic state variables:
  V_m        double - Membrane potential in mV
  g_ex       double - Excitatory synaptic conductance in nS.
  g_in       double - Inhibitory synaptic conductance in nS.
  w          double - Spike-adaptation current in pA.

Membrane Parameters:
  C_m        double - Capacity of the membrane in pF
  t_ref      double - Duration of refractory period in ms. 
  V_reset    double - Reset value for V_m after a spike. In mV.
  E_L        double - Leak reversal potential in mV. 
  g_L        double - Leak conductance in nS.
  I_e        double - Constant external input current in pA.

Spike adaptation parameters:
  a          double - Subthreshold adaptation in nS.
  b          double - Spike-triggered adaptation in pA.
  Delta_T    double - Slope factor in mV
  tau_w      double - Adaptation time constant in ms
  V_t        double - Spike initiation threshold in mV
  V_peak     double - Spike detection threshold in mV.

Synaptic parameters
  E_ex       double - Excitatory reversal potential in mV.
  tau_syn_ex double - Rise time of excitatory synaptic conductance in ms (exp function).
  E_in       double - Inhibitory reversal potential in mV.
  tau_syn_in double - Rise time of the inhibitory synaptic conductance in ms (exp function).

Integration parameters
  gsl_error_tol  double - This parameter controls the admissible error of the GSL integrator.
                          Reduce it if NEST complains about numerical instabilities.

Author: Adapted from aeif_cond_alpha by Lyle Muller

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

References: Brette R and Gerstner W (2005) Adaptive Exponential Integrate-and-Fire Model as 
            an Effective Description of Neuronal Activity. J Neurophysiol 94:3637-3642

SeeAlso: iaf_cond_exp, aeif_cond_alpha
*/

namespace nest
{
  extern "C"
  int aeif_cond_exp_dynamics (double, const double*, double*, void*);

  class aeif_cond_exp:
    public Archiving_Node
  {
    
  public:        
    
    aeif_cond_exp();
    aeif_cond_exp(const aeif_cond_exp&);
    ~aeif_cond_exp();

    using Node::handle;
    using Node::handles_test_event;

    port send_test_event(Node&, rport, synindex, bool);
    
    void handle(SpikeEvent &);
    void handle(CurrentEvent &);
    void handle(DataLoggingRequest &);
    
    port handles_test_event(SpikeEvent&, rport);
    port handles_test_event(CurrentEvent&, rport);
    port handles_test_event(DataLoggingRequest&, rport);

    void get_status(DictionaryDatum &) const;
    void set_status(const DictionaryDatum &);
    
  private:
    
    void init_state_(const Node &proto);
    void init_buffers_();
    void calibrate();
    void update(const Time &, const long_t, const long_t);

    // END Boilerplate function declarations ----------------------------

    // Friends --------------------------------------------------------

    // make dynamics function quasi-member
    friend int aeif_cond_exp_dynamics(double, const double*, double*, void*);

    // The next two classes need to be friends to access the State_ class/member
    friend class RecordablesMap<aeif_cond_exp>;
    friend class UniversalDataLogger<aeif_cond_exp>;

  private:
    // ---------------------------------------------------------------- 

    struct Parameters_
    {
      double_t V_peak_;     
      double_t V_reset_;    
      double_t t_ref_;      

      double_t g_L;         
      double_t C_m;         
      double_t E_ex;        
      double_t E_in;        
      double_t E_L;         
      double_t Delta_T;     
      double_t tau_w;       
      double_t a;           
      double_t b;           
      double_t V_th;        
      double_t t_ref;       
      double_t tau_syn_ex;  
      double_t tau_syn_in;  
      double_t I_e;         
  
      double_t gsl_error_tol;   
  
      Parameters_();  

      void get(DictionaryDatum &) const;  
      void set(const DictionaryDatum &);  
    };

  public:
    // ---------------------------------------------------------------- 

    struct State_
    {
      enum StateVecElems
      {
    V_M   = 0,
    G_EXC    ,  // 1
    G_INH    ,  // 2
    W        ,  // 3
    STATE_VEC_SIZE
      };

      double_t y_[STATE_VEC_SIZE];  
      int_t    r_;           

      State_(const Parameters_ &);  
      State_(const State_ &);
      State_& operator = (const State_ &);

      void get(DictionaryDatum &) const;
      void set(const DictionaryDatum &, const Parameters_ &);
    };

    // ---------------------------------------------------------------- 

    struct Buffers_
    {
      Buffers_(aeif_cond_exp &);                    
      Buffers_(const Buffers_ &, aeif_cond_exp &);  

      UniversalDataLogger<aeif_cond_exp> 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_;  
      
      // IntergrationStep_ 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_t step_;             
      double   IntegrationStep_;  

      double_t I_stim_;
    };

    // ---------------------------------------------------------------- 

    struct Variables_
    {
      int_t RefractoryCounts_;
    };

    // Access functions for UniversalDataLogger -------------------------------
    
    template <State_::StateVecElems elem>
    double_t get_y_elem_() const { return S_.y_[elem]; }

    // ---------------------------------------------------------------- 

    Parameters_ P_;
    State_      S_;
    Variables_  V_;
    Buffers_    B_;

    static RecordablesMap<aeif_cond_exp> recordablesMap_;
  };

  inline
  port aeif_cond_exp::send_test_event(Node& target, rport receptor_type, synindex, bool)
  {
    SpikeEvent e;
    e.set_sender(*this);

    return target.handles_test_event(e, receptor_type);
  }

  inline
  port aeif_cond_exp::handles_test_event(SpikeEvent&, rport receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return 0;
  }

  inline
  port aeif_cond_exp::handles_test_event(CurrentEvent&, rport receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return 0;
  }

  inline
  port aeif_cond_exp::handles_test_event(DataLoggingRequest &dlr,
                    rport receptor_type)
  {
    if (receptor_type != 0)
      throw UnknownReceptorType(receptor_type, get_name());
    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
  }
 
  inline
  void aeif_cond_exp::get_status(DictionaryDatum &d) const
  {
    P_.get(d);
    S_.get(d);
    Archiving_Node::get_status(d);

    (*d)[names::recordables] = recordablesMap_.get_list();
  }

  inline
  void aeif_cond_exp::set_status(const DictionaryDatum &d)
  {
    Parameters_ ptmp = P_;  // temporary copy in case of errors
    ptmp.set(d);            // throws if BadProperty
    State_      stmp = S_;  // temporary copy in case of errors
    stmp.set(d, ptmp);      // 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.
    Archiving_Node::set_status(d);

    // if we get here, temporaries contain consistent set of properties
    P_ = ptmp;
    S_ = stmp;
  }
  
} // namespace

#endif // HAVE_GSL_1_11
#endif // AEIF_COND_EXP_H