iaf_cond_alpha_mc.h

Go to the documentation of this file.

/*
 *  iaf_cond_alpha_mc.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 IAF_COND_ALPHA_MC_H
#define IAF_COND_ALPHA_MC_H

#include "config.h"

#ifdef HAVE_GSL

#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 "dictdatum.h"
#include "name.h"

#include <vector>

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

/* BeginDocumentation
Name: iaf_cond_alpha_mc - PROTOTYPE Multi-compartment conductance-based leaky integrate-and-fire neuron model.

Description:
THIS MODEL IS A PROTOTYPE FOR ILLUSTRATION PURPOSES. IT IS NOT YET
FULLY TESTED. USE AT YOUR OWN PERIL!

iaf_cond_alpha_mc is an implementation of a multi-compartment spiking
neuron using IAF dynamics with conductance-based synapses. It serves
mainly to illustrate the implementation of multicompartment models in
NEST.

The model has three compartments: soma, proximal and distal dendrite,
labeled as s, p, and d, respectively. Compartments are connected through 
passive conductances as follows

C_m.s d/dt V_m.s = ... - g_sp ( V_m.s - V_m.p )

C_m.p d/dt V_m.p = ... - g_sp ( V_m.p - V_m.s ) - g_pd ( V_m.p - V_m.d )

C_m.d d/dt V_m.d = ...                          - g_pd ( V_m.d - V_m.p )

A spike is fired when the somatic membrane potential exceeds threshold,
V_m.s >= V_th. After a spike, somatic membrane potential is clamped to
a reset potential, V_m.s == V_reset, for the refractory period. Dendritic
membrane potentials are not manipulated after a spike.

There is one excitatory and one inhibitory conductance-based synapse
onto each compartment, with alpha-function time course. The alpha
function is normalised such that an event of weight 1.0 results in a
peak current of 1 nS at t = tau_syn. Each compartment can also receive
current input from a current generator, and an external (rheobase)
current can be set for each compartment.

Synapses, including those for injection external currents, are addressed through 
the receptor types given in the receptor_types entry of the state dictionary. Note
that in contrast to the single-compartment iaf_cond_alpha model, all synaptic
weights must be positive numbers!


Parameters: 
The following parameters can be set in the status dictionary. Parameters 
for each compartment are collected in a sub-dictionary; these sub-dictionaries
are called "soma", "proximal", and "distal", respectively. In the list below,
these parameters are marked with an asterisk.

V_m*         double - Membrane potential in mV 
E_L*         double - Leak reversal potential in mV.
C_m*         double - Capacity of the membrane in pF
E_ex*        double - Excitatory reversal potential in mV.
E_in*        double - Inhibitory reversal potential in mV.
g_L*         double - Leak conductance in nS;
tau_syn_ex*  double - Rise time of the excitatory synaptic alpha function in ms.
tau_syn_in*  double - Rise time of the inhibitory synaptic alpha function in ms.
I_e*         double - Constant input current in pA.

g_sp         double - Conductance connecting soma and proximal dendrite, in nS.
g_pd         double - Conductance connecting proximal and distal dendrite, in nS.
t_ref        double - Duration of refractory period in ms. 
V_th         double - Spike threshold in mV.
V_reset      double - Reset potential of the membrane in mV.

Example:
See examples/nest/mc_neuron.py.

Remark: 
This is a prototype for illustration which has undergone only limited testing. 
Details of the implementation and user-interface will likely change.
USE AT YOUR OWN PERIL!

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

References: 

Meffin, H., Burkitt, A. N., & Grayden, D. B. (2004). An analytical
model for the large, fluctuating synaptic conductance state typical of
neocortical neurons in vivo. J.  Comput. Neurosci., 16, 159–175.

Bernander, O ., Douglas, R. J., Martin, K. A. C., & Koch, C. (1991).
Synaptic background activity influences spatiotemporal integration in
single pyramidal cells.  Proc. Natl. Acad. Sci. USA, 88(24),
11569–11573.

Author: Plesser

SeeAlso: iaf_cond_alpha
*/

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

  class iaf_cond_alpha_mc : public Archiving_Node
  {

    // Boilerplate function declarations --------------------------------
    
  public:
    
    iaf_cond_alpha_mc();
    iaf_cond_alpha_mc(const iaf_cond_alpha_mc&);
    ~iaf_cond_alpha_mc();

    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(Time const &, const long_t, const long_t);

    // Enumerations and constants specifying structure and properties ----

    enum Compartments_ { SOMA = 0, PROX, DIST, NCOMP };

    static const port MIN_SPIKE_RECEPTOR = 1;

    enum SpikeSynapseTypes { SOMA_EXC=MIN_SPIKE_RECEPTOR, SOMA_INH, 
                 PROX_EXC  , PROX_INH, 
                 DIST_EXC  , DIST_INH, 
                 SUP_SPIKE_RECEPTOR };

    static const size_t NUM_SPIKE_RECEPTORS = SUP_SPIKE_RECEPTOR - MIN_SPIKE_RECEPTOR;

    static const port MIN_CURR_RECEPTOR = SUP_SPIKE_RECEPTOR;

    enum CurrentSynapseTypes { I_SOMA=MIN_CURR_RECEPTOR, I_PROX, I_DIST, 
                   SUP_CURR_RECEPTOR };

    static const size_t NUM_CURR_RECEPTORS = SUP_CURR_RECEPTOR - MIN_CURR_RECEPTOR;

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

    friend int iaf_cond_alpha_mc_dynamics (double, const double*, double*, void*);

    friend class RecordablesMap<iaf_cond_alpha_mc>;
    friend class UniversalDataLogger<iaf_cond_alpha_mc>;


    // Parameters ------------------------------------------------------ 

    struct Parameters_ {
      double_t V_th;       
      double_t V_reset;    
      double_t t_ref;      
      double_t g_conn[NCOMP-1];    
      double_t g_L[NCOMP];         
      double_t C_m[NCOMP];         
      double_t E_ex[NCOMP];        
      double_t E_in[NCOMP];        
      double_t E_L[NCOMP];         
      double_t tau_synE[NCOMP];    
      double_t tau_synI[NCOMP];    
      double_t I_e[NCOMP];         
  
      Parameters_();  
      Parameters_(const Parameters_&); 
      Parameters_& operator=(const Parameters_&); 

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

    // State variables  ------------------------------------------------------ 

  public:
    struct State_ {
      
      enum StateVecElems_ { V_M = 0,
                DG_EXC, G_EXC,
                DG_INH, G_INH,
                STATE_VEC_COMPS };

      static const size_t STATE_VEC_SIZE = STATE_VEC_COMPS * NCOMP;

      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_&);

      static size_t idx(size_t comp, StateVecElems_ elem)
      { return comp * STATE_VEC_COMPS + elem; }

    };    
  private:

    // Internal buffers -------------------------------------------------------- 

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

      UniversalDataLogger<iaf_cond_alpha_mc> logger_;

      std::vector<RingBuffer> spikes_;
      std::vector<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_[NCOMP];      
    };

    // Internal variables --------------------------------------------- 
    
    struct Variables_ { 
      double_t PSConInit_E_[NCOMP]; 
      
      double_t PSConInit_I_[NCOMP];    
      
      int_t    RefractoryCounts_;
    };

    // Access functions for UniversalDataLogger -------------------------------

    template <State_::StateVecElems_ elem, Compartments_ comp>
    double_t get_y_elem_() const { return S_.y_[S_.idx(comp, elem)]; }
    
    double_t get_r_() const { return Time::get_resolution().get_ms() * S_.r_; }
    
    // Data members ---------------------------------------------------- 

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

    static std::vector<Name> comp_names_;

    // static DictionaryDatum receptor_dict_;

    static RecordablesMap<iaf_cond_alpha_mc> recordablesMap_;
  };
  
  inline
  port iaf_cond_alpha_mc::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 iaf_cond_alpha_mc::handles_test_event(SpikeEvent&, rport receptor_type)
  {
    if ( receptor_type < MIN_SPIKE_RECEPTOR || receptor_type >= SUP_SPIKE_RECEPTOR )
    {
      if ( receptor_type < 0 || receptor_type >= SUP_CURR_RECEPTOR )
    throw UnknownReceptorType(receptor_type, get_name());
      else
    throw IncompatibleReceptorType(receptor_type, get_name(), "SpikeEvent");
    }
    return receptor_type - MIN_SPIKE_RECEPTOR;
  }
 
  inline
  port iaf_cond_alpha_mc::handles_test_event(CurrentEvent&, rport receptor_type)
  {
    if ( receptor_type < MIN_CURR_RECEPTOR || receptor_type >= SUP_CURR_RECEPTOR )
    {
      if ( receptor_type >= 0 && receptor_type < MIN_CURR_RECEPTOR ) 
    throw IncompatibleReceptorType(receptor_type, get_name(), "CurrentEvent");
      else
    throw UnknownReceptorType(receptor_type, get_name());
    }
    return receptor_type - MIN_CURR_RECEPTOR;
  }
 
  inline
  port iaf_cond_alpha_mc::handles_test_event(DataLoggingRequest& dlr, 
                         rport receptor_type)
  {
    if ( receptor_type != 0 )
    {
      if ( receptor_type < 0 || receptor_type >= SUP_CURR_RECEPTOR )
    throw UnknownReceptorType(receptor_type, get_name());
      else
    throw IncompatibleReceptorType(receptor_type, get_name(), "DataLoggingRequest");
    }
    return B_.logger_.connect_logging_device(dlr, recordablesMap_);
  }

  inline
  void iaf_cond_alpha_mc::get_status(DictionaryDatum &d) const
  {
    P_.get(d);
    S_.get(d);
    Archiving_Node::get_status(d);
    
    (*d)[names::recordables] = recordablesMap_.get_list();

    DictionaryDatum receptor_dict_ = new Dictionary();
    (*receptor_dict_)[Name("soma_exc")]  = SOMA_EXC;
    (*receptor_dict_)[Name("soma_inh")]  = SOMA_INH;
    (*receptor_dict_)[Name("soma_curr")] = I_SOMA;
    
    (*receptor_dict_)[Name("proximal_exc")]  = PROX_EXC;
    (*receptor_dict_)[Name("proximal_inh")]  = PROX_INH;
    (*receptor_dict_)[Name("proximal_curr")] = I_PROX;
    
    (*receptor_dict_)[Name("distal_exc")]  = DIST_EXC;
    (*receptor_dict_)[Name("distal_inh")]  = DIST_INH;
    (*receptor_dict_)[Name("distal_curr")] = I_DIST;

    (*d)[names::receptor_types] = receptor_dict_;
  }

  inline
  void iaf_cond_alpha_mc::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
#endif //IAF_COND_ALPHA_MC_H