diff --git a/pyvcell/_internal/simdata/zarr_writer.py b/pyvcell/_internal/simdata/zarr_writer.py index 090cc4d..3f8bc93 100644 --- a/pyvcell/_internal/simdata/zarr_writer.py +++ b/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 + # ``::``) 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) diff --git a/pyvcell/vcml/__init__.py b/pyvcell/vcml/__init__.py index dda84f9..c4b9148 100644 --- a/pyvcell/vcml/__init__.py +++ b/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", diff --git a/pyvcell/vcml/models_app.py b/pyvcell/vcml/models_app.py index 41953b8..1802d30 100644 --- a/pyvcell/vcml/models_app.py +++ b/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]: diff --git a/pyvcell/vcml/models_math.py b/pyvcell/vcml/models_math.py index 6e2badd..34664c3 100644 --- a/pyvcell/vcml/models_math.py +++ b/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 + ```` ("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) diff --git a/pyvcell/vcml/vcml_reader.py b/pyvcell/vcml/vcml_reader.py index 3d24d29..549badd 100644 --- a/pyvcell/vcml/vcml_reader.py +++ b/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 ```` 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) diff --git a/pyvcell/vcml/vcml_writer.py b/pyvcell/vcml/vcml_writer.py index 4423620..b1f06f1 100644 --- a/pyvcell/vcml/vcml_writer.py +++ b/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``; ```` 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) diff --git a/tests/conftest.py b/tests/conftest.py index 4ba4c4e..38587f5 100644 --- a/tests/conftest.py +++ b/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, ) diff --git a/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml b/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml new file mode 100644 index 0000000..554a68a --- /dev/null +++ b/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml @@ -0,0 +1,311 @@ + + + + + + + 1.0 + 0.5 + 1.0 + 1.0 + 1.0 + 1.0 + + + s0 + + + s1 + + + s2 + + + s3 + + + + + + + + + + + + + ((Kf_r0 * s0) - (1.0 * Kr_r0 * s1)) + + + + + + + ((Kf_r1 * s0) - (Kr_r1 * s2)) + 0.0 + 1.0 + + + + + + + ((Kf_r2 * s2) - (Kr_r2 * s3)) + 0.0 + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ((pow(( - 5.0 + x),2.0) + pow(( - 5.0 + y),2.0) + pow(( - 5.0 + z),2.0)) < 16.0) + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + (1.0 + sin(x)) + 1.0E-4 + + + + (1.0 + cos(x)) + 1.0E-4 + + + + (1.0 + (sin(x) * cos(y))) + 1.0E-4 + + + + (1.0 + cos(y)) + 1.0000000000000002E-6 + + + + + + + + 96485.3321 + 9.64853321E-5 + 1.0E-9 + 6.02214179E11 + 3.141592653589793 + 8314.46261815 + 300.0 + 1.0 + 1000.0 + 1.0 + 1.0 + 1.0 + 0.001660538783162726 + 0.5 + 1.0 + 1.0 + 1.0 + 1.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 1.0E-4 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 1.0E-4 + 0.0 + 1.0 + 1.0 + + + + + ((Kf_r0 * s0) - (1.0 * Kr_r0 * s1)) + ((Kf_r1 * s0) - (Kr_r1 * s2)) + ((Kf_r2 * s2) - (Kr_r2 * s3)) + (AreaPerUnitArea_m0 / VolumePerUnitVolume_c0) + (AreaPerUnitArea_m0 / VolumePerUnitVolume_c1) + (1.0 + sin(x)) + (1.0 + cos(x)) + (1.0 + cos(y)) + (1.0 + (sin(x) * cos(y))) + (VolumePerUnitVolume_c0 * vcRegionVolume('subdomain0')) + (VolumePerUnitVolume_c1 * vcRegionVolume('subdomain1')) + (AreaPerUnitArea_m0 * vcRegionArea('subdomain0_subdomain1_membrane')) + vcRegionArea('subdomain0_subdomain1_membrane') + vcRegionVolume('subdomain0') + vcRegionVolume('subdomain1') + + + + + + + + + + 0.0 + s3_diffusionRate + s3_init_umol_l_1 + + + + + + + + + + + + - J_r0 + s0_diffusionRate + s0_init_umol_l_1 + + + 0.0 + J_r0 + s1_init_umol_l_1 + + + + + + + + + + + + 0.0 + (J_r1 - J_r2) + s2_init_umol_dm_2 + + + + ( - 1.0 * KFlux_m0_c0 * J_r1) + 0.0 + + + 0.0 + 0.0 + + + 0.0 + (KFlux_m0_c1 * J_r2) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/fixtures/vcell_model_fixtures.py b/tests/fixtures/vcell_model_fixtures.py index d6e063c..526fe29 100644 --- a/tests/fixtures/vcell_model_fixtures.py +++ b/tests/fixtures/vcell_model_fixtures.py @@ -19,6 +19,13 @@ def vcml_spatial_small_3d_path() -> Path: return FIXTURE_DATA_DIR / "SmallSpatialProject_3D.vcml" +@pytest.fixture +def vcml_spatial_well_mixed_3d_path() -> Path: + """SmallSpatialProject_3D with s1 (volume) and s2 (membrane) well-mixed, math regenerated by + libvcell: s1 becomes a VolumeRegionVariable, s2 a MembraneRegionVariable.""" + return FIXTURE_DATA_DIR / "SmallSpatialProject_3D_well_mixed.vcml" + + @pytest.fixture def vcml_spatial_bunny_3d_path() -> Path: return FIXTURE_DATA_DIR / "Bunny.vcml" diff --git a/tests/sim_results/test_region_variables.py b/tests/sim_results/test_region_variables.py new file mode 100644 index 0000000..b7d3973 --- /dev/null +++ b/tests/sim_results/test_region_variables.py @@ -0,0 +1,71 @@ +"""Well-mixed species in a spatial application: region variables through math, VCML and results. + +A well-mixed species has one value per connected region. VCell generates it as a VolumeRegionVariable / +MembraneRegionVariable with a region equation, and fvsolver writes it as one value per region. The zarr +writer replicates a volume-region variable onto every voxel of its region, so it reads like a volume +variable: a channel of its own, and a binding for the functions that use it. +""" + +from pathlib import Path + +import numpy as np + +import pyvcell.vcml as vc +from pyvcell._internal.simdata.simdata_models import VariableType + + +def test_species_mapping_flags_roundtrip(vcml_spatial_well_mixed_3d_path: Path) -> None: + """WellMixed / ForceConstant / ForceContinuous survive a read -> write -> read (they used to be dropped, + so a well-mixed species came back spatial, and a clamped species came back unclamped).""" + biomodel = vc.VcmlReader.biomodel_from_file(vcml_spatial_well_mixed_3d_path) + flags = {m.species_name: m.well_mixed for m in biomodel.applications[0].species_mappings} + assert flags == {"s0": False, "s1": True, "s2": True, "s3": False} + + mapping = next(m for m in biomodel.applications[0].species_mappings if m.species_name == "s0") + mapping.force_constant = True + vcml_str = vc.VcmlWriter().write_vcml(document=vc.VCMLDocument(biomodel=biomodel)) + assert 'LocalizedCompoundRef="s1" ForceConstant="false" WellMixed="true"' in vcml_str + roundtripped = vc.VcmlReader.biomodel_from_str(vcml_str) + assert roundtripped == biomodel + + +def test_volume_region_variable_in_results(vcml_spatial_well_mixed_3d_path: Path) -> None: + biomodel = vc.load_vcml_file(vcml_spatial_well_mixed_3d_path) + biomodel.applications[0].add_sim(name="sim", duration=1.0, output_time_step=0.25, mesh_size=(12, 12, 12)) + result = vc.simulate(biomodel, "sim") + + labels = [channel.label for channel in result.channel_data] + # s1 (well-mixed, subdomain0) is a channel after the volume variables; J_r0 = f(s0, s1) now evaluates + assert labels == [ + "region_mask", + "t", + "x", + "y", + "z", + "s0", + "s3", + "s1", + "J_r0", + "s0_init_uM", + "s1_init_uM", + "s3_init_uM", + ] + + header = next(h for h in result.pde_dataset.variables_block_headers() if h.var_info.var_name == "subdomain0::s1") + assert header.var_info.variable_type == VariableType.VOLUME_REGION + + zarr_data = result.zarr_dataset + region_map = np.asarray(zarr_data[0, 0]) + in_subdomain0 = np.isin(region_map, list(result.mesh.get_volume_region_ids("subdomain0"))) + s1_index = labels.index("s1") + times = result.pde_dataset.times() + for t_index, t in enumerate(times): + s1 = np.asarray(zarr_data[t_index, s1_index]) + per_region = np.asarray(result.pde_dataset.get_data(header.var_info, t)) + # every voxel holds its region's value: constant across subdomain0, equal to the solver's value + assert np.allclose(s1, per_region[region_map.astype(int)]) + assert np.ptp(s1[in_subdomain0]) == 0.0 + # the species actually evolves (so the replication is not trivially constant in time) + assert not np.isclose( + np.asarray(zarr_data[0, s1_index])[in_subdomain0][0], np.asarray(zarr_data[-1, s1_index])[in_subdomain0][0] + ) diff --git a/tests/vcml/test_math_description.py b/tests/vcml/test_math_description.py index 402f848..2392082 100644 --- a/tests/vcml/test_math_description.py +++ b/tests/vcml/test_math_description.py @@ -57,6 +57,35 @@ def test_math_description_membrane_jump_condition(vcml_spatial_bunny_3d_path: Pa assert jump_condition.out_flux is not None +def test_math_description_region_equations(vcml_spatial_well_mixed_3d_path: Path) -> None: + """Well-mixed species in a spatial application become region variables with region equations.""" + biomodel = vc.VcmlReader.biomodel_from_file(vcml_spatial_well_mixed_3d_path) + math = biomodel.applications[0].math_description + assert math is not None + var_types = {v.name: v.var_type for v in math.variables} + assert var_types["s1"] == vc.MathVariableType.volume_region + assert var_types["s2"] == vc.MathVariableType.membrane_region + + (subdomain0,) = (c for c in math.compartment_subdomains if c.name == "subdomain0") + assert subdomain0.volume_region_equations == [ + vc.VolumeRegionEquation( + name="s1", uniform_rate="0.0", volume_rate="J_r0", initial="s1_init_umol_l_1", solution_type="unknown" + ) + ] + assert "s1" not in [eq.name for eq in subdomain0.pde_equations + subdomain0.ode_equations] + + (membrane,) = math.membrane_subdomains + assert membrane.membrane_region_equations == [ + vc.MembraneRegionEquation( + name="s2", + uniform_rate="0.0", + membrane_rate="(J_r1 - J_r2)", + initial="s2_init_umol_dm_2", + solution_type="unknown", + ) + ] + + def test_math_description_ode_nonspatial(vcml_nonspatial_ode_path: Path) -> None: biomodel = vc.VcmlReader.biomodel_from_file(vcml_nonspatial_ode_path) math = biomodel.applications[0].math_description @@ -115,10 +144,11 @@ def test_math_description_roundtrip( vcml_spatial_particle_path: Path, vcml_spatial_model_1d_path: Path, vcml_spatial_bunny_3d_path: Path, + vcml_spatial_well_mixed_3d_path: Path, ) -> None: """Full biomodel write -> read preserves the math description (and everything else). - Covers nonspatial (ODE / stochastic) and spatial (PDE / membrane / particle) + Covers nonspatial (ODE / stochastic) and spatial (PDE / membrane / particle / region) geometries. """ paths = [ @@ -127,6 +157,7 @@ def test_math_description_roundtrip( vcml_spatial_particle_path, vcml_spatial_model_1d_path, vcml_spatial_bunny_3d_path, + vcml_spatial_well_mixed_3d_path, ] for path in paths: biomodel = vc.VcmlReader.biomodel_from_file(path)