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
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
bootstrap-module-100725
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_alpha_bias.cpp
 *
 *  This file is part of NEST
 *
 *  Copyright (C) 2005-2009 by
 *  The NEST Initiative
 *
 *  See the file AUTHORS for details.
 *
 *  Permission is granted to compile and modify
 *  this file for non-commercial use.
 *  See the file LICENSE for details.
 *
 */


#include "iaf_cond_alpha_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 "universal_data_logger_impl.h"
#include <limits>

#include <iomanip>
#include <iostream>
#include <cstdio>

#include <iostream>
using namespace std;

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

nest::RecordablesMap<mynest::iaf_cond_alpha_bias> mynest::iaf_cond_alpha_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_alpha_bias>::create()
  {
    // use standard names whereever you can for consistency!
    insert_(names::V_m, 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::V_M>);
    insert_(names::g_ex, 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::G_EXC>);
    insert_(names::g_in, 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::G_INH>);
    insert_(Name("z_j"), 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::Z_J>);
    insert_(Name("e_j"), 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::E_J>);
    insert_(Name("p_j"), 
	    &mynest::iaf_cond_alpha_bias::get_y_elem_<mynest::iaf_cond_alpha_bias::State_::P_J>);

    insert_(names::t_ref_remaining, 
	    &mynest::iaf_cond_alpha_bias::get_r_);

    insert_(Name("bias"), 
	    &mynest::iaf_cond_alpha_bias::get_bias_);

    insert_(Name("epsilon"), 
	    &mynest::iaf_cond_alpha_bias::get_epsilon_);

    insert_(Name("kappa"), 
	    &mynest::iaf_cond_alpha_bias::get_kappa_);
  }
}

/* ---------------------------------------------------------------- 
 * Iteration function
 * ---------------------------------------------------------------- */

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

  // get access to node so we can almost work as in a member function
  assert(pnode);
  const mynest::iaf_cond_alpha_bias& node =  *(reinterpret_cast<mynest::iaf_cond_alpha_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 nest::double_t I_syn_exc = y[S::G_EXC] * ( y[S::V_M] - node.P_.E_ex );
  const nest::double_t I_syn_inh = y[S::G_INH] * ( y[S::V_M] - node.P_.E_in );
  const nest::double_t I_leak    = node.P_.g_L * ( y[S::V_M] - node.P_.E_L  );
  const nest::double_t I_bias    = node.P_.gain * std::log(y[S::P_J]);
  
  // dV_m/dt
  f[0] = (-I_leak - I_syn_exc - I_syn_inh + node.B_.I_stim_ + node.P_.I_e + I_bias) / node.P_.C_m;

  // d dg_exc/dt, dg_exc/dt
  f[1] = -y[S::DG_EXC] / node.P_.tau_synE;
  f[2] =  y[S::DG_EXC] - (y[S::G_EXC]/node.P_.tau_synE); 

  // d dg_exc/dt, dg_exc/dt
  f[3] = -y[S::DG_INH] / node.P_.tau_synI;
  f[4] =  y[S::DG_INH] - (y[S::G_INH]/node.P_.tau_synI); 

  f[5] = (- y[S::Z_J] + node.P_.epsilon) / node.P_.tau_j;
  f[6] = (y[S::Z_J] - y[S::E_J]) / node.P_.tau_e;
  f[7] = node.P_.kappa * (y[S::E_J] - y[S::P_J]) / node.P_.tau_p;

  return GSL_SUCCESS;
 }

/* ---------------------------------------------------------------- 
 * Default constructors defining default parameters and state
 * ---------------------------------------------------------------- */
    
mynest::iaf_cond_alpha_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
    kappa   (1.0      ), // dopamine
    fmax    (20.0     ), 
    gain    (1.0     ), 
    bias    (0.0      ),
    epsilon (0.001     ) 
{
  recordablesMap_.create();
}

mynest::iaf_cond_alpha_bias::State_::State_(const Parameters_& p)
  : r(0), bias(0)
{
  y[V_M] = p.E_L;  // initialize to reversal potential
  for ( size_t i = 1 ; i < 5 ; ++i )
    y[i] = 0;
  y[Z_J] = 0.01;
  y[E_J] = 0.01;
  y[P_J] = 0.01;
  //y[BIAS] = 0.0;
}

mynest::iaf_cond_alpha_bias::State_::State_(const State_& s)
  : r(s.r),bias(s.bias)
{
  for ( size_t i = 0 ; i < STATE_VEC_SIZE ; ++i )
    y[i] = s.y[i];
}

mynest::iaf_cond_alpha_bias::State_& mynest::iaf_cond_alpha_bias::State_::operator=(const State_& s)
{
  if ( this == &s )  // avoid assignment to self
    return *this;

  for ( size_t i = 0 ; i < STATE_VEC_SIZE ; ++i )
    y[i] = s.y[i];

  r = s.r;
  return *this;
}

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

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

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

void mynest::iaf_cond_alpha_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, "kappa",      kappa);
  def<nest::double_t>(dd, "bias",      bias);
  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_alpha_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, "kappa",      kappa);
  updateValue<nest::double_t>(dd, "gain",      gain);
  updateValue<nest::double_t>(dd, "bias",      bias);
  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_alpha_bias::State_::get(DictionaryDatum &dd) const
{
  def<double>(dd, nest::names::V_m, y[V_M]); // Membrane potential
}

void mynest::iaf_cond_alpha_bias::State_::set(const DictionaryDatum& dd, const Parameters_&)
{
  updateValue<double>(dd, nest::names::V_m, y[V_M]);
}


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

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

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

mynest::iaf_cond_alpha_bias::~iaf_cond_alpha_bias()
{
  // GSL structs only allocated by init_nodes_(), 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_alpha_bias::init_node_(const Node& proto)
{
  const iaf_cond_alpha_bias& pr = downcast<iaf_cond_alpha_bias>(proto);
  P_ = pr.P_;
  S_ = pr.S_;
}

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

void mynest::iaf_cond_alpha_bias::init_buffers_()
{
  Archiving_Node::clear_history();

  B_.spike_exc_.clear();       // includes resize
  B_.spike_inh_.clear();       // includes resize
  B_.currents_.clear();        // includes resize

  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_alpha_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_alpha_bias::calibrate()
{
  B_.logger_.init();  // ensures initialization in case mm connected after Simulate

  V_.PSConInit_E  = 1.0 * numerics::e / P_.tau_synE;
  V_.PSConInit_I  = 1.0 * numerics::e / P_.tau_synI;
  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_alpha_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);
    }

    // refractoriness and spike generation 
    if ( S_.r )
    {// neuron is absolute refractory
	    --S_.r;
	    S_.y[State_::V_M] = P_.V_reset;  // clamp potential
    }
    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; /* 10k = 1000 * 10 timesteps... */
        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_.kappa * (S_.y[State_::E_J] - S_.y[State_::P_J]) * B_.step_ / P_.tau_p;


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

    // add incoming spikes
    S_.y[State_::DG_EXC] += B_.spike_exc_.get_value(lag) * V_.PSConInit_E;
    S_.y[State_::DG_INH] += B_.spike_inh_.get_value(lag) * V_.PSConInit_I;

    S_.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_alpha_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_alpha_bias::handle(nest::CurrentEvent& e)
{
  assert(e.get_delay() > 0);

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

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

#endif //HAVE_GSL

Loading data, please wait...