Add support for different geometries on different sections of the sam…

…e region. (#1000)

* Add support for different geometries on different sections of the same region.
Fix an issue with multicompartment reactions with regions with some different sections.

* Add test data.

* Fix test for python2.

* Fix missing `sec` following @ramcdougal's review

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

@@ -18,7 +18,7 @@
 except:
   pass
 from .rangevar import RangeVar
-from .geometry import membrane, inside, Shell, FractionalVolume, FixedCrossSection, FixedPerimeter, ScalableBorder, DistributedBoundary
+from .geometry import membrane, inside, Shell, FractionalVolume, FixedCrossSection, FixedPerimeter, ScalableBorder, DistributedBoundary, MultipleGeometry
 from .plugins import set_solver
 # deprecated:
 # from geometry import ConstantArea, ConstantVolume

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

@@ -416,3 +416,77 @@ def volumes1d(self, sec):
             vols[iseg] = volume
 
         return vols
+
+class MultipleGeometry(RxDGeometry):
+    """ support for different geometries on different sections of a region.
+    
+    Example use:
+    - for radial diffusion in a dendrite (dend) with longitudinal diffusion
+      from a spine (spine). The region for the outer shell of the dendrite
+      (0.8,1] should include the while spine [0,1];
+      MultipleGeometry(secs=[dend,spine], geos=[Shell(0.8,1), rxd.inside])
+
+    Args:   
+    sections (list, optional) a list or list-of-lists of sections where the
+        corresponding geometry should be used. If None the same geometry used
+        for all sections, otherwise the list must be the same length as the
+        geos list.
+        If None is in the list, then the corresponding geometry in geos is used
+        as a default value for any section not included in the lists.
+    geometries (list) a list of geometries that are used for the corresponding
+        list of sections in secs. All geometries must be volumes or all
+        all geometries must be areas.
+
+    """
+    def __init__(self, secs=None, geos=None):
+        self._secs = {}
+        self._default = None
+        if not secs:
+            if isinstance(geos,list):
+                self._default = geos[0]
+            elif isinstance(geos,RxDGeometry):
+                self._default = geos
+            else:
+                raise RxDException("MultipleGeometry requires a list-of-lists of sections and their corresponding geometry")
+        else:
+            assert(len(secs) == len(geos))
+            if all([g.is_area() for g in geos]):
+                self.is_area = _always_true
+                self.is_volume = _always_false
+            elif all([g.is_volume() for g in geos]):
+                self.is_area = _always_false
+                self.is_volume = _always_true
+            else:
+                raise RxDException("MultipleGeometry requires all geometries are areas or all geometries are volumes")
+            for s, g in zip(secs,geos):
+                if not s:
+                    self._default = g
+                elif isinstance(s,list):
+                    self._secs[h.SectionList(s)] = g
+                else:
+                    self._secs[h.SectionList([s])] = g
+    def __repr__(self):
+        secs = [[s for s in sl] for sl in self._secs]
+        geos = [self._secs[sl] for sl in self._secs]
+        return 'MultipleGeometry(secs=%r, geos=%r)' % (secs, geos)
+
+    def _get_geo(self, sec):
+        if not isinstance(sec, nrn.Section):
+            sec = sec._sec
+        for sl in self._secs:
+            if sec in sl:
+                geo = self._secs[sl]
+                break
+        else:
+            if self._default:
+                geo = self._default
+            else:
+                raise RxDException('MultipleGeometry is not defined on section %r' % sec)
+        return geo
+
+    def volumes1d(self, sec):
+        return self._get_geo(sec).volumes1d(sec)
+    def surface_areas1d(self, sec):
+        return self._get_geo(sec).surface_areas1d(sec)
+    def neighbor_areas1d(self, sec):
+        return self._get_geo(sec).neighbor_areas1d(sec)

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

@@ -68,6 +68,11 @@ def __init__(self, regions):
                 _c_region_lookup[rptr] = [self]
 
     def add_reaction(self, rptr, region):
+        # for multicompartment reaction -- check all regions are present
+        if rptr() and hasattr(rptr(),'_changing_species'):
+            for sptr in rptr()._changing_species:
+                if sptr() and hasattr(sptr(),'_region') and sptr()._region not in self._regions:
+                    return
         if rptr in self._react_regions:
             if region not in self._react_regions[rptr]:
                 self._react_regions[rptr].append(region)

test/rxd/test_multicompartment_reactions_aligment.py

@@ -0,0 +1,70 @@
+from testutils import compare_data, tol
+
+def test_multicompartment_reactions_aligment(neuron_instance):
+    """ A test for multicompartment reactions where one regions has more
+    sections than the other. 
+    """
+
+    h, rxd, data, save_path = neuron_instance
+
+    #Dendrite and spine
+    dend = h.Section(name='dend')
+    dend.nseg = 101
+    dend.pt3dclear()
+    dend.pt3dadd(-10, 0, 0, 1)
+    dend.pt3dadd(10, 0, 0, 1)
+
+    spine= h.Section(name='spine')
+    spine.nseg = 11
+    spine.pt3dclear()
+    spine.pt3dadd(0, 0, 0, 0.2)
+    spine.pt3dadd(0, 1, 0, 0.2)
+    spine.connect(dend(0.5))
+
+
+    # Where -- define the shells and borders
+    N = 5                   # number of shells -- must be >= 2
+
+    # Where -- shells and border between them
+    shells = []
+    border = []
+    for i in range(N-1):
+        shells.append(rxd.Region([dend], name='shell%i' % i,
+                      geometry=rxd.Shell(float(i)/N, (1.0+i)/N)))
+        border.append(rxd.Region([dend], name='border%i' % i,
+                      geometry=rxd.FixedPerimeter(2.0*h.PI*(1.0+i))))
+
+    shells.append(rxd.Region([dend,spine], nrn_region='i', name='shell%i' % (N-1),
+                  geometry=rxd.MultipleGeometry(secs=[dend, spine],
+                                                geos=[rxd.Shell((N-1.0)/N,1),
+                                                      rxd.inside])))
+
+    # Who -- calcium with an inhomogeneous initial condition
+    Dca = 0.6   #um^2/ms
+
+    # scale factor so the flux (Dca/dr)*Ca has units molecules/um^2/ms (where dr is the thickness of the shell)
+    mM_to_mol_per_um = 6.0221409e+23 * 1e-18
+    ca = rxd.Species(shells, d=Dca, name='ca', charge=2, initial=lambda nd: 1000e-6 if nd.segment in spine else 100e-6)
+
+
+    # What -- use reactions to setup diffusion between the shells
+    cas = []    # calcium on the shells
+    for reg in shells:
+        cas.append(ca[reg])
+
+    # create the multi-compartment reactions between the pairs of shells
+    diffusions = []
+    for i in range(N-1):
+        diffusions.append(rxd.MultiCompartmentReaction(cas[i], cas[i+1], 
+                          mM_to_mol_per_um*Dca, 
+                          mM_to_mol_per_um*Dca,
+                          border=border[i]))
+
+
+    h.finitialize(-65)
+    h.continuerun(2.5)
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol
+
+