An electrophysiological model of GABAergic double bouquet cells (Chrysanthidis et al. 2019)

 Download zip file 
Help downloading and running models
Accession:257610
We present an electrophysiological model of double bouquet cells (DBCs) and integrate them into an established cortical columnar microcircuit model that implements a BCPNN (Bayesian Confidence Propagation Neural Network) learning rule. The proposed architecture effectively solves the problem of duplexed learning of inhibition and excitation by replacing recurrent inhibition between pyramidal cells in functional columns of different stimulus selectivity with a plastic disynaptic pathway. The introduction of DBCs improves the biological plausibility of our model, without affecting the model's spiking activity, basic operation, and learning abilities.
Reference:
1 . Chrysanthidis N, Fiebig F, Lansner A (2019) Introducing double bouquet cells into a modular cortical associative memory model Journal of Computational Neuroscience
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Neocortex U1 interneuron basket PV GABA cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Abstract integrate-and-fire adaptive exponential (AdEx) neuron; Neocortex layer 2-3 interneuron; Neocortex bitufted interneuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEST;
Model Concept(s): Learning;
Implementer(s): Chrysanthidis, Nikolaos [nchr at kth.se]; Fiebig, Florian [fiebig at kth.se]; Lansner, Anders [ala at kth.se];
Search NeuronDB for information about:  Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex U1 interneuron basket PV GABA cell;
/
ChrysanthidisEtAl2019
BCPNN_NEST_Module
module-100725
autom4te.cache
libltdl
sli
aclocal.m4 *
aeif_cond_exp_multisynapse.cpp *
aeif_cond_exp_multisynapse.h *
bcpnn_connection.cpp *
bcpnn_connection.h *
bcpnn_connection_backup.cpp *
bcpnn_connection_backup.h *
bootstrap.sh *
compile *
config.guess *
config.sub *
configure *
configure.ac *
depcomp *
iaf_cond_alpha_bias.cpp *
iaf_cond_alpha_bias.h *
iaf_cond_exp_bias.cpp *
iaf_cond_exp_bias.h *
install-sh *
ltmain.sh *
Makefile.am *
Makefile.in *
missing *
pt_module.cpp *
pt_module.h *
pt_module_config.h.in *
pt_module_names.cpp *
pt_module_names.h *
                            
/*
 *  iaf_cond_exp_bias.cpp
 *
 *  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/>.
 *
 */

#include "iaf_cond_exp_bias.h"

#ifdef HAVE_GSL

#include "exceptions.h"
#include "network.h"
#include "dict.h"
#include "integerdatum.h"
#include "doubledatum.h"
#include "dictutils.h"
#include "numerics.h"
#include <limits>

#include "universal_data_logger_impl.h"
#include "event.h"

#include <iomanip>
#include <iostream>
using namespace std;
#include <cstdio>

/* ---------------------------------------------------------------- 
 * Recordables map
 * ---------------------------------------------------------------- */

nest::RecordablesMap<mynest::iaf_cond_exp_bias> mynest::iaf_cond_exp_bias::recordablesMap_;

namespace nest  // template specialization must be placed in namespace
{
  // Override the create() method with one call to RecordablesMap::insert_() 
  // for each quantity to be recorded.
  template <>
  void RecordablesMap<mynest::iaf_cond_exp_bias>::create()
  {
    // use standard names whereever you can for consistency!
    insert_(names::V_m, 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::V_M>);
    insert_(names::g_ex, 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::G_EXC>);
    insert_(names::g_in, 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::G_INH>);
    insert_(Name("z_j"), 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::Z_J>);
    insert_(Name("e_j"), 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::E_J>);
    insert_(Name("p_j"), 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::P_J>);
    insert_(Name("I_bias"), 
	    &mynest::iaf_cond_exp_bias::get_y_elem_<mynest::iaf_cond_exp_bias::State_::I_BIAS>);
    insert_(Name("epsilon"), 
	    &mynest::iaf_cond_exp_bias::get_epsilon_);
    insert_(Name("K"), 
	    &mynest::iaf_cond_exp_bias::get_K_);
  }
}

