Biologically Constrained Basal Ganglia model (BCBG model) (Lienard, Girard 2014)

 Download zip file 
Help downloading and running models
Accession:150206
We studied the physiology and function of the basal ganglia through the design of mean-field models of the whole basal ganglia. The parameterizations are optimized with multi-objective evolutionary algorithm to respect best a collection of numerous anatomical data and electrophysiological data. The main outcomes of our study are: • The strength of the GPe to GPi/SNr connection does not support opposed activities in the GPe and GPi/SNr. • STN and MSN target more the GPe than the GPi/SNr. • Selection arises from the structure of the basal ganglia, without properly segregated direct and indirect pathways and without specific inputs from pyramidal tract neurons of the cortex. Selection is enhanced when the projection from GPe to GPi/SNr has a diffuse pattern.
Reference:
1 . Liénard J, Girard B (2014) A biologically constrained model of the whole basal ganglia addressing the paradoxes of connections and selection. J Comput Neurosci 36:445-68 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Axon; Channel/Receptor; Dendrite; Neural mass;
Brain Region(s)/Organism: Basal ganglia;
Cell Type(s): Neostriatum medium spiny direct pathway GABA cell; Neostriatum medium spiny indirect pathway GABA cell; Neocortex U1 L5B pyramidal pyramidal tract GLU cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Substantia nigra pars reticulata principal GABA cell; Subthalamus nucleus projection neuron; Globus pallidus neuron; Neostriatum fast spiking interneuron; Neostriatum spiny neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: C or C++ program;
Model Concept(s): Parameter Fitting; Pathophysiology; Parkinson's; Winner-take-all; Action Selection/Decision Making;
Implementer(s): Girard, Benoit [girard at isir.upmc.fr]; Liénard, Jean [lienard at isir.upmc.fr];
Search NeuronDB for information about:  Neostriatum medium spiny direct pathway GABA cell; Neostriatum medium spiny indirect pathway GABA cell; Neocortex U1 L5B pyramidal pyramidal tract GLU cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Substantia nigra pars reticulata principal GABA cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
#ifndef BCBG2_HPP
#define BCBG2_HPP

#include "constants.hpp"

#include <iostream>
#include "math.h"
#include "string.h"
#include <vector>
#include "stdio.h"
#include <stdlib.h>
#include <boost/random.hpp>

  template < typename T >
T **Allocate2DArray( int nRows, int nCols)
{
  T **ppi;
  T *pool;
  T *curPtr;
  //(step 1) allocate memory for array of elements of column

  ppi = new T*[nRows];

  //(step 2) allocate memory for array of elements of each row
  pool = new T [nRows * nCols];

  // Now point the pointers in the right place
  curPtr = pool;
  for( int i = 0; i < nRows; i++)
  {
    *(ppi + i) = curPtr;
    curPtr += nCols;
  }
  return ppi;
}

  template < typename T >
void Free2DArray(T** Array)
{
  delete [] *Array;
  delete [] Array;
}

typedef struct _MemorySCN
{
  std::vector <float> S_backup;
  std::vector <std::vector <float> > H_in_RK4_backup;
  std::vector <std::vector <float> > Hp_in_RK4_backup;
} MemorySCN;

typedef struct _MemoryMCN
{
  std::vector <std::vector <float> > S_backup;
  std::vector <std::vector <std::vector <float> > > H_in_RK4_backup;
  std::vector <std::vector <std::vector <float> > > Hp_in_RK4_backup;
} MemoryMCN;

typedef struct _MemoryBCBG2
{
  int tmod_backup;
  int previous_tmod_backup;
  std::vector <MemorySCN > scns;
  std::vector <MemoryMCN > mcns;
} MemoryBCBG2;

