Region variables: VolumeRegion/MembraneRegion equations in the math model, VCML, and results

pyvcell/_internal/simdata/zarr_writer.py

@@ -20,10 +20,20 @@ def write_zarr(pde_dataset: PdeDataSet, data_functions: DataFunctions, mesh: Car
     volume_data_vars: list[DataBlockHeader] = [
         v for v in pde_dataset.variables_block_headers() if v.var_info.variable_type == VariableType.VOLUME
     ]
+    # Volume-region state variables (well-mixed species: one value per connected volume region, as
+    # ``<domain>::<name>``) are replicated onto every voxel of their region, so they read like volume
+    # variables — a channel each, and a numexpr binding for the functions that use them. The solver's
+    # built-in region quantities (``vcRegionVolume``, no domain) are not state and stay out.
+    volume_region_vars: list[DataBlockHeader] = [
+        v
+        for v in pde_dataset.variables_block_headers()
+        if v.var_info.variable_type == VariableType.VOLUME_REGION and "::" in v.var_info.var_name
+    ]
     volume_functions: list[NamedFunction] = [
         f for f in data_functions.named_functions if f.variable_type == VariableType.VOLUME
     ]
-    num_channels = len(volume_data_vars) + len(volume_functions) + 5  # 5 extra channels for region map, t, x, y, z
+    # 5 extra channels for region map, t, x, y, z
+    num_channels = len(volume_data_vars) + len(volume_region_vars) + len(volume_functions) + 5
     num_t: int = len(pde_dataset.times())
     times: list[float] = pde_dataset.times()
     header = pde_dataset.first_data_zip_file_metadata().file_header
@@ -142,6 +152,28 @@ def _get_domain_mask(dname: str) -> np.ndarray:
             channel_metadata[c]["mean_values"].append(float(np.mean(masked)))
             c = c + 1
 
+        # add volume-region state variables, replicated per region onto the voxels
+        for v in volume_region_vars:
+            region_values: NDArray[np.float64] = np.asarray(pde_dataset.get_data(v.var_info, times[t]))
+            var_data = region_values[region_map_flat].reshape((num_z, num_y, num_x))
+            z1[t, c, :, :, :] = var_data
+            domain_name, var_name = v.var_info.var_name.split("::")
+            bindings[var_name] = var_data
+            if t == 0:
+                channel_metadata.append({
+                    "index": c,
+                    "label": var_name,
+                    "domain_name": domain_name,
+                    "min_values": [],
+                    "max_values": [],
+                    "mean_values": [],
+                })
+            masked = var_data[_get_domain_mask(domain_name)]
+            channel_metadata[c]["min_values"].append(float(np.min(masked)))
+            channel_metadata[c]["max_values"].append(float(np.max(masked)))
+            channel_metadata[c]["mean_values"].append(float(np.mean(masked)))
+            c = c + 1
+
         # add volumetric functions
         for f in volume_functions:
             evaluated_data = f.evaluate(variable_bindings=bindings)

pyvcell/vcml/__init__.py

@@ -56,6 +56,7 @@
     MathFunction,
     MathVariable,
     MathVariableType,
+    MembraneRegionEquation,
     MembraneSubDomain,
     OdeEquation,
     ParticleInitialCount,
@@ -64,6 +65,7 @@
     PdeEquation,
     VariableInitialCount,
     Velocity,
+    VolumeRegionEquation,
 )
 from pyvcell.vcml.vcml_reader import VcmlReader
 
@@ -211,6 +213,7 @@ def __dir__() -> list[str]:
     "MathFunction",
     "MathVariable",
     "MathVariableType",
+    "MembraneRegionEquation",
     "MembraneSubDomain",
     "Model",
     "ModelParameter",
@@ -240,6 +243,7 @@ def __dir__() -> list[str]:
     "VcmlWriter",
     "Velocity",
     "Version",
+    "VolumeRegionEquation",
     "cartesian_mesh_from_geometry",
     "connect",
     "field_data_refs",