extern "C"
inline int mynest::iaf_cond_exp_bias_dynamics(double, const double y[], double f[], void* pnode)
{ 
  // a shorthand
  typedef mynest::iaf_cond_exp_bias::State_ S;

  // get access to node so we can almost work as in a member function
  assert(pnode);
  const mynest::iaf_cond_exp_bias& node =  *(reinterpret_cast<mynest::iaf_cond_exp_bias*>(pnode));

  // y[] here is---and must be---the state vector supplied by the integrator,
  // not the state vector in the node, node.S_.y[]. 

  // The following code is verbose for the sake of clarity. We assume that a
  // good compiler will optimize the verbosity away ...
  const double I_syn_exc = y[S::G_EXC] * (y[S::V_M] - node.P_.E_ex); 
  const double I_syn_inh = y[S::G_INH] * (y[S::V_M] - node.P_.E_in); 
  const double I_L       = node.P_.g_L * ( y[S::V_M] - node.P_.E_L );

  //V dot
  f[0]= ( - I_L + node.B_.I_stim_ + node.P_.I_e - I_syn_exc - I_syn_inh + y[S::I_BIAS]) / node.P_.C_m;

  f[1] = -y[S::G_EXC] / node.P_.tau_synE;
  f[2] = -y[S::G_INH] / node.P_.tau_synI;

  f[3] = (- y[S::Z_J] + node.P_.epsilon) / node.P_.tau_j;
  f[4] = (y[S::Z_J] - y[S::E_J]) / node.P_.tau_e;
  f[5] = node.P_.K * (y[S::E_J] - y[S::P_J]) / node.P_.tau_p;
  f[6] = node.P_.gain * std::log(y[S::P_J]);
  return GSL_SUCCESS;
}

/* ---------------------------------------------------------------- 
 * Default constructors defining default parameters and state
 * ---------------------------------------------------------------- */
    
mynest::iaf_cond_exp_bias::Parameters_::Parameters_()
  : V_th_      (-55.0    ),  // mV
    V_reset_   (-60.0    ),  // mV
    t_ref_     (  2.0    ),  // ms
    g_L        ( 16.6667 ),  // nS
    C_m        (250.0    ),  // pF
    E_ex       (  0.0    ),  // mV
    E_in       (-85.0    ),  // mV
    E_L        (-70.0    ),  // mV
    tau_synE   (  0.2    ),  // ms
    tau_synI   (  2.0    ),  // ms
    I_e        (  0.0    ),   // pA
    tau_j   ( 10.0    ),  // ms
    tau_e   (100.0    ),  // ms
    tau_p   (1000.0   ),  // ms
    K   (1.0      ), // dopamine
    fmax    (20.0     ), 
    gain    (1.0     ), 
    epsilon (0.001     ) 
{
}

mynest::iaf_cond_exp_bias::State_::State_(const Parameters_& p)
  : r_(0)
{
  y_[V_M] = p.E_L;
  y_[G_EXC] = y_[G_INH] = 0;
  y_[Z_J] = 0.01;
  y_[E_J] = 0.01;
  y_[P_J] = 0.01;
  y_[I_BIAS] = 0.0;
}

mynest::iaf_cond_exp_bias::State_::State_(const State_& s)
  : r_(s.r_)
{
  for ( size_t i = 0 ; i < STATE_VEC_SIZE ; ++i )
    y_[i] = s.y_[i];
}

mynest::iaf_cond_exp_bias::State_& mynest::iaf_cond_exp_bias::State_::operator=(const State_& s)
{
  assert(this != &s);  // would be bad logical error in program
  
  for ( size_t i = 0 ; i < STATE_VEC_SIZE ; ++i )
    y_[i] = s.y_[i];
  r_ = s.r_;
  return *this;
}