class BCBG2
{
  public:
#ifdef TSIROGIANNIS_2010
    BCBG2() : _uni_gaussian(base_generator_t(42u), boost::normal_distribution<>(3.0f, 1.0f)) {}
#else
    //BCBG2() : _uni_one_gaussian(base_generator_t(42u), boost::normal_distribution<>(1.0f, 0.5f)) {} // version with noise
    BCBG2() : _uni_one_gaussian(base_generator_t(42u), boost::normal_distribution<>(1.0f, 0.0f)) {} // no noise
#endif
    void initialize(int nucleus_type, int max_tau, float dt);
    void update(float dt);
    void updateSingleChannelNucleus(int steps);
    void updateSingleChannelNucleus(int steps, int integration_method);
    void updateSingleChannelNucleusWithNoise(int steps, int integration_method);
    void updateSingleChannelNucleusTsiro(int steps, float value);
    float updateSingleChannelNucleusTsiro(int steps); // returns the random cortical firing rate used
    void updateMultiChannelsNucleus(int steps);
    void updateMultiChannelsNucleus_basic_test(int steps, float time);
    void updateMultiChannelsNucleus_humphries_test(int steps, float time);
    float scalevar(float m);
    float scalevarreduc(float m, float reduc);
    void updateMultiChannelsNucleus_exploexplo_test(int steps, float time, int i);
    void hammerizeSingleChannelNucleus(int steps);
    void hammerizeSingleChannelNucleusTsiro(int steps);
    void updateSingleChannelNucleusVana(float input, int bunch);
    void updateSingleChannelNucleusDebug(float time_of_peak_in_s, float dt, float nb_of_s);
    void testPSPSingleChannelNucleus(int pulse, int steps);
    void stabilize_all(int steps);
    bool testTimingResponse();

    float basic_test(int ch, float time);
    float exploexplo_test(int ch, float time, int i);
    float exploexplo_testbis(int ch, float time, float debut_time, float end_time);
    float visualization_humphries_et_al_2006(int ch, float time);

    void updateNucleusCells(int steps);
    void autopiloteNucleusCells(int steps, float* values);

    void save_all();
    void load_all();
    void save_all(MemoryBCBG2& mem);
    void load_all(MemoryBCBG2& mem);

    class SingleChannelNucleus;
    friend class SingleChannelNucleus;
    class MultiChannelsNucleus;
    friend class MultiChannelsNucleus;
    class NucleusCells;
    friend class NucleusCells;

    SingleChannelNucleus* add_single_channel_nucleus(float Smax, float Sini, float vh, float k, int damping, float* dists, float* diams, int compartment_nb, const char* id);
    MultiChannelsNucleus* add_multi_channels_nucleus(int channels_number, float Smax, float Sini, float vh, float k, int connectivity_type, float* dists, float* diams, int compartment_nb, const char* id);
    NucleusCells* add_nucleus_cells(int cells_number, float Smax, float Sini, float vh, float k, float V_rest, float R, float Rm, float Ri, float* dists, float* diams, int compartment_nb, const char* id);

    SingleChannelNucleus* get_single_channel_nucleus(int i) {return SCN[i];};
    MultiChannelsNucleus* get_multi_channels_nucleus(int i) {return MCN[i];};
    NucleusCells* get_nucleus_cells(int i) {return NC[i];};

  protected:
    // general variables of the circuit
    int tmod;
    int previous_tmod;
    int max_tau;
    int nucleus_type;
    float dt;
    float treal;

    // list of nuclei
    std::vector <SingleChannelNucleus *> SCN;
    std::vector <MultiChannelsNucleus *> MCN;
    std::vector <NucleusCells *> NC;
    int n_nuclei; // counter used in the addition of nuclei

    // random generator
    typedef boost::mt19937 base_generator_t;
    typedef boost::variate_generator<base_generator_t, boost::normal_distribution<> > rand_t;
    rand_t _uni_one_gaussian;

    // for use in store/restore function
    int tmod_backup;
    int previous_tmod_backup;
};



class Constants {
  public:
    // defining some constants
    static const float A_AMPA = 0.001 * 5.43656365692;        // 2*e
    static const float D_AMPA = 0.65238763883017081;  // (1/4)*e
    static const float A_NMDA = 0.001 * 0.27182818284590454;  // 0.1*e
    static const float D_NMDA = 0.027182818284590453; // (1/100)*e
    static const float A_GABA = 0.001 * 5.4365636569180902;   // 2*e
    static const float D_GABA = 0.45304697140984085;  // (1/6)*e
};