pyvcell/vcml/models_app.py

@@ -97,6 +97,13 @@ class SpeciesMapping(VcmlNode):
     velocity_y: float | str | None = None
     velocity_z: float | str | None = None
     boundary_values: list[float | str | None] = Field(default_factory=list)
+    # VCell's LocalizedCompoundSpec flags. ``force_constant``: a clamped species (its initial condition for
+    # all time). ``well_mixed``: in a spatial application, one value per region — the math becomes a
+    # VolumeRegionVariable / MembraneRegionVariable. ``force_continuous``: keep a species deterministic in a
+    # hybrid stochastic application.
+    force_constant: bool = False
+    well_mixed: bool = False
+    force_continuous: bool = False
 
     @property
     def expressions(self) -> list[str]:

pyvcell/vcml/models_math.py

@@ -108,6 +108,38 @@ class PdeEquation(VcmlNode):
     velocity: Velocity | None = None
 
 
+class VolumeRegionEquation(VcmlNode):
+    """An ODE for a volume-region variable: one value per connected region of the compartment.
+
+    ``d/dt u_R = uniform_rate + (1/|R|) ∫_R volume_rate dV`` over each region ``R`` (plus the fluxes the
+    membranes bounding ``R`` contribute through their jump conditions). ``solution_type`` is VCell's
+    ``<SolutionType Type=...>`` ("unknown" or "exact"); an exact solution is ``solution``.
+    """
+
+    name: str
+    uniform_rate: str | None = None
+    volume_rate: str | None = None
+    initial: str | None = None
+    solution_type: str | None = None
+    solution: str | None = None
+
+
+class MembraneRegionEquation(VcmlNode):
+    """An ODE for a membrane-region variable: one value per connected region of the membrane.
+
+    ``d/dt u_R = uniform_rate + (1/|R|) ∫_R membrane_rate dA`` over each membrane region ``R`` — most often
+    the membrane potential of an electrically coupled model. ``solution_type`` / ``solution`` as for
+    :class:`VolumeRegionEquation`.
+    """
+
+    name: str
+    uniform_rate: str | None = None
+    membrane_rate: str | None = None
+    initial: str | None = None
+    solution_type: str | None = None
+    solution: str | None = None
+
+
 class JumpCondition(VcmlNode):
     """A membrane flux jump condition for a volume variable."""
 
@@ -181,6 +213,7 @@ class CompartmentSubDomain(VcmlNode):
     boundary_types: list[MathBoundaryType] = Field(default_factory=list)
     ode_equations: list[OdeEquation] = Field(default_factory=list)
     pde_equations: list[PdeEquation] = Field(default_factory=list)
+    volume_region_equations: list[VolumeRegionEquation] = Field(default_factory=list)
     variable_initial_counts: list[VariableInitialCount] = Field(default_factory=list)
     jump_processes: list[JumpProcess] = Field(default_factory=list)
     particle_jump_processes: list[ParticleJumpProcess] = Field(default_factory=list)
@@ -196,6 +229,7 @@ class MembraneSubDomain(VcmlNode):
     boundary_types: list[MathBoundaryType] = Field(default_factory=list)
     ode_equations: list[OdeEquation] = Field(default_factory=list)
     pde_equations: list[PdeEquation] = Field(default_factory=list)
+    membrane_region_equations: list[MembraneRegionEquation] = Field(default_factory=list)
     jump_conditions: list[JumpCondition] = Field(default_factory=list)
     particle_jump_processes: list[ParticleJumpProcess] = Field(default_factory=list)
     particle_properties: list[ParticleProperties] = Field(default_factory=list)

pyvcell/vcml/vcml_reader.py

@@ -339,6 +339,12 @@ def _parse_compartment_subdomain(self, element: _Element) -> vcm.CompartmentSubD
                 subdomain.ode_equations.append(self._parse_ode_equation(child))
             elif tag == "PdeEquation":
                 subdomain.pde_equations.append(self._parse_pde_equation(child))
