From 0121b75628d26717321b63a890bad5d2dd482c2a Mon Sep 17 00:00:00 2001 From: Eric Cramer <13970720+emcramer@users.noreply.github.com> Date: Thu, 23 Jul 2026 14:37:36 -0700 Subject: [PATCH] feat: PhysiCell substrate and microenvironment access MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Close the main functional gap from the v0.2.0 defect report: the microenvironment .mat was parseable but unreachable from the high-level reader classes, so substrate-coupled analysis — the common case for tumour ABM work — could not be expressed through the API. Add two complementary views of "how much substrate": - Environmental concentration at a location. PhysiCellTimeStep gains substrates, voxel_positions, and substrate_at(name, x, y, z=None), which returns the concentration in the nearest voxel via a cached cKDTree. Verified exactly equal to a brute-force argmin nearest-voxel lookup on every cell of the example frame. - Per-cell internalized amount. internalized_substrates() returns PhysiCell's internalized_total_substrates field as a DataFrame keyed by substrate name, reflecting the framework's actual uptake dynamics rather than approximating uptake from the voxel the cell sits in. The microenvironment file is located from the XML's microenvironment/domain/data/filename element (falling back to the default naming), resolved lazily inside the timestep. discover_physicell_ timesteps keeps its 3-tuple return type, so no existing caller breaks. Adds tests/test_physicell_microenvironment.py (22 tests). --- CHANGELOG.md | 43 +++ spatialtissuepy/synthetic/physicell/parser.py | 12 + spatialtissuepy/synthetic/physicell/reader.py | 258 +++++++++++++++++- tests/test_physicell_microenvironment.py | 237 ++++++++++++++++ 4 files changed, 549 insertions(+), 1 deletion(-) create mode 100644 tests/test_physicell_microenvironment.py diff --git a/CHANGELOG.md b/CHANGELOG.md index b215a03..c9bb2cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,49 @@ All notable changes to spatialtissuepy are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +Substrate / microenvironment access for PhysiCell simulations — the main +functional gap called out in the v0.2.0 defect report. Additive; nothing +existing changes behavior. + +### Added + +- **Environmental substrate sampling on `PhysiCellTimeStep`.** The + microenvironment `.mat` (oxygen and other diffusible fields) is now reachable + through the high-level API, having previously been parseable only by calling + `parse_microenvironment_mat` directly: + - `substrates` → `{name: (n_voxels,) array}` of concentration fields + - `voxel_positions` → `(n_voxels, 3)` voxel-center coordinates + - `substrate_at(name, x, y, z=None)` → concentration in the voxel nearest each + query point, via a cached `scipy.spatial.cKDTree`. Verified to match a + brute-force `argmin` nearest-voxel lookup exactly on every cell of the + bundled example frame. + - `substrate_names` on both `PhysiCellTimeStep` and `PhysiCellSimulation`. +- **Per-cell internalized substrates.** `PhysiCellTimeStep.internalized_substrates()` + returns a DataFrame of PhysiCell's `internalized_total_substrates` field, one + column per substrate (named), one row per cell, in `positions` order. This is + the amount each cell has actually taken up — accounting for PhysiCell's uptake + dynamics — as opposed to the environmental concentration at the cell's + location that `substrate_at` gives. (Populated only when the simulation + enabled `track_internalized_substrates_in_each_agent`; otherwise present but + zero, which the method documents.) +- `parse_physicell_xml` now records the microenvironment field-data filename at + `metadata.extra['microenvironment_file']`, read from the XML's + `microenvironment/domain/data/filename` element so relocated or renamed + outputs still resolve. The reader falls back to the default + `output{index:08d}_microenvironment0.mat` naming. +- Regression tests (`tests/test_physicell_microenvironment.py`, 22 tests) + covering nearest-voxel equivalence, tie handling, 2-D/3-D and scalar/array + queries, the internalized-vs-raw-matrix mapping, empty-microenvironment + frames, file resolution, and `discover_physicell_timesteps` arity. + +### Notes + +- `discover_physicell_timesteps` keeps its 3-tuple return type unchanged; the + microenvironment file is resolved lazily inside `PhysiCellTimeStep`, so no + caller that unpacks the tuple is affected. + ## [0.3.0] - 2026-07-23 PhysiCell reader correctness fixes, reported against v0.2.0 by the diff --git a/spatialtissuepy/synthetic/physicell/parser.py b/spatialtissuepy/synthetic/physicell/parser.py index 59fc12f..58e9e85 100644 --- a/spatialtissuepy/synthetic/physicell/parser.py +++ b/spatialtissuepy/synthetic/physicell/parser.py @@ -134,6 +134,17 @@ def parse_physicell_xml(xml_path: Path) -> PhysiCellMetadata: name = var.get('name', 'unknown') substrate_names.append(name) + # Extract the microenvironment field-data filename. PhysiCell records it at + # microenvironment/domain/data/filename; this is distinct from the mesh file + # under mesh/voxels/filename. Reading it lets the reader locate the + # substrate .mat even when it does not follow the default naming convention. + microenvironment_file = None + me_elem = root.find('.//microenvironment') + if me_elem is not None: + data_fn = me_elem.find('.//data/filename') + if data_fn is not None and data_fn.text: + microenvironment_file = data_fn.text.strip() + # Extract cell type names and IDs cell_type_names = [] cell_type_ids = [] @@ -170,6 +181,7 @@ def parse_physicell_xml(xml_path: Path) -> PhysiCellMetadata: extra = { 'custom_labels': custom_labels, 'xml_path': str(xml_path), + 'microenvironment_file': microenvironment_file, } return PhysiCellMetadata( diff --git a/spatialtissuepy/synthetic/physicell/reader.py b/spatialtissuepy/synthetic/physicell/reader.py index 6ca8132..839b3b4 100644 --- a/spatialtissuepy/synthetic/physicell/reader.py +++ b/spatialtissuepy/synthetic/physicell/reader.py @@ -11,7 +11,7 @@ import xml.etree.ElementTree as ET from dataclasses import dataclass, field from pathlib import Path -from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union +from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple, Union import numpy as np import pandas as pd @@ -24,9 +24,54 @@ PhysiCellMetadata, get_cell_type_mapping, parse_cells_mat, + parse_microenvironment_mat, parse_physicell_xml, ) + +def _find_microenvironment_mat( + xml_path: Path, + microenvironment_file: Optional[str] = None, +) -> Optional[Path]: + """ + Locate the microenvironment ``.mat`` file for a frame. + + Resolution order: + + 1. The filename declared in the frame's XML (``microenvironment/domain/ + data/filename``), if given -- this survives renamed or relocated output. + 2. The default naming convention ``output{index:08d}_microenvironment0.mat`` + next to the XML. + + Parameters + ---------- + xml_path : Path + Path to the frame's ``output*.xml``. + microenvironment_file : str, optional + Filename from the XML, as parsed into + ``metadata.extra['microenvironment_file']``. + + Returns + ------- + Path or None + Path to the microenvironment ``.mat`` if found, else ``None``. + """ + xml_path = Path(xml_path) + + if microenvironment_file: + candidate = xml_path.parent / microenvironment_file + if candidate.exists(): + return candidate + + match = re.search(r'output(\d+)', xml_path.stem) + if match: + index = int(match.group(1)) + candidate = xml_path.parent / f'output{index:08d}_microenvironment0.mat' + if candidate.exists(): + return candidate + + return None + # ----------------------------------------------------------------------------- # PhysiCell TimeStep # ----------------------------------------------------------------------------- @@ -69,8 +114,13 @@ class PhysiCellTimeStep(ABMTimeStep): cells_mat_path: Path = None cell_type_mapping: Dict[int, str] = field(default_factory=dict) include_dead_cells: bool = False + microenvironment_mat_path: Optional[Path] = None _cell_data: Optional[Dict[str, np.ndarray]] = field(default=None, repr=False) _physicell_metadata: Optional[PhysiCellMetadata] = field(default=None, repr=False) + _microenvironment: Optional[Dict[str, Any]] = field(default=None, repr=False) + _me_loaded: bool = field(default=False, repr=False) + _voxel_tree: Any = field(default=None, repr=False) + _voxel_tree_ndim: Optional[int] = field(default=None, repr=False) def _load_cell_data(self) -> Dict[str, np.ndarray]: """Load cell data from MAT file (cached).""" @@ -99,6 +149,200 @@ def _load_metadata(self) -> PhysiCellMetadata: self._physicell_metadata = parse_physicell_xml(self.source_path) return self._physicell_metadata + def _load_microenvironment(self) -> Dict[str, Any]: + """ + Load the substrate field data for this frame (cached). + + Resolves the microenvironment ``.mat`` from the XML (or the default + naming convention) if a path was not supplied. A frame with no + microenvironment file yields empty substrates rather than raising, so + callers can probe ``substrates`` unconditionally. + """ + if self._me_loaded: + return self._microenvironment + + empty = { + 'voxel_positions': np.empty((0, 3)), + 'concentrations': {}, + 'raw_data': None, + } + + path = self.microenvironment_mat_path + if path is None: + try: + me_file = self._load_metadata().extra.get('microenvironment_file') + except (OSError, ET.ParseError): + me_file = None + path = _find_microenvironment_mat(self.source_path, me_file) + self.microenvironment_mat_path = path + + if path is None: + self._microenvironment = empty + else: + try: + substrate_names = self._load_metadata().substrate_names + except (OSError, ET.ParseError): + substrate_names = None + self._microenvironment = parse_microenvironment_mat( + path, substrate_names + ) + + self._me_loaded = True + return self._microenvironment + + @property + def substrate_names(self) -> List[str]: + """Names of the diffusible substrates, in matrix order.""" + return list(self._load_metadata().substrate_names) + + @property + def substrates(self) -> Dict[str, np.ndarray]: + """ + Substrate concentration fields, keyed by name. + + Each value is a ``(n_voxels,)`` array aligned with + :attr:`voxel_positions`. Empty if the frame has no microenvironment + file. + """ + return self._load_microenvironment()['concentrations'] + + @property + def voxel_positions(self) -> np.ndarray: + """``(n_voxels, 3)`` array of voxel-center coordinates.""" + return self._load_microenvironment()['voxel_positions'] + + def substrate_at( + self, + name: str, + x: Union[float, np.ndarray], + y: Union[float, np.ndarray], + z: Optional[Union[float, np.ndarray]] = None, + ) -> np.ndarray: + """ + Sample a substrate's environmental concentration at spatial points. + + Returns the concentration in the voxel nearest each query point, found + with a KD-tree over the voxel centers. This is the concentration of the + substrate *in the space around* a location; for the amount a cell has + actually taken up, use :meth:`internalized_substrates`. + + Parameters + ---------- + name : str + Substrate name; must be one of :attr:`substrate_names`. + x, y : float or array-like + Query coordinates. Scalars or equal-length arrays. + z : float or array-like, optional + Query z. If given, nearest-voxel search is 3-D; otherwise it uses + the x-y plane (appropriate for 2-D simulations). + + Returns + ------- + np.ndarray + Concentration at each query point, matching the input shape. + """ + me = self._load_microenvironment() + concentrations = me['concentrations'] + if name not in concentrations: + available = ', '.join(concentrations) or '(none)' + raise ValueError( + f"Unknown substrate {name!r}. Available: {available}" + ) + + voxels = me['voxel_positions'] + if voxels.shape[0] == 0: + raise ValueError( + f"No microenvironment voxels available for {self.source_path}" + ) + + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + if z is None: + query = np.column_stack([np.ravel(x), np.ravel(y)]) + ndim = 2 + else: + z = np.asarray(z, dtype=float) + query = np.column_stack([np.ravel(x), np.ravel(y), np.ravel(z)]) + ndim = 3 + + tree = self._voxel_kdtree(ndim) + _, idx = tree.query(query) + result = concentrations[name][idx] + + return result if x.ndim else result[0] + + def _voxel_kdtree(self, ndim: int): + """Build (and cache) a KD-tree over voxel centers for `ndim` axes.""" + if self._voxel_tree is None or self._voxel_tree_ndim != ndim: + from scipy.spatial import cKDTree + + voxels = self._load_microenvironment()['voxel_positions'] + self._voxel_tree = cKDTree(voxels[:, :ndim]) + self._voxel_tree_ndim = ndim + return self._voxel_tree + + def internalized_substrates( + self, + include_dead_cells: Optional[bool] = None, + ) -> pd.DataFrame: + """ + Per-cell internalized substrate amounts, one column per substrate. + + These come from PhysiCell's ``internalized_total_substrates`` field on + each cell, so they reflect the framework's actual uptake dynamics rather + than the environmental concentration at the cell's location. Use this + instead of :meth:`substrate_at` when you need what a cell has taken up, + not what surrounds it. + + Rows are in the same order as :attr:`positions`. + + .. note:: + PhysiCell only accumulates this field when the simulation enables + ``track_internalized_substrates_in_each_agent``. If that option was + off, the column is present but uniformly zero -- that is the model's + recorded value, not a parsing artifact. For the environmental + concentration, which is always available, use :meth:`substrate_at`. + + Parameters + ---------- + include_dead_cells : bool, optional + Whether to include dead cells. Defaults to the instance attribute. + + Returns + ------- + pd.DataFrame + Columns named by substrate; one row per cell. + + Raises + ------ + ValueError + If the frame's cell data does not record internalized substrates + (the model did not write the ``internalized_total_substrates`` + field, or the XML labels were unavailable). + """ + data = self._load_cell_data() + columns = data.get('columns') or {} + names = self.substrate_names + + field_keys = [f'internalized_total_substrates_{i}' for i in range(len(names))] + if not names or not all(k in columns for k in field_keys): + raise ValueError( + "This frame does not record internalized substrates. It " + "requires the PhysiCell 'internalized_total_substrates' field " + "and readable XML ; neither is optional here." + ) + + if include_dead_cells is None: + include_dead_cells = self.include_dead_cells + if include_dead_cells: + mask = np.ones(len(data['cell_types']), dtype=bool) + else: + mask = data['dead_flags'] == 0 + + return pd.DataFrame( + {name: columns[key][mask] for name, key in zip(names, field_keys)} + ) + @property def n_cells(self) -> int: """Number of cells at this time step.""" @@ -390,6 +634,18 @@ def n_timesteps(self) -> int: """Number of time steps.""" return len(self._timestep_files) + @property + def substrate_names(self) -> List[str]: + """ + Names of the diffusible substrates in this simulation. + + Read from the first frame's XML; empty if there are no time steps. + """ + if not self._timestep_files: + return [] + first_xml = self._timestep_files[0][1] + return list(parse_physicell_xml(first_xml).substrate_names) + @property def times(self) -> np.ndarray: """Array of simulation times.""" diff --git a/tests/test_physicell_microenvironment.py b/tests/test_physicell_microenvironment.py new file mode 100644 index 0000000..a60abc1 --- /dev/null +++ b/tests/test_physicell_microenvironment.py @@ -0,0 +1,237 @@ +""" +Tests for PhysiCell microenvironment / substrate access (v0.4.0). + +Covers two distinct notions of "how much substrate": + +- environmental concentration at a location, sampled from the voxel field via + ``PhysiCellTimeStep.substrate_at`` (nearest-voxel, KD-tree); +- the amount a cell has actually internalized, read per-cell from PhysiCell's + ``internalized_total_substrates`` field via + ``PhysiCellTimeStep.internalized_substrates``. + +The bundled example simulation is used for the equivalence checks; small +synthetic inputs cover ties, missing files, and error paths. +""" + +from pathlib import Path + +import numpy as np +import pytest +from scipy.io import loadmat, savemat + +from spatialtissuepy.synthetic.physicell import ( + PhysiCellTimeStep, + discover_physicell_timesteps, + read_physicell_simulation, + read_physicell_timestep, +) +from spatialtissuepy.synthetic.physicell.reader import _find_microenvironment_mat + + +@pytest.fixture +def example_physicell_dir(): + path = ( + Path(__file__).parent.parent + / 'examples' / 'sample_data' / 'example_physicell_sim' + ) + if not path.exists(): + pytest.skip(f"Example PhysiCell data not found at {path}") + return path + + +@pytest.fixture +def frame_xml(example_physicell_dir): + xmls = sorted(example_physicell_dir.glob('output*.xml')) + if not xmls: + pytest.skip("No frames in example simulation") + return xmls[len(xmls) // 2] + + +def _injected_timestep(voxels, concentrations): + """A PhysiCellTimeStep with a microenvironment injected, no files read.""" + ts = PhysiCellTimeStep(time=0.0, time_index=0, source_path=Path('none.xml')) + ts._microenvironment = { + 'voxel_positions': np.asarray(voxels, dtype=float), + 'concentrations': {k: np.asarray(v, dtype=float) + for k, v in concentrations.items()}, + 'raw_data': None, + } + ts._me_loaded = True + return ts + + +class TestSubstrateFields: + def test_substrate_names(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert 'oxygen' in ts.substrate_names + + def test_substrates_shape(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + oxygen = ts.substrates['oxygen'] + assert oxygen.ndim == 1 + assert oxygen.shape[0] == ts.voxel_positions.shape[0] + + def test_voxel_positions_are_3d(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert ts.voxel_positions.shape[1] == 3 + + +class TestSubstrateAt: + def test_matches_bruteforce_argmin(self, frame_xml): + """The report's acceptance criterion: exact match to nearest-voxel.""" + ts = read_physicell_timestep(frame_xml) + coords = ts.to_spatial_data().coordinates + xs, ys = coords[:, 0], coords[:, 1] + + got = ts.substrate_at('oxygen', xs, ys) + + vox = ts.voxel_positions[:, :2] + field = ts.substrates['oxygen'] + brute = np.array([ + field[np.argmin((vox[:, 0] - x) ** 2 + (vox[:, 1] - y) ** 2)] + for x, y in zip(xs, ys) + ]) + np.testing.assert_array_equal(got, brute) + + def test_scalar_query_returns_scalar(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + val = ts.substrate_at('oxygen', 0.0, 0.0) + assert np.isscalar(val) or val.ndim == 0 + + def test_array_query_preserves_length(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + out = ts.substrate_at('oxygen', np.array([0.0, 10.0, 20.0]), + np.array([0.0, 10.0, 20.0])) + assert out.shape == (3,) + + def test_unknown_substrate_raises(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + with pytest.raises(ValueError, match="Unknown substrate"): + ts.substrate_at('unobtanium', 0.0, 0.0) + + def test_tie_is_deterministic_and_valid(self): + """A point equidistant between two voxels returns one of their values. + + KD-tree tie-breaking need not match np.argmin, so this asserts the + result is one of the tied voxels and is stable, rather than pinning a + specific choice. + """ + voxels = [[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]] + ts = _injected_timestep(voxels, {'oxygen': [1.0, 2.0]}) + + first = ts.substrate_at('oxygen', 5.0, 0.0) + second = ts.substrate_at('oxygen', 5.0, 0.0) + + assert first in (1.0, 2.0) + assert first == second # deterministic + + def test_3d_query(self): + voxels = [[0, 0, 0], [0, 0, 10]] + ts = _injected_timestep(voxels, {'o2': [5.0, 9.0]}) + assert ts.substrate_at('o2', 0.0, 0.0, 1.0) == 5.0 + assert ts.substrate_at('o2', 0.0, 0.0, 9.0) == 9.0 + + +class TestInternalizedSubstrates: + def test_columns_match_raw_matrix(self, frame_xml, example_physicell_dir): + """internalized_total_substrates maps to rows 92+ in matrix order.""" + ts = read_physicell_timestep(frame_xml, include_dead_cells=True) + matrix = loadmat(str(ts.cells_mat_path))['cells'] + if matrix.shape[0] < 97: + pytest.skip("Frame lacks internalized_total_substrates rows") + + inte = ts.internalized_substrates(include_dead_cells=True) + for i, name in enumerate(ts.substrate_names): + np.testing.assert_array_equal(inte[name].to_numpy(), matrix[92 + i, :]) + + def test_columns_named_by_substrate(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + inte = ts.internalized_substrates() + assert list(inte.columns) == ts.substrate_names + + def test_row_count_follows_dead_filter(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + assert len(ts.internalized_substrates()) == ts.n_cells + assert len(ts.internalized_substrates(include_dead_cells=True)) == \ + ts.n_cells_total + + def test_row_order_matches_positions(self, frame_xml): + ts = read_physicell_timestep(frame_xml) + inte = ts.internalized_substrates() + assert len(inte) == ts.positions.shape[0] + + def test_raises_without_field(self, tmp_path): + """A frame whose cell data lacks the field raises informatively.""" + # Wider than tall so orientation is unambiguous; no , so the + # parser exposes no 'columns' and no substrate names are declared. + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((5, 20))}) + (tmp_path / 'output00000000.xml').write_text( + "" + ) + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + with pytest.raises(ValueError, match="internalized"): + ts.internalized_substrates() + + +class TestFileResolution: + def test_resolves_from_xml_filename(self, frame_xml, example_physicell_dir): + expected = example_physicell_dir / f'{frame_xml.stem}_microenvironment0.mat' + got = _find_microenvironment_mat( + frame_xml, f'{frame_xml.stem}_microenvironment0.mat' + ) + assert got == expected + + def test_falls_back_to_naming_convention(self, frame_xml): + # No filename hint -> should still find it by convention. + got = _find_microenvironment_mat(frame_xml, None) + assert got is not None and got.exists() + + def test_returns_none_when_absent(self, tmp_path): + (tmp_path / 'output00000000.xml').write_text("") + assert _find_microenvironment_mat( + tmp_path / 'output00000000.xml', None + ) is None + + +class TestEmptyMicroenvironment: + def test_missing_file_yields_empty(self, tmp_path): + """No microenvironment file -> empty substrates, no exception.""" + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((30, 4))}) + (tmp_path / 'output00000000.xml').write_text( + "" + ) + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + + assert ts.substrates == {} + assert ts.voxel_positions.shape == (0, 3) + + def test_substrate_at_on_empty_raises(self, tmp_path): + savemat(str(tmp_path / 'output00000000_cells.mat'), + {'cells': np.zeros((30, 4))}) + (tmp_path / 'output00000000.xml').write_text("") + ts = read_physicell_timestep(tmp_path / 'output00000000.xml') + with pytest.raises(ValueError): + ts.substrate_at('oxygen', 0.0, 0.0) + + +class TestSimulationIntegration: + def test_substrate_names(self, example_physicell_dir): + sim = read_physicell_simulation(example_physicell_dir) + assert 'oxygen' in sim.substrate_names + + def test_timesteps_have_substrate_access(self, example_physicell_dir): + sim = read_physicell_simulation(example_physicell_dir) + ts = sim.get_timestep(0) + assert ts.voxel_positions.shape[0] > 0 + + +class TestBackwardCompatibility: + def test_discover_still_returns_3_tuples(self, example_physicell_dir): + """Adding microenvironment support must not change discover's shape.""" + steps = discover_physicell_timesteps(example_physicell_dir) + assert steps and all(len(s) == 3 for s in steps) + index, xml_path, mat_path = steps[0] # unpacks as before + assert xml_path.suffix == '.xml' + assert mat_path.suffix == '.mat'