/* ---------------------------------------------------------------- 
 * Parameter and state extractions and manipulation functions
 * ---------------------------------------------------------------- */

void mynest::iaf_cond_exp_bias::Parameters_::get(DictionaryDatum &dd) const
{
  def<double>(dd,nest::names::V_th,         V_th_);
  def<double>(dd,nest::names::V_reset,      V_reset_);
  def<double>(dd,nest::names::t_ref,        t_ref_);
  def<double>(dd,nest::names::g_L,          g_L);
  def<double>(dd,nest::names::E_L,          E_L); 
  def<double>(dd,nest::names::E_ex,         E_ex);
  def<double>(dd,nest::names::E_in,         E_in);
  def<double>(dd,nest::names::C_m,          C_m);
  def<double>(dd,nest::names::tau_syn_ex,   tau_synE);
  def<double>(dd,nest::names::tau_syn_in,   tau_synI);
  def<double>(dd,nest::names::I_e,          I_e);
  def<nest::double_t>(dd, "tau_j",      tau_j);
  def<nest::double_t>(dd, "tau_e",      tau_e);
  def<nest::double_t>(dd, "tau_p",      tau_p);
  def<nest::double_t>(dd, "K",      K);
  def<nest::double_t>(dd, "gain",      gain);
  def<nest::double_t>(dd, "fmax",      fmax);
  def<nest::double_t>(dd, "epsilon",      epsilon);
}

void mynest::iaf_cond_exp_bias::Parameters_::set(const DictionaryDatum& dd)
{
  // allow setting the membrane potential
  updateValue<double>(dd,nest::names::V_th,    V_th_);
  updateValue<double>(dd,nest::names::V_reset, V_reset_);
  updateValue<double>(dd,nest::names::t_ref,   t_ref_);
  updateValue<double>(dd,nest::names::E_L,     E_L);
  
  updateValue<double>(dd,nest::names::E_ex,    E_ex);
  updateValue<double>(dd,nest::names::E_in,    E_in);
  
  updateValue<double>(dd,nest::names::C_m,     C_m);
  updateValue<double>(dd,nest::names::g_L,     g_L);

  updateValue<double>(dd,nest::names::tau_syn_ex, tau_synE);
  updateValue<double>(dd,nest::names::tau_syn_in, tau_synI);

  updateValue<double>(dd,nest::names::I_e,     I_e);
  updateValue<nest::double_t>(dd, "tau_j",      tau_j);
  updateValue<nest::double_t>(dd, "tau_e",      tau_e);
  updateValue<nest::double_t>(dd, "tau_p",      tau_p);
  updateValue<nest::double_t>(dd, "K",      K);
  updateValue<nest::double_t>(dd, "gain",      gain);
  updateValue<nest::double_t>(dd, "fmax",      fmax);
  updateValue<nest::double_t>(dd, "epsilon",      epsilon);

  if ( V_reset_ >= V_th_ )
    throw nest::BadProperty("Reset potential must be smaller than threshold.");
    
  if ( C_m <= 0 )
    throw nest::BadProperty("Capacitance must be strictly positive.");
    
  if ( t_ref_ < 0 )
    throw nest::BadProperty("Refractory time cannot be negative.");
      
  if ( tau_synE <= 0 || tau_synI <= 0 )
    throw nest::BadProperty("All time constants must be strictly positive.");
}

void mynest::iaf_cond_exp_bias::State_::get(DictionaryDatum &dd) const
{
  def<double>(dd, nest::names::V_m, y_[V_M]); // Membrane potential
  def<double>(dd, "I_bias", y_[I_BIAS]); 
}

void mynest::iaf_cond_exp_bias::State_::set(const DictionaryDatum& dd, const Parameters_&)
{
  updateValue<double>(dd, nest::names::V_m, y_[V_M]);
  updateValue<double>(dd, "I_bias", y_[I_BIAS]);
}

