Fix RxD MultiCompartment Reactions with different number of source and destination terms.

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

@@ -379,27 +379,31 @@ def intersection(los):
                 ] + [areas / volumes[di] / molecules_per_mM_um3 for di in dests_indices]
             # TODO: check for multicompartment reaction within the ECS
             elif sources_ecs and not dests_ecs:
+                # only using the first dests_indices for scaling ECS sources
+                # each dest must be in same location so give the same alpha
                 self._mult = [
                     -areas
                     / (
                         numpy.prod(s()._extracellular()._dx)
-                        * s().alpha_by_location(xyz_by_index(di))
+                        * s().alpha_by_location(xyz_by_index(dests_indices[0]))
                     )
                     / molecules_per_mM_um3
-                    for s, di in zip(sources_ecs, dests_indices)
+                    for s in sources_ecs
                 ] + [areas / volumes[di] / molecules_per_mM_um3 for di in dests_indices]
             elif not sources_ecs and dests_ecs:
+                # only using the first sources_indices for scaling ECS dests
+                # each source must be in same location so give the same alpha
                 self._mult = [
                     -areas / volumes[si] / molecules_per_mM_um3
                     for si in sources_indices
                 ] + [
                     areas
                     / (
                         numpy.prod(s()._extracellular()._dx)
-                        * s().alpha_by_location(xyz_by_index(si))
+                        * s().alpha_by_location(xyz_by_index(sources_indices[0]))
                     )
                     / molecules_per_mM_um3
-                    for s, si in zip(dests_ecs, sources_indices)
+                    for s in dests_ecs
                 ]
             elif self._trans_membrane:
                 # An ecs <-> ecs reaction that use the membrane area and intracellular concentration in the rates

test/rxd/test_multicompartment_reactions.py

@@ -89,3 +89,101 @@ def test_multicompartment_reactions(neuron_instance):
     if not save_path:
         max_err = compare_data(data)
         assert max_err < tol
+
+
+def test_mcr_multiple_dests(neuron_instance):
+    """Test that a MultiCompartmentReaction with destinations does not corrupt
+    the multiplier array for other reactions."""
+    h, rxd, data, save_path = neuron_instance
+
+    dend = h.Section("dend")
+    cyt = rxd.Region([dend], name="cyt", nrn_region="i")
+    mem = rxd.Region([dend], name="cell_mem", geometry=rxd.membrane())
+    ecs = rxd.Extracellular(-100, -100, -100, 100, 100, 100, dx=100)
+
+    a = rxd.Species(
+        [cyt, mem, ecs],
+        name="a",
+        charge=1,
+        d=0,
+        initial=10.0,
+        ecs_boundary_conditions=140.0,
+    )
+    # Extra species for the MCR with multiple ECS destinations
+    x = rxd.Species([cyt, ecs], name="x", charge=1, d=0, initial=1.0)
+    y = rxd.Species([cyt, ecs], name="y", charge=1, d=0, initial=2.0)
+    z = rxd.Species([cyt, ecs], name="z", charge=1, d=0, initial=3.0)
+
+    mcr = rxd.MultiCompartmentReaction(
+        x[cyt] + y[cyt],
+        y[ecs] + y[ecs] + z[ecs],
+        10,
+        mass_action=False,
+        membrane=mem,
+        membrane_flux=False,
+    )
+
+    gl = 0.0003 * 1e-14 / 1.60217662e-19  # same scale as HH leak
+    leak = rxd.MultiCompartmentReaction(
+        a[cyt],
+        a[ecs],
+        gl * (rxd.v - (-65.0)),
+        mass_action=False,
+        membrane=mem,
+        membrane_flux=True,
+    )
+
+    h.finitialize(-70)
+    h.continuerun(100)
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol
+
+
+def test_mcr_multiple_sources(neuron_instance):
+    """Test that a MultiCompartmentReaction with sources does not corrupt the
+    multiplier array for other reactions."""
+
+    h, rxd, data, save_path = neuron_instance
+    dend = h.Section("dend")
+    cyt = rxd.Region([dend], name="cyt", nrn_region="i")
+    mem = rxd.Region([dend], name="cell_mem", geometry=rxd.membrane())
+    ecs = rxd.Extracellular(-100, -100, -100, 100, 100, 100, dx=100)
+
+    a = rxd.Species(
+        [cyt, mem, ecs],
+        name="a",
+        charge=1,
+        d=0,
+        initial=10.0,
+        ecs_boundary_conditions=140.0,
+    )
+    # Extra species for the MCR with multiple ECS sources
+    x = rxd.Species([cyt, ecs], name="x", charge=1, d=0, initial=1.0)
+    y = rxd.Species([cyt, ecs], name="y", charge=1, d=0, initial=2.0)
+    z = rxd.Species([cyt, ecs], name="z", charge=1, d=0, initial=3.0)
+
+    mcr = rxd.MultiCompartmentReaction(
+        x[cyt] + y[cyt] + z[cyt],
+        x[ecs] + y[ecs],
+        10,
+        mass_action=False,
+        membrane=mem,
+        membrane_flux=False,
+    )
+
+    gl = 0.0003 * 1e-14 / 1.60217662e-19  # same scale as HH leak
+    leak = rxd.MultiCompartmentReaction(
+        a[cyt],
+        a[ecs],
+        gl * (rxd.v - (-65.0)),
+        mass_action=False,
+        membrane=mem,
+        membrane_flux=True,
+    )
+
+    h.finitialize(-70)
+    h.continuerun(100)
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol