Make rxd.Species keep the corresponding SpeciesOnRegion alive. (#1270)

* Make rxd.Species keep corresponding SpeciesOnRegion alive.

* Add comment.

share/lib/python/neuron/rxd/species.py

@@ -1253,6 +1253,7 @@ def __init__(self, regions=None, d=0, name=None, charge=0, initial=None, atolsca
         self.initial = initial
         self._atolscale = atolscale
         self._ecs_boundary_conditions = ecs_boundary_conditions
+        self._species_on_region = weakref.WeakKeyDictionary()
         if represents and not _ontology_id.match(represents):
             raise RxDException("The represents=%s is not valid CURIE" % represents)
         else:
@@ -1529,12 +1530,16 @@ def __getitem__(self, r):
         The resulting object is a SpeciesOnRegion.
         This is useful for defining reaction schemes for MultiCompartmentReaction."""
         if isinstance(r, region.Region) and r in self._regions:
-            return SpeciesOnRegion(self, r)
+            if r not in self._species_on_region:
+                 self._species_on_region[r] = SpeciesOnRegion(self, r)
+            return self._species_on_region[r]
         elif isinstance(r, region.Extracellular):
             if not hasattr(self,'_extracellular_instances'):
                 initializer._do_init()
             if r in self._extracellular_instances:
-                return SpeciesOnExtracellular(self, self._extracellular_instances[r])
+                if r not in self._species_on_region:
+                    self._species_on_region[r] = SpeciesOnExtracellular(self, self._extracellular_instances[r])
+                return self._species_on_region[r]
         raise RxDException('no such region')
 
     def _update_node_data(self):
@@ -1991,12 +1996,16 @@ def __getitem__(self, r):
         The resulting object is a ParameterOnRegion or ParameterOnExtracellular.
         This is useful for defining reaction schemes for MultiCompartmentReaction."""
         if isinstance(r, region.Region) and r in self._regions:
-            return ParameterOnRegion(self, r)
+            if r not in self._species_on_region:
+                self._species_on_region[r] = ParameterOnRegion(self, r)
+            return self._species_on_region[r]
         elif isinstance(r, region.Extracellular):
             if not hasattr(self,'_extracellular_instances'):
                 initializer._do_init()
             if r in self._extracellular_instances:
-                return ParameterOnExtracellular(self, self._extracellular_instances[r])
+                if r not in self._species_on_region:
+                    self._species_on_region[r] = ParameterOnExtracellular(self, self._extracellular_instances[r])
+                return self._species_on_region[r]
         raise RxDException('no such region')
 
 class ParameterOnRegion(SpeciesOnRegion):

test/rxd/test_multicompartment_reactions.py

@@ -53,9 +53,8 @@ def test_multicompartment_reactions(neuron_instance):
 
     ip3 = rxd.Species(cyt, d=ip3Diff, name="ip3", initial=ip3_init)
     ip3r_gate_state = rxd.Species(cyt_er_membrane, name="gate", initial=0.8)
-    h_gate = ip3r_gate_state[cyt_er_membrane]
     minf = ip3[cyt] * 1000.0 * ca[cyt] / (ip3[cyt] + kip3) / (1000.0 * ca[cyt] + kact)
-    k = gip3r * (minf * h_gate) ** 3
+    k = gip3r * (minf * ip3r_gate_state[cyt_er_membrane]) ** 3
 
     ip3r = rxd.MultiCompartmentReaction(ca[er], ca[cyt], k, k, membrane=cyt_er_membrane)
 
@@ -69,8 +68,10 @@ def test_multicompartment_reactions(neuron_instance):
     leak = rxd.MultiCompartmentReaction(
         ca[er], ca[cyt], gleak, gleak, membrane=cyt_er_membrane
     )
-
-    ip3rg = rxd.Rate(h_gate, (1.0 / (1 + 1000.0 * ca[cyt] / (0.3)) - h_gate) / ip3rtau)
+    # test the SpeciesOnRegion remains in scope
+    ip3rg = rxd.Rate(ip3r_gate_state[cyt_er_membrane], 
+                     (1.0 / (1 + 1000.0 * ca[cyt] / (0.3)) - 
+                     ip3r_gate_state[cyt_er_membrane]) / ip3rtau)
 
     h.finitialize(-65)