mynest::iaf_cond_exp_bias::Buffers_::Buffers_(iaf_cond_exp_bias& n)
  : logger_(n),
    s_(0),
    c_(0),
    e_(0)
{
  // Initialization of the remaining members is deferred to
  // init_buffers_().
}

mynest::iaf_cond_exp_bias::Buffers_::Buffers_(const Buffers_&, iaf_cond_exp_bias& n)
  : logger_(n),
    s_(0),
    c_(0),
    e_(0)
{
  // Initialization of the remaining members is deferred to
  // init_buffers_().
}

/* ---------------------------------------------------------------- 
 * Default and copy constructor for node, and destructor
 * ---------------------------------------------------------------- */

mynest::iaf_cond_exp_bias::iaf_cond_exp_bias()
  : Archiving_Node(), 
    P_(), 
    S_(P_),
    B_(*this)
{
  recordablesMap_.create();
}

mynest::iaf_cond_exp_bias::iaf_cond_exp_bias(const iaf_cond_exp_bias& n)
  : Archiving_Node(n), 
    P_(n.P_), 
    S_(n.S_),
    B_(n.B_, *this)
{
}

mynest::iaf_cond_exp_bias::~iaf_cond_exp_bias()
{
  // GSL structs may not have been allocated, so we need to protect destruction
  if ( B_.s_ ) gsl_odeiv_step_free(B_.s_);
  if ( B_.c_ ) gsl_odeiv_control_free(B_.c_);
  if ( B_.e_ ) gsl_odeiv_evolve_free(B_.e_);
}

/* ---------------------------------------------------------------- 
 * Node initialization functions
 * ---------------------------------------------------------------- */

void mynest::iaf_cond_exp_bias::init_state_(const Node& proto)
{
  const iaf_cond_exp_bias& pr = downcast<iaf_cond_exp_bias>(proto);
  S_ = pr.S_;
}

void mynest::iaf_cond_exp_bias::init_buffers_()
{
  B_.spike_exc_.clear();          // includes resize
  B_.spike_inh_.clear();          // includes resize
  B_.currents_.clear();           // includes resize
  Archiving_Node::clear_history();

  B_.logger_.reset();

  B_.step_ = nest::Time::get_resolution().get_ms();
  B_.IntegrationStep_ = B_.step_;

  static const gsl_odeiv_step_type* T1 = gsl_odeiv_step_rkf45;
  
  if ( B_.s_ == 0 )
    B_.s_ = gsl_odeiv_step_alloc (T1, State_::STATE_VEC_SIZE);
  else 
    gsl_odeiv_step_reset(B_.s_);
    
  if ( B_.c_ == 0 )  
    B_.c_ = gsl_odeiv_control_y_new (1e-3, 0.0);
  else
    gsl_odeiv_control_init(B_.c_, 1e-3, 0.0, 1.0, 0.0);
    
  if ( B_.e_ == 0 )  
    B_.e_ = gsl_odeiv_evolve_alloc(State_::STATE_VEC_SIZE);
  else 
    gsl_odeiv_evolve_reset(B_.e_);
  
  B_.sys_.function  = iaf_cond_exp_bias_dynamics; 
  B_.sys_.jacobian  = NULL;
  B_.sys_.dimension = State_::STATE_VEC_SIZE;
  B_.sys_.params    = reinterpret_cast<void*>(this);

  B_.I_stim_ = 0.0;
}

void mynest::iaf_cond_exp_bias::calibrate()
{
  B_.logger_.init();  // ensures initialization in case mm connected after Simulate

  V_.RefractoryCounts_ = nest::Time(nest::Time::ms(P_.t_ref_)).get_steps();
  assert(V_.RefractoryCounts_ >= 0);  // since t_ref_ >= 0, this can only fail in error
}

/* ---------------------------------------------------------------- 
 * Update and spike handling functions
 * ---------------------------------------------------------------- */