class BCBG2::MultiChannelsNucleus
{
  public:
    friend class BCBG2;
    MultiChannelsNucleus(BCBG2& bg) : bg_(bg) {}
    void initialize_multi_channels_nucleus(int ch_n, float Smax, float Sini, float vh, float k, int connectivity_type, float* dists, float* diams, int compartment_nb, const char* id);
    void update_multi_channels_nucleus_stabilize(int steps);
    void update_multi_channels_nucleus_evo();
    void set_afferent(float A, float D, int Sign, float C, int T, float distance, MultiChannelsNucleus* N);
    void set_afferent(float A, float D, int Sign, float C, int T, float connexion_scheme, float distance, MultiChannelsNucleus* N);
    float compute_distance_factor(float distance);
    void initialize_new_afferent();
    void display(std::ostream&);
    void set_dt(float value);

    // various getter/setter
    float get_S(int tau_in, int ch_n) { return S[tau_in][ch_n]; }
    // Warning ! get_S() returns S for _last step done_ not for  _current step_
    float get_S(int ch_n) { return S[(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau][ch_n]; } 
    void set_S(float value, int tau_in, int ch_n) { S[tau_in][ch_n] = value; }
    void set_S(float value, int ch_n) { S[bg_.tmod][ch_n] = value; }
    int get_afferents_number() { return n; } 
    int get_channels_number() { return ch_n; } 
    char* get_name() { return id; }

    void save(MemoryMCN& mem);
    void load(MemoryMCN& mem);


  protected:
    // to access the constants of the circuit
    BCBG2& bg_;

    // general variables of the multi_channels_nucleus
    int ch_n; // number of concurrent channels
    std::vector <std::vector <float> > S; //output frequency
    float Smax;
    float vh;
    float k;
    float kvh;
    float dt;
    char *id;

    float Rm;                      // Membrane resistance
    float Ri;                      // Internal resistivity
    float membrane_area;           // Surface of the neuron
    std::vector <float> distances; // used to reconstruct a simple model of dendrites
    std::vector <float> diameters; // used to reconstruct a simple model of dendrites

    // afferents
    float k1,k2,k3,k4,k5,k6,k7;
    std::vector <std::vector <float> > H_in;
    std::vector <std::vector <std::vector <float> > > H_in_RK4;
    std::vector <std::vector <std::vector <float> > > Hp_in_RK4;
    std::vector <std::vector <std::vector <float> > > Hpp_in_RK4;
    std::vector <float> A_in;
    std::vector <float> D_in;
    std::vector <float> C_in;
    std::vector <float> ADC_in;
    std::vector <float> DD_in;
    std::vector <float> D2_in;
    std::vector <float> dtADC_in;
    std::vector <float> dtDD_in;
    std::vector <float> dtD2_in;
    std::vector <int> T_in;
    std::vector <int> Sign_in;
    std::vector <float> conn_scheme_in;
    std::vector <MultiChannelsNucleus *> N_in;
    int n; // counter used in the addition of afferents

    // to store/restore the state of a circuit
    std::vector <std::vector <float> > S_backup; //output frequency
    std::vector <std::vector <std::vector <float> > > H_in_RK4_backup;
    std::vector <std::vector <std::vector <float> > > Hp_in_RK4_backup;
};

class BCBG2::NucleusCells
{
  // this class implements the cell-level approximation
  // NOT TO BE USED - THIS IS NOT FINISHED
  // Use instead the classes SingleChannelNucleus (single channel allows quick simulation) and MultiChannelsNucleus (the multi channels version)
  public:

    struct Target {
      BCBG2::NucleusCells * target_nucleus;
      int cell_number;
      int ion_channel_number;
      float T; 
      float distance_factor;
    };
    
    struct IonChannel {
      float A;           // | 
      float D;           // | 
      float Vrev;        // | 
      bool is_nmda;      // | these informations are duplicated for each cell
      float time_to_add; // |
      float to_add;      // |
    };

    struct Spike {
      int target_number;
      float t;
      float t_expiration; // set to the same value as IonChannel
    };
    

    friend class BCBG2;
    NucleusCells(BCBG2& bg) : bg_(bg), is_firing(rng) {}