+            elif tag == "VolumeRegionEquation":
+                subdomain.volume_region_equations.append(
+                    vcm.VolumeRegionEquation.model_validate(
+                        self._parse_region_equation(child, "VolumeRate", "volume_rate")
+                    )
+                )
             elif tag in ("VariableInitialCount", "VariableInitialPoissonExpectedCount"):
                 subdomain.variable_initial_counts.append(
                     vcm.VariableInitialCount(
@@ -368,6 +374,12 @@ def _parse_membrane_subdomain(self, element: _Element) -> vcm.MembraneSubDomain:
                 subdomain.ode_equations.append(self._parse_ode_equation(child))
             elif tag == "PdeEquation":
                 subdomain.pde_equations.append(self._parse_pde_equation(child))
+            elif tag == "MembraneRegionEquation":
+                subdomain.membrane_region_equations.append(
+                    vcm.MembraneRegionEquation.model_validate(
+                        self._parse_region_equation(child, "MembraneRate", "membrane_rate")
+                    )
+                )
             elif tag == "JumpCondition":
                 subdomain.jump_conditions.append(self._parse_jump_condition(child))
             elif tag in ("ParticleJumpProcess", "LangevinParticleJumpProcess"):
@@ -376,6 +388,25 @@ def _parse_membrane_subdomain(self, element: _Element) -> vcm.MembraneSubDomain:
                 subdomain.particle_properties.append(self._parse_particle_properties(child))
         return subdomain
 
+    def _parse_region_equation(self, element: _Element, rate_tag: str, rate_field: str) -> dict[str, str | None]:
+        """The fields of a ``VolumeRegionEquation`` / ``MembraneRegionEquation``: ``UniformRate``, the
+        region's own rate (``rate_tag``), ``Initial``, and ``<SolutionType Type="unknown|exact">`` whose
+        text is the exact solution."""
+        fields: dict[str, str | None] = {"name": element.get("Name", default="unnamed")}
+        for child in element:
+            tag = strip_namespace(child.tag)
+            if tag == "UniformRate":
+                fields["uniform_rate"] = _text(child)
+            elif tag == rate_tag:
+                fields[rate_field] = _text(child)
+            elif tag == "Initial":
+                fields["initial"] = _text(child)
+            elif tag == "SolutionType":
+                fields["solution_type"] = child.get("Type")
+                if child.get("Type") == "exact":
+                    fields["solution"] = _text(child)
+        return fields
+
     def _parse_ode_equation(self, element: _Element) -> vcm.OdeEquation:
         equation = vcm.OdeEquation(
             name=element.get("Name", default="unnamed"), solution_type=element.get("SolutionType")
@@ -613,7 +644,12 @@ def visit_BoundariesTypes(self, element: _Element, node: vc.CompartmentMapping)
 
     def visit_LocalizedCompoundSpec(self, element: _Element, node: vc.Application) -> None:
         species_name: str = element.get("LocalizedCompoundRef", default="unnamed")
-        species_mapping = vc.SpeciesMapping(species_name=species_name)
+        species_mapping = vc.SpeciesMapping(
+            species_name=species_name,
+            force_constant=element.get("ForceConstant") == "true",
+            well_mixed=element.get("WellMixed") == "true",
+            force_continuous=element.get("ForceContinuous") == "true",
+        )
         node.species_mappings.append(species_mapping)
         self.generic_visit_children(element, species_mapping)
 

pyvcell/vcml/vcml_writer.py

@@ -26,11 +26,13 @@
     JumpProcess,
     MathBoundaryType,
     MathDescription,
+    MembraneRegionEquation,
     MembraneSubDomain,
     OdeEquation,
     ParticleJumpProcess,
     ParticleProperties,
     PdeEquation,
+    VolumeRegionEquation,
 )
 
 
@@ -218,7 +220,13 @@ def write_application(self, application: Application, parent: _Element) -> None:
         reaction_context_element = Element("ReactionContext")
         parent.append(reaction_context_element)
         for species_mapping in application.species_mappings:
-            mapping_element = Element("LocalizedCompoundSpec", LocalizedCompoundRef=species_mapping.species_name)
+            mapping_element = Element(
+                "LocalizedCompoundSpec",
+                LocalizedCompoundRef=species_mapping.species_name,
+                ForceConstant=str(species_mapping.force_constant).lower(),
+                WellMixed=str(species_mapping.well_mixed).lower(),
+                ForceContinuous=str(species_mapping.force_continuous).lower(),
+            )
             reaction_context_element.append(mapping_element)
             self.write_species_mapping(species_mapping, mapping_element)
         for reaction_mapping in application.reaction_mappings:
@@ -466,6 +474,14 @@ def write_compartment_subdomain(self, subdomain: CompartmentSubDomain, parent: _
             self.write_ode_equation(ode_equation, subdomain_element)
         for pde_equation in subdomain.pde_equations:
             self.write_pde_equation(pde_equation, subdomain_element)
+        for volume_region_equation in subdomain.volume_region_equations:
+            self.write_region_equation(
+                "VolumeRegionEquation",
+                volume_region_equation,
+                "VolumeRate",
+                volume_region_equation.volume_rate,
+                subdomain_element,
+            )
         for variable_initial_count in subdomain.variable_initial_counts:
             tag = "VariableInitialPoissonExpectedCount" if variable_initial_count.poisson else "VariableInitialCount"
             count_element = Element(tag, Name=variable_initial_count.name)
@@ -491,6 +507,14 @@ def write_membrane_subdomain(self, subdomain: MembraneSubDomain, parent: _Elemen
             self.write_ode_equation(ode_equation, subdomain_element)
         for pde_equation in subdomain.pde_equations:
             self.write_pde_equation(pde_equation, subdomain_element)
+        for membrane_region_equation in subdomain.membrane_region_equations:
+            self.write_region_equation(
+                "MembraneRegionEquation",
+                membrane_region_equation,
+                "MembraneRate",
+                membrane_region_equation.membrane_rate,
+                subdomain_element,
+            )
         for jump_condition in subdomain.jump_conditions:
             self.write_jump_condition(jump_condition, subdomain_element)
         for particle_jump_process in subdomain.particle_jump_processes:
@@ -542,6 +566,26 @@ def write_pde_equation(self, equation: PdeEquation, parent: _Element) -> None:
             equation_element.append(Element("Velocity", attrib=velocity_attrs))
         self._append_text_element(equation_element, "Solution", equation.solution)
 
+    def write_region_equation(
+        self,
+        tag: str,
+        equation: VolumeRegionEquation | MembraneRegionEquation,
+        rate_tag: str,
+        rate: str | None,
+        parent: _Element,
+    ) -> None:
+        """A ``VolumeRegionEquation`` / ``MembraneRegionEquation`` as VCell's Xmlproducer writes it: the
+        rates default to ``0.0``; ``<SolutionType Type=...>`` carries an exact solution as its text."""
+        equation_element = Element(tag, Name=equation.name)
+        parent.append(equation_element)
+        self._append_text_element(equation_element, "UniformRate", equation.uniform_rate or "0.0")
+        self._append_text_element(equation_element, rate_tag, rate or "0.0")
+        self._append_text_element(equation_element, "Initial", equation.initial)
+        solution_type = Element("SolutionType", Type=equation.solution_type or "unknown")
+        if solution_type.get("Type") == "exact" and equation.solution is not None:
+            solution_type.text = equation.solution
+        equation_element.append(solution_type)
+
     def write_jump_condition(self, condition: JumpCondition, parent: _Element) -> None:
         condition_element = Element("JumpCondition", Name=condition.name)
         parent.append(condition_element)

tests/conftest.py

@@ -29,6 +29,7 @@
     vcml_spatial_model_1d_path,
     vcml_spatial_particle_path,
     vcml_spatial_small_3d_path,
+    vcml_spatial_well_mixed_3d_path,
     vcml_tutorial_multiapp_pde_path,
 )