void mynest::iaf_cond_exp_bias::update(nest::Time const & origin, const nest::long_t from, const nest::long_t to)
{
   
  assert(to >= 0 && (nest::delay) from < nest::Scheduler::get_min_delay());
  assert(from < to);

  for ( nest::long_t lag = from ; lag < to ; ++lag )
  {
    
    double t = 0.0;

    // numerical integration with adaptive step size control:
    // ------------------------------------------------------
    // gsl_odeiv_evolve_apply performs only a single numerical
    // integration step, starting from t and bounded by step;
    // the while-loop ensures integration over the whole simulation
    // step (0, step] if more than one integration step is needed due
    // to a small integration step size;
    // note that (t+IntegrationStep > step) leads to integration over
    // (t, step] and afterwards setting t to step, but it does not
    // enforce setting IntegrationStep to step-t; this is of advantage
    // for a consistent and efficient integration across subsequent
    // simulation intervals
    while ( t < B_.step_ )
    {
      const int status = gsl_odeiv_evolve_apply(B_.e_, B_.c_, B_.s_, 
			   &B_.sys_,             // system of ODE
			   &t,                   // from t
			    B_.step_,            // to t <= step
			   &B_.IntegrationStep_, // integration step size
			    S_.y_); 	         // neuronal state

      if ( status != GSL_SUCCESS )
        throw nest::GSLSolverFailure(get_name(), status);
    }

    S_.y_[State_::I_BIAS] = P_.gain * std::log(S_.y_[State_::P_J]);

    S_.y_[State_::G_EXC] += B_.spike_exc_.get_value(lag);
    S_.y_[State_::G_INH] += B_.spike_inh_.get_value(lag);

    // absolute refractory period
    if ( S_.r_ )
    {// neuron is absolute refractory
      --S_.r_; 
      S_.y_[State_::V_M] = P_.V_reset_; 
    }
    else
      // neuron is not absolute refractory
      if ( S_.y_[State_::V_M] >= P_.V_th_ )
	    {
	      S_.r_              = V_.RefractoryCounts_;
	      S_.y_[State_::V_M] = P_.V_reset_;

              S_.y_[State_::Z_J] += (1000.0/(P_.fmax*B_.step_) - S_.y_[State_::Z_J] + P_.epsilon) * B_.step_ / P_.tau_j;
              S_.y_[State_::E_J] += (S_.y_[State_::Z_J] - S_.y_[State_::E_J]) * B_.step_ / P_.tau_e;
              S_.y_[State_::P_J] += P_.K * (S_.y_[State_::E_J] - S_.y_[State_::P_J]) * B_.step_ / P_.tau_p;

	      set_spiketime(nest::Time::step(origin.get_steps()+lag+1));
	  
	      nest::SpikeEvent se;
	      network()->send(*this, se, lag);
	    }

    S_.y_[State_::I_BIAS] = P_.gain * std::log(S_.y_[State_::P_J]);
    
    // set new input current
    B_.I_stim_ = B_.currents_.get_value(lag);

    // log state data
    B_.logger_.record_data(origin.get_steps() + lag);
  }
}

void mynest::iaf_cond_exp_bias::handle(nest::SpikeEvent & e)
{
  assert(e.get_delay() > 0);

  if(e.get_weight() > 0.0)
    B_.spike_exc_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
			    e.get_weight() * e.get_multiplicity() );
  else
    B_.spike_inh_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
			    -e.get_weight() * e.get_multiplicity() );  // ensure conductance is positive
}

void mynest::iaf_cond_exp_bias::handle(nest::CurrentEvent& e)
{
  assert(e.get_delay() > 0);

  const double_t c=e.get_current();
  const double_t w=e.get_weight();

  // add weighted current; HEP 2002-10-04
  B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()), 
		      w *c);
}

void mynest::iaf_cond_exp_bias::handle(nest::DataLoggingRequest& e)
{
  B_.logger_.handle(e);
}

#endif //HAVE_GSL