    void initialize_cells(int cells_n, float Smax, float Sini, float V_threshold, float k, float V_rest, float R, float Rm, float Ri, float* dists, float* diams, int compartment_nb, const char* id, int nucleus_number);
    int affect_ion_channel(float A, float D, float Vrev, bool is_nmda);
    void add_efferent_axon(NucleusCells* nc, float cs, int T, float C, float distance_factor, int ion_channel_n);
    void add_afferent_axon(NucleusCells* nc, float cs, int T, float C, float distance_factor, int ion_channel_n);
    void add_afferent_synapse(NucleusCells* source_nucleus, float cs, int T, float C, float distance_factor, int ion_channel_n);
    bool set_afferent(float A, float D, int S, float C, int T, float cs, float distance, NucleusCells* N);
    bool add(float A, float D, float Vrev, float C, int T, float cs, float distance, NucleusCells* N);
    bool nmda(float A, float D, float Vrev, float C, int T, float cs, float distance, NucleusCells* N);
    bool add_ampa_gaba_nmda(float A, float D, float Vrev, float C, int T, float cs, float distance, NucleusCells* N, bool is_nmda);
    void set_dt(float value);
    void fast_add_g(int cell_i, int ion_channel_i, float distance_factor, float t_sent);
    void add_g(int cell_i, int ion_channel_i, float distance_factor, float t_sent);
    void update_cells();
    void set_rate(float rate, bool reset);
    void set_rate(float rate); // wrapper around previous function, with reset=true
    void prepare_next_iteration();
    void display(std::ostream& os);
    void display_summary(std::ostream& os);

    // various getter/setter
    int get_nucleus_number() {
      return this->nucleus_number;
    }

    float get_S(int cell_i) { float hz=0; for (int i=0; i<times_of_spikes[cell_i].size(); i++) {if (bg_.treal - times_of_spikes[cell_i][i] < 1.) {hz+=1.;} } return (hz); } 
    float get_meanS()
    {
      float hz=0; 
      for (int cell_i=0; cell_i<cells_n; cell_i++) { 
        for (int i=0; i<times_of_spikes[cell_i].size(); i++) {
          if (bg_.treal - times_of_spikes[cell_i][i] < 1.) {
            hz+=1.;
          }
        }
      }
      return (hz/((float)cells_n));
    }
    void debugg()
    { 
      float now = bg_.treal; 
      std::cout << "NOW=" << now << std::endl;
      std::cout << "DT=" << bg_.dt << std::endl;
      for (int cell_i=0; cell_i<cells_n; cell_i++) {
        std::cout << cell_i << " : " << last_spiked[cell_i] << std::endl;
      }
    } 
    float get_instantaneous_S3()
    {
      float hz=0; 
      for (int cell_i=0; cell_i<cells_n; cell_i++) { 
        for (int i=0; i<times_of_spikes[cell_i].size(); i++) {
          if (bg_.treal - times_of_spikes[cell_i][i] < 0.1) {
            hz+=1.;
          }
        }
      }
      return (10.*hz/((float)cells_n));
    }
    float get_instantaneous_S2()
    { 
      float now = bg_.treal; 
      int nb_discharging = 0;
      for (int cell_i=0; cell_i<cells_n; cell_i++) {
        if (now - last_spiked[cell_i] < absolute_refractory_period-dt) { // refractory period dependent of maximum rate
          nb_discharging++;
        } 
      }
      return (((float)nb_discharging) / ((float)cells_n)) * (absolute_refractory_period/dt);
    } 
    float get_instantaneous_S()
    { 
      float now = bg_.treal; 
      float isi=1./Smax; 
      float no_spiking_cells = 0.;
      for (int cell_i=0; cell_i<cells_n; cell_i++) {
        if (now - last_spiked[cell_i] >= 1./Smax + 2.*dt) { 
          isi += now - last_spiked[cell_i]; 
          no_spiking_cells = no_spiking_cells + 1.;
        } 
      }
      if (no_spiking_cells == 0.) {
        return 0.;
      }
      return (no_spiking_cells)/isi; 
    } 
    float get_meanV() { float v=0; for (int cell_i=0; cell_i<cells_n; cell_i++) { v+=V[cell_i];} return (v/((float)cells_n)); } // TODO TODO TODO TODO TODO original
    float get_meanV_no_silenced() { float v=0; for (int cell_i=0; cell_i<cells_n; cell_i++) { if (V[cell_i] != V_rest) {v+=V[cell_i];}} return (v/((float)cells_n)); } // corrected to exclude silenced
    float get_meanG() { float gg=0; for (int cell_i=0; cell_i<cells_n; cell_i++) { for (int ion_channel_i=0; ion_channel_i<ion_channels.size(); ion_channel_i++) {gg+=g[cell_i][ion_channel_i];}} return (gg/((float)cells_n*((float)ion_channels.size()))); } 
    float get_mean_inc_sum_v() { float mean_inc_sum_v=0; for (int cell_i=0; cell_i<cells_n; cell_i++) { mean_inc_sum_v+=inc_sum_v[cell_i];} return (mean_inc_sum_v/((float)cells_n)); } 
    float get_V(int cells_i) { return V[cells_i]; } 
    int get_afferents_number() { return n; } 
    char* get_name() { return id; }
    int get_cells_number() {return cells_n; }
    float get_ion_channel_time_to_add(int ion_channel_i) {return ion_channels[ion_channel_i].time_to_add; } // TODO this is dirty, all cells should share ion_channels (see struct IonChannel)

