Merge pull request #2773 from WillemWybo/user_configurable_compartmen…

…tal_model_initializtion

User configurable voltage initialization for the compartmental model

models/cm_default.cpp

@@ -223,6 +223,8 @@ nest::cm_default::set_status( const DictionaryDatum& statusdict )
 void
 nest::cm_default::add_compartment_( DictionaryDatum& dd )
 {
+  dd->clear_access_flags();
+
   if ( dd->known( names::params ) )
   {
     c_tree_.add_compartment(
@@ -232,10 +234,14 @@ nest::cm_default::add_compartment_( DictionaryDatum& dd )
   {
     c_tree_.add_compartment( getValue< long >( dd, names::parent_idx ) );
   }
+
+  ALL_ENTRIES_ACCESSED( *dd, "cm_default::add_compartment_", "Unread dictionary entries: " );
 }
 void
 nest::cm_default::add_receptor_( DictionaryDatum& dd )
 {
+  dd->clear_access_flags();
+
   const long compartment_idx = getValue< long >( dd, names::comp_idx );
   const std::string receptor_type = getValue< std::string >( dd, names::receptor_type );
 
@@ -257,6 +263,8 @@ nest::cm_default::add_receptor_( DictionaryDatum& dd )
   {
     compartment->compartment_currents.add_synapse( receptor_type, syn_idx );
   }
+
+  ALL_ENTRIES_ACCESSED( *dd, "cm_default::add_receptor_", "Unread dictionary entries: " );
 }
 
 void

models/cm_default.h

@@ -124,19 +124,28 @@ variables, use the receptor index ``{state_variable_name}{receptor_index}``:
 
 .. code-block:: Python
 
-    mm = nest.Create('multimeter', 1, {'record_from': ['v_comp0'}, ...})
+    mm = nest.Create('multimeter', 1, {'record_from': ['v_comp0', ...]})
 
 Current generators can be connected to the model. In this case, the receptor
 type is the compartment index:
 
 .. code-block:: Python
 
     dc = nest.Create('dc_generator', {...})
-    nest.Connect(dc, cm, syn_spec={..., 'receptor_type': 0}
+    nest.Connect(dc, cm, syn_spec={..., 'receptor_type': 0})
 
 Parameters
 ++++++++++
 
+Note that the compartmental model does not explicitly ensure that units are consistent.
+Therefore, it is on the user to ensure that units are consistent throughout the model.
+The quantities that have fixed units are membrane voltage [mV] and time [ms].
+Other units need to be consistent: if e.g. conductances are in uS, that means
+that the associated currents will be uS*mV = nA. By consequence, the capacitance needs to
+be in nF to ensure that the capacitive current is also in nA. This further means
+that the connection weights to receptors are in uS, and that the amplitudes of current
+injectors are in nA.
+
 The following parameters can be set in the status dictionary.
 
 =========== ======= ===========================================================
@@ -146,10 +155,11 @@ The following parameters can be set in the status dictionary.
 The following parameters can be used when adding compartments using ``SetStatus()``
 
 =========== ======= ===============================================================
- C_m        uF      Capacitance of compartment (default: 1 uF)
+ C_m        nF      Capacitance of compartment (default: 1 nF)
  g_C        uS      Coupling conductance with parent compartment (default: 0.01 uS)
  g_L        uS      Leak conductance of the compartment (default: 0.1 uS)
  e_L        mV      Leak reversal of the compartment (default: -70. mV)
+ v_comp     mV      Initialization voltage of the compartment (default: -75. mV)
 =========== ======= ===============================================================
 
 Ion channels and receptor types for the default model are hardcoded.

models/cm_tree.cpp

@@ -28,11 +28,11 @@ nest::Compartment::Compartment( const long compartment_index, const long parent_
   , comp_index( compartment_index )
   , p_index( parent_index )
   , parent( nullptr )
-  , v_comp( 0.0 )
   , ca( 1.0 )
   , gc( 0.01 )
   , gl( 0.1 )
   , el( -70. )
+  , v_comp( el )
   , gg0( 0.0 )
   , ca__div__dt( 0.0 )
   , gl__div__2( 0.0 )
@@ -42,11 +42,11 @@ nest::Compartment::Compartment( const long compartment_index, const long parent_
   , gg( 0.0 )
   , hh( 0.0 )
   , n_passed( 0 )
+  , compartment_currents( v_comp )
 {
-  v_comp = el;
-
-  compartment_currents = CompartmentCurrents();
+  compartment_currents = CompartmentCurrents( v_comp );
 }
+
 nest::Compartment::Compartment( const long compartment_index,
   const long parent_index,
   const DictionaryDatum& compartment_params )
@@ -55,11 +55,11 @@ nest::Compartment::Compartment( const long compartment_index,
   , comp_index( compartment_index )
   , p_index( parent_index )
   , parent( nullptr )
-  , v_comp( 0.0 )
   , ca( 1.0 )
   , gc( 0.01 )
   , gl( 0.1 )
   , el( -70. )
+  , v_comp( el )
   , gg0( 0.0 )
   , ca__div__dt( 0.0 )
   , gl__div__2( 0.0 )
@@ -69,15 +69,19 @@ nest::Compartment::Compartment( const long compartment_index,
   , gg( 0.0 )
   , hh( 0.0 )
   , n_passed( 0 )
+  , compartment_currents( v_comp )
 {
+  compartment_params->clear_access_flags();
+
   updateValue< double >( compartment_params, names::C_m, ca );
   updateValue< double >( compartment_params, names::g_C, gc );
   updateValue< double >( compartment_params, names::g_L, gl );
   updateValue< double >( compartment_params, names::e_L, el );
+  updateValue< double >( compartment_params, names::v_comp, v_comp );
 
-  v_comp = el;
+  compartment_currents = CompartmentCurrents( v_comp, compartment_params );
 
-  compartment_currents = CompartmentCurrents( compartment_params );
+  ALL_ENTRIES_ACCESSED( *compartment_params, "compartment_params", "Unread dictionary entries: " );
 }
 
 void

models/cm_tree.h

@@ -51,44 +51,46 @@ namespace nest
 class Compartment
 {
 private:
-  // aggragators for numerical integration
+  //! aggragators for numerical integration
   double xx_;
   double yy_;
 
 public:
-  // compartment index
+  //! compartment index
   long comp_index;
-  // parent compartment index
+  //! parent compartment index
   long p_index;
-  // tree structure indices
+  //! tree structure indices
   Compartment* parent;
   std::vector< Compartment > children;
-  // vector for synapses
-  CompartmentCurrents compartment_currents;
 
-  // buffer for currents
+  //! buffer for currents
   RingBuffer currents;
-  // voltage variable
-  double v_comp;
-  // electrical parameters
+  //! electrical parameters
   double ca; // compartment capacitance [uF]
   double gc; // coupling conductance with parent (meaningless if root) [uS]
   double gl; // leak conductance of compartment [uS]
   double el; // leak current reversal potential [mV]
-  // auxiliary variables for efficienchy
+
+  //! voltage variable
+  double v_comp;
+
+  //! auxiliary variables for efficienchy
   double gg0;
   double ca__div__dt;
   double gl__div__2;
   double gc__div__2;
   double gl__times__el;
-  // for numerical integration
+  //! for numerical integration
   double ff;
   double gg;
   double hh;
-  // passage counter for recursion
+  //! passage counter for recursion
   int n_passed;
 
-  // constructor, destructor
+  //! vector for synapses
+  CompartmentCurrents compartment_currents;
+
   Compartment( const long compartment_index, const long parent_index );
   Compartment( const long compartment_index, const long parent_index, const DictionaryDatum& compartment_params );
   ~Compartment() {};
@@ -97,10 +99,10 @@ class Compartment
   void pre_run_hook();
   std::map< Name, double* > get_recordables();
 
-  // matrix construction
+  //! matrix construction
   void construct_matrix_element( const long lag );
 
-  // maxtrix inversion
+  //! maxtrix inversion
   inline void gather_input( const std::pair< double, double >& in );
   inline std::pair< double, double > io();
   inline double calc_v( const double v_in );
@@ -148,7 +150,7 @@ nest::Compartment::calc_v( const double v_in )
 class CompTree
 {
 private:
-  /*
+  /**
   structural data containers for the compartment model
   */
   mutable Compartment root_;
@@ -158,30 +160,38 @@ class CompTree
 
   long size_ = 0;
 
-  // recursion functions for matrix inversion
+  //! recursion functions for matrix inversion
   void solve_matrix_downsweep( Compartment* compartment_ptr, std::vector< Compartment* >::iterator leaf_it );
   void solve_matrix_upsweep( Compartment* compartment, double vv );
 
-  // functions for pointer initialization
+  //! functions for pointer initialization
   void set_parents();
   void set_compartments();
   void set_leafs();
 
 public:
-  // constructor, destructor
   CompTree();
   ~CompTree() {};
 
-  // initialization functions for tree structure
+  //! add a compartment to the tree structure
   void add_compartment( const long parent_index );
   void add_compartment( const long parent_index, const DictionaryDatum& compartment_params );
   void add_compartment( Compartment* compartment, const long parent_index );
+
+  //! initialize the tree for simulation
   void pre_run_hook();
+
+  //! fix all pointers in the tree, the tree structure should not be modified between calling
+  //! this function and starting the simulation
   void init_pointers();
+
+  //! associate each receptor with a spike buffer
   void set_syn_buffers( std::vector< RingBuffer >& syn_buffers );
+
+  //! make all state variables accessible for recording
   std::map< Name, double* > get_recordables();
 
-  // get a compartment pointer from the tree
+  //! get a compartment pointer from the tree
   Compartment* get_compartment( const long compartment_index ) const;
   Compartment* get_compartment( const long compartment_index, Compartment* compartment, const long raise_flag ) const;
   Compartment* get_compartment_opt( const long compartment_indx ) const;
@@ -191,23 +201,23 @@ class CompTree
     return &root_;
   };
 
-  // get tree size (number of compartments)
+  //! get tree size (number of compartments)
   long
   get_size() const
   {
     return size_;
   };
 
-  // get voltage values
+  //! get voltage values
   std::vector< double > get_voltage() const;
   double get_compartment_voltage( const long compartment_index );
 
-  // construct the numerical integration matrix and vector
+  //! construct the numerical integration matrix and vector
   void construct_matrix( const long lag );
-  // solve the matrix equation for next timestep voltage
+  //! solve the matrix equation for next timestep voltage
   void solve_matrix();
 
-  // print function
+  //! print function
   void print_tree() const;
 }; // CompTree
 

nestkernel/nest_names.cpp

@@ -622,6 +622,7 @@ const Name use_wfr( "use_wfr" );
 const Name v( "v" );
 const Name V_act_NMDA( "V_act_NMDA" );
 const Name V_clamp( "V_clamp" );
+const Name v_comp( "v_comp" );
 const Name V_epsp( "V_epsp" );
 const Name V_m( "V_m" );
 const Name V_min( "V_min" );

nestkernel/nest_names.h

@@ -650,6 +650,7 @@ extern const Name use_wfr;
 extern const Name v;
 extern const Name V_act_NMDA;
 extern const Name V_clamp;
+extern const Name v_comp;
 extern const Name V_epsp;
 extern const Name V_m;
 extern const Name V_min;

pynest/examples/compartmental_model/receptors_and_current.py

@@ -41,12 +41,14 @@
     "g_C": 0.0,  # soma has no parent
     "g_L": 1.0,  # [nS] Leak conductance
     "e_L": -70.0,  # [mV] leak reversal
+    "v_comp": -70.0,  # [mV] voltage initialization
 }
 dend_params = {
     "C_m": 0.1,  # [pF] Capacitance
     "g_C": 0.1,  # [nS] Coupling conductance to parent (soma here)
     "g_L": 0.1,  # [nS] Leak conductance
     "e_L": -70.0,  # [mV] leak reversal
+    "v_comp": -70.0,  # [mV] voltage initialization
 }
 
 ###############################################################################

pynest/examples/compartmental_model/two_comps.py

@@ -42,6 +42,7 @@
     "g_C": 0.0,  # soma has no parent
     "g_L": 8.924572508,  # [nS] Leak conductance
     "e_L": -75.0,  # [mV] leak reversal
+    "v_comp": -75.0,  # [mV] compartment initialization
     # ion channel params
     "gbar_Na": 4608.698576715,  # [nS] Na maximal conductance
     "e_Na": 60.0,  # [mV] Na reversal
@@ -58,13 +59,15 @@
     "g_C": 1.255439494,
     "g_L": 0.192992878,
     "e_L": -75.0,
+    "v_comp": -75.0,
 }
 dend_params_active = {
     # passive parameters
     "C_m": 1.929929,  # [pF] Capacitance
     "g_C": 1.255439494,  # [nS] Coupling conductance to parent (soma here)
     "g_L": 0.192992878,  # [nS] Leak conductance
     "e_L": -70.0,  # [mV] leak reversal
+    "v_comp": -70.0,  # [mV] compartment initialization
     # ion channel params
     "gbar_Na": 17.203212493,  # [nS] Na maximal conductance
     "e_Na": 60.0,  # [mV] Na reversal