    void save();
    void restore();

  protected:
    // to access the constants of the circuit
    BCBG2& bg_;
    boost::mt19937 rng;
    boost::uniform_01<boost::mt19937> is_firing;

    // general variables of the nucleus_cells
    int cells_n;                   // number of cells
    std::vector <std::vector <float> > S; //output frequency
    float Smax;
    std::vector <float> V_threshold;
    std::vector <float> V_reset;                 // Reset potential
    std::vector <float> V;         // Mean voltage of the neurons
    std::vector <float> inc_sum_v;     // Mean input voltage of the neurons from ionic channels
    std::vector <std::vector <float> > g;
    std::vector <std::vector <float> > next_g;
    float V_rest;                  // Mean resting voltage
    float R;                       // Membrane input resistance
    float Rm;                      // Membrane resistance
    float Ri;                      // Internal resistivity
    float membrane_area;           // Surface of the neuron
    std::vector <float> distances; // used to reconstruct a simple model of dendrites
    std::vector <float> diameters; // used to reconstruct a simple model of dendrites
    float absolute_refractory_period;
    float time_to_vanish;
    float dt;
    char *id;
    int nucleus_number;

    // afferents
    float k1,k2,k3,k4,k5,k6,k7;
    std::vector <float> Vrev_in;            // Reversal Potential associated with synapse
    std::vector <float> Distance_factor_in; // Distance factor computed as a single or multi compartment with sealed end
    std::vector <NucleusCells *> N_in;
    int n; // counter used in the addition of afferents

    std::vector <float> last_spiked; // TODO
    std::vector <struct IonChannel> ion_channels; // TODO
    std::vector <std::vector <struct Target> > targets; // TODO
    std::vector <std::vector <struct Spike> > spikes; // TODO
    std::vector <std::vector <float> > times_of_spikes; // TODO
};


class BCBG2::SingleChannelNucleus
{
  public:
    friend class BCBG2;
    SingleChannelNucleus(BCBG2& bg) : bg_(bg) {}
    void initialize_single_channel_nucleus(float Smax, float Sini, float vh, float k, int damping, float* dists, float* diams, int compartment_nb, const char* id);
    void update_single_channel_nucleus_stabilize(int steps);
    void update_single_channel_nucleus_evo();
    float get_all_inputs()  // for debugging purpose
    {
      float r = 0;
      r += A_in[0] * C_in[0] * N_in[0]->get_S((bg_.tmod - T_in[0] - 1 + bg_.max_tau ) % bg_.max_tau);
      r += A_in[2] * C_in[2] * N_in[2]->get_S((bg_.tmod - T_in[2] - 1 + bg_.max_tau ) % bg_.max_tau);
      return r;      
    }
    void update_single_channel_nucleus_euler();
    void update_single_channel_nucleus_rk3();
    void update_single_channel_nucleus_evo_Hp();
    void update_single_channel_nucleus_vana();
    void update_single_channel_nucleus_damping_vana();
    void update_single_channel_nucleus();
    void set_afferent(float A, float D, int Sign, float C, int T, float distance, SingleChannelNucleus* N);
    void set_afferent(float A, float D, int Sign, float C, int T, float to_be_ignored, float distance, SingleChannelNucleus* N); // Dummy function, calls other "set_afferent"
    void set_afferent(float nu, int T, SingleChannelNucleus* N);
    float compute_distance_factor(float distance);
    void initialize_new_afferent();
    void display(std::ostream&);
    void display_Hs() {for (int i=0;i<n;i++) std::cout << "H=" << H_in_RK4[i][0] << " Hp=" << Hp_in_RK4[i][0] << std::endl; std::cout <<"-------------" << std::endl;}; // for debug
    void set_dt(float value);

    void save();
    void load();
    void save(MemorySCN& mem);
    void load(MemorySCN& mem);

    // various getter/setter
    float get_S(int tau_in) { if (is_damped) {return Phi_in_RK4[tau_in];} else {return S[tau_in];} }
    // Warning ! get_S() returns S for _last step done_ not for  _current step_
    float get_S() { if (is_damped) {return Phi_in_RK4[(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau];} else {return S[(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau];} } 
    void set_S(float value, int tau_in) { S[tau_in] = value; }
    void set_S(float value) { S[bg_.tmod] = value; }
    float get_H(int i) { return H_in_RK4[i][(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau]; } 
    float get_V() {float V=0; for (int i=0;i<n;i++) V+= H_in_RK4[i][(bg_.tmod  + bg_.max_tau) % bg_.max_tau] * Sign_in[i]; return V; } 
    float get_Hp(int i) { return Hp_in_RK4[i][(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau]; } 
    void scale_Hp(int i, float value) { Hp_in_RK4[i][(bg_.tmod - 1 + bg_.max_tau) % bg_.max_tau] *= value; } 
    int get_afferents_number() { return n; } 
    char* get_name() { return id; }


  protected:
    // to access the constants of the circuit
    BCBG2& bg_;

    // general variables of the multi_channels_nucleus
    std::vector <float> S; //output frequency
    float Smax;
    float vh;
    float k;
    float kvh;
    float dt;
    char *id;
    int is_damped;

    float Rm;                      // Membrane resistance
    float Ri;                      // Internal resistivity
    float membrane_area;           // Surface of the neuron
    std::vector <float> distances; // used to reconstruct a simple model of dendrites
    std::vector <float> diameters; // used to reconstruct a simple model of dendrites

    // afferents
    std::vector <float> H_in;
    float k1,k2,k3,k4,k5,k6,k7;
    float k1bis,k2bis,k3bis;
    std::vector <float> Phi_in_RK4;
    std::vector <float> Phip_in_RK4;
    std::vector <std::vector <float> > H_in_RK4;
    std::vector <std::vector <float> > Hp_in_RK4;
    std::vector <std::vector <float> > Hpp_in_RK4;
    std::vector <float> H_in_old;
    std::vector <float> Hp_in;
    std::vector <float> Hpp_in;
    std::vector <float> Hp_in_old;
    std::vector <float> Hpp_in_old;
    std::vector <float> A_in;
    std::vector <float> D_in;
    std::vector <float> C_in;
    std::vector <float> ADC_in;
    std::vector <float> DD_in;
    std::vector <float> D2_in;
    std::vector <float> nu_in;
    std::vector <float> dtADC_in;
    std::vector <float> dtDD_in;
    std::vector <float> dtD2_in;
    std::vector <int> T_in;
    std::vector <int> Sign_in;
    std::vector <SingleChannelNucleus *> N_in;
    int n; // counter used in the addition of afferents

    // to store/restore the state of a circuit
    std::vector <float> S_backup; //output frequency
    std::vector <std::vector <float> > H_in_RK4_backup;
    std::vector <std::vector <float> > Hp_in_RK4_backup;
};


#endif

Loading data, please wait...