diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8af0a81..2872ce6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,78 @@ 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]
+
+PhysiCell reader correctness fixes, reported against v0.2.0 by the
+llm-abm-consistency harness (OHSU ChangLab) and reproduced against both their
+PhysiCell 1.14.2 output and this repository's bundled example simulation.
+
+### Fixed
+
+- **Cell matrix orientation is now resolved from the frame's declared variable
+  count** rather than by comparing the matrix dimensions. The old heuristic
+  (`if shape[0] > shape[1]: transpose`) assumed every frame holds more cells
+  than the model has variables, and silently returned transposed garbage for
+  any frame that does not — sparse initial conditions, early time steps, small
+  explants, and simulations approaching extinction. A frame that matches
+  neither orientation now raises `ValueError` instead of returning data. When
+  no labels are available, the old heuristic still applies but emits a
+  `UserWarning` naming the ambiguity.
+- **`PhysiCellTimeStep.to_dataframe()` now honors `include_dead_cells`.** It
+  previously ignored the flag while `n_cells`, `positions`, `to_spatial_data()`
+  and `cell_counts_by_type()` all respected it, so the same object reported
+  different populations depending on the accessor and row indices did not align
+  between `to_dataframe()` and `positions`.
+- **`parse_microenvironment_mat` validates matrix orientation** against the
+  declared substrate count when `substrate_names` is supplied, converting a
+  would-be silent misread into a clear error.
+- **`__version__` is now read from the installed distribution metadata**, so it
+  cannot drift from `pyproject.toml`. It reported `0.1.0` at the `v0.2.0` tag,
+  which misattributed results in downstream provenance capture. A regression
+  test asserts the two agree.
+
+### Added
+
+- **`parse_cells_mat(..., labels=...)`**: derive column positions from the
+  MultiCellDS `<labels>` block that `parse_physicell_xml` already parses. The
+  parsed labels were previously stored at `metadata.extra['custom_labels']` and
+  never read anywhere in the package. Precedence is `index_mapping` >
+  `labels` > the existing row-count autodetect, so behavior is unchanged when
+  labels are not supplied.
+- **`'columns'` key in the `parse_cells_mat` result**: every labelled column,
+  including standard PhysiCell fields absent from the hard-coded index tables
+  (`is_motile`, `migration_speed`) and model-specific custom variables. The
+  return dict was previously a closed set, so anything else the model wrote was
+  unreachable.
+- **`'orientation'` key** reporting the on-disk layout of `raw_data`, either
+  `'variables_x_cells'` or `'cells_x_variables'`.
+- **`expand_cell_labels()` and `declared_variable_count()`** helpers for
+  working with `<labels>` blocks. Vector labels expand as `{name}_x/_y/_z` for
+  `size == 3` and `{name}_0 … {name}_{n-1}` otherwise.
+- **`to_dataframe(extra_columns=True)`** appends every labelled column. Off by
+  default, since a full PhysiCell frame carries 150+ variables.
+- **`to_dataframe(include_dead_cells=...)`** to override the instance attribute
+  per call.
+- Regression tests (`tests/test_physicell_labels.py`, 25 tests) covering label
+  expansion, all four orientation cases, cross-model label isolation, the
+  microenvironment guard, dead-cell filtering, and version consistency.
+
+### Changed
+
+- **`raw_data` from `parse_cells_mat` is now returned exactly as loaded from
+  disk.** It was previously the reoriented matrix, so callers could not use it
+  to recover the truth when the orientation heuristic misfired. Callers relying
+  on the transposed `raw_data` should consult the new `'orientation'` key.
+- `PhysiCellTimeStep._load_cell_data()` and `to_trajectory_dataframe()` now
+  resolve columns from each frame's own XML, so models with differing variable
+  counts parse correctly in the same session.
+
+### Known gaps
+
+- `parse_microenvironment_mat` is implemented and exported but still not wired
+  into `PhysiCellTimeStep` / `PhysiCellSimulation`; substrate fields remain
+  reachable only through the parser. Planned for v0.3.0.
+
 ## [0.2.0] - 2026-06-26
 
 First beta release. This release adds AI agent access via an MCP server, a
diff --git a/spatialtissuepy/__init__.py b/spatialtissuepy/__init__.py
index 3e090af..75c116f 100644
--- a/spatialtissuepy/__init__.py
+++ b/spatialtissuepy/__init__.py
@@ -21,7 +21,16 @@
 io : Input/output utilities
 """
 
-__version__ = "0.1.0"
+from importlib.metadata import PackageNotFoundError
+from importlib.metadata import version as _pkg_version
+
+try:
+    # Single source of truth: the version declared in pyproject.toml, as
+    # recorded in the installed distribution metadata.
+    __version__ = _pkg_version("spatialtissuepy")
+except PackageNotFoundError:  # pragma: no cover - running from a source tree
+    __version__ = "0.0.0.dev0"
+
 __author__ = "spatialtissuepy developers"
 
 from spatialtissuepy.core.spatial_data import SpatialTissueData
diff --git a/spatialtissuepy/synthetic/physicell/parser.py b/spatialtissuepy/synthetic/physicell/parser.py
index 874ba67..59fc12f 100644
--- a/spatialtissuepy/synthetic/physicell/parser.py
+++ b/spatialtissuepy/synthetic/physicell/parser.py
@@ -7,6 +7,7 @@
 
 from __future__ import annotations
 
+import warnings
 import xml.etree.ElementTree as ET
 from dataclasses import dataclass
 from pathlib import Path
@@ -287,10 +288,119 @@ def get_cell_type_mapping(
 }
 
 
+def expand_cell_labels(labels: Dict[int, Tuple[str, int]]) -> Dict[int, str]:
+    """
+    Expand a MultiCellDS ``<labels>`` block into a flat column mapping.
+
+    Vector-valued labels occupy consecutive columns in the cell matrix. A label
+    of ``size == 3`` is expanded to ``{name}_x/_y/_z`` (the convention PhysiCell
+    uses for spatial vectors); any other ``size > 1`` is expanded positionally
+    to ``{name}_0 ... {name}_{size-1}``.
+
+    Parameters
+    ----------
+    labels : dict
+        Mapping of ``{start_index: (name, size)}``, as produced by
+        :func:`parse_physicell_xml` at ``metadata.extra['custom_labels']``.
+
+    Returns
+    -------
+    dict
+        Mapping of ``{column_index: column_name}``.
+
+    Examples
+    --------
+    >>> expand_cell_labels({0: ('ID', 1), 1: ('position', 3)})
+    {0: 'ID', 1: 'position_x', 2: 'position_y', 3: 'position_z'}
+    """
+    expanded: Dict[int, str] = {}
+
+    for start_index, (name, size) in labels.items():
+        if size == 1:
+            expanded[start_index] = name
+        elif size == 3:
+            for offset, axis in enumerate('xyz'):
+                expanded[start_index + offset] = f'{name}_{axis}'
+        else:
+            for offset in range(size):
+                expanded[start_index + offset] = f'{name}_{offset}'
+
+    return expanded
+
+
+def declared_variable_count(labels: Dict[int, Tuple[str, int]]) -> int:
+    """
+    Number of matrix rows implied by a ``<labels>`` block.
+
+    This accounts for the width of the final label, so a trailing vector label
+    is not undercounted the way ``max(labels) + 1`` would be.
+
+    Parameters
+    ----------
+    labels : dict
+        Mapping of ``{start_index: (name, size)}``.
+
+    Returns
+    -------
+    int
+        Expected number of variables (rows) per cell.
+    """
+    return max(index + size for index, (_, size) in labels.items())
+
+
+def _orient_cell_matrix(
+    cell_matrix: np.ndarray,
+    n_expected: Optional[int],
+    mat_path: Path,
+) -> Tuple[np.ndarray, str]:
+    """
+    Normalize a cell matrix to ``(n_variables, n_cells)``.
+
+    When the declared variable count is known, orientation is resolved by
+    matching it against the array shape -- unambiguous even when a frame holds
+    fewer cells than the model has variables. Without it, fall back to the
+    legacy magnitude heuristic and warn, since that assumption silently
+    corrupts small frames.
+
+    Returns
+    -------
+    (np.ndarray, str)
+        The matrix in ``(n_variables, n_cells)`` orientation, and the
+        orientation of the input as loaded from disk.
+    """
+    n_rows, n_cols = cell_matrix.shape
+
+    if n_expected is not None:
+        # Prefer axis 0 as variables, so a square matrix resolves the way
+        # MultiCellDS actually writes it.
+        if n_rows == n_expected:
+            return cell_matrix, 'variables_x_cells'
+        if n_cols == n_expected:
+            return cell_matrix.T, 'cells_x_variables'
+        raise ValueError(
+            f"cell matrix {cell_matrix.shape} matches neither orientation for "
+            f"{n_expected} declared variables in {mat_path}"
+        )
+
+    if n_rows > n_cols:
+        warnings.warn(
+            f"Guessing cell matrix orientation for {mat_path} from its shape "
+            f"{cell_matrix.shape} because no <labels> block was supplied. This "
+            f"is wrong for any frame holding fewer cells than the model has "
+            f"variables. Pass labels= to resolve orientation unambiguously.",
+            UserWarning,
+            stacklevel=3,
+        )
+        return cell_matrix.T, 'cells_x_variables'
+
+    return cell_matrix, 'variables_x_cells'
+
+
 def parse_cells_mat(
     mat_path: Path,
     cell_type_mapping: Optional[Dict[int, str]] = None,
-    index_mapping: Optional[Dict[str, int]] = None
+    index_mapping: Optional[Dict[str, int]] = None,
+    labels: Optional[Dict[int, Tuple[str, int]]] = None,
 ) -> Dict[str, np.ndarray]:
     """
     Parse a PhysiCell cells_physicell.mat file.
@@ -303,6 +413,13 @@ def parse_cells_mat(
         Mapping from cell type IDs to names.
     index_mapping : dict, optional
         Custom row index mapping for different PhysiCell versions.
+    labels : dict, optional
+        The ``<labels>`` block from the frame's XML, as
+        ``{start_index: (name, size)}``. When supplied, column positions and
+        matrix orientation are resolved from the frame's own declaration rather
+        than from a hard-coded table, and every labelled column is returned
+        under ``'columns'``. Precedence for column positions is
+        ``index_mapping`` > ``labels`` > row-count autodetect.
 
     Returns
     -------
@@ -314,8 +431,18 @@ def parse_cells_mat(
         - 'volumes': (n_cells,) array of total volumes
         - 'radii': (n_cells,) array of cell radii
         - 'phases': (n_cells,) array of cell cycle phases
+        - 'dead_flags': (n_cells,) array of dead flags
         - 'ids': (n_cells,) array of cell IDs
-        - 'raw_data': Full matrix from file
+        - 'columns': dict of every labelled column, keyed by name. Empty when
+          no ``labels`` were supplied.
+        - 'raw_data': Full matrix exactly as loaded from disk
+        - 'orientation': layout of ``raw_data``, either ``'variables_x_cells'``
+          or ``'cells_x_variables'``
+
+    Notes
+    -----
+    ``raw_data`` is returned unmodified, so it round-trips to what
+    ``scipy.io.loadmat`` returned; use ``orientation`` to interpret it.
     """
     from scipy.io import loadmat
 
@@ -340,13 +467,10 @@ def parse_cells_mat(
     if cell_matrix is None:
         raise ValueError(f"Could not find cell data in {mat_path}")
 
-    # Ensure cells are columns
-    if cell_matrix.shape[0] > cell_matrix.shape[1]:
-        cell_matrix = cell_matrix.T
+    raw_data = cell_matrix
 
-    n_cells = cell_matrix.shape[1]
-
-    if n_cells == 0:
+    # An empty frame carries no orientation signal; report it as written.
+    if 0 in cell_matrix.shape:
         return {
             'positions': np.empty((0, 3)),
             'cell_types': np.array([], dtype=str),
@@ -354,15 +478,31 @@ def parse_cells_mat(
             'volumes': np.array([]),
             'radii': np.array([]),
             'phases': np.array([], dtype=int),
+            'dead_flags': np.array([], dtype=int),
             'ids': np.array([], dtype=int),
-            'raw_data': cell_matrix,
+            'columns': {},
+            'raw_data': raw_data,
+            'orientation': 'variables_x_cells',
         }
 
-    # Use provided or default index mapping
+    label_columns = expand_cell_labels(labels) if labels else {}
+    n_expected = declared_variable_count(labels) if labels else None
+
+    cell_matrix, orientation = _orient_cell_matrix(
+        cell_matrix, n_expected, mat_path
+    )
+
+    n_cells = cell_matrix.shape[1]
+
+    # Column positions: explicit index_mapping wins, then the frame's own
+    # labels, then the legacy row-count autodetect.
     if index_mapping is None:
-        # Auto-detect based on matrix shape
-        # PhysiCell 1.10+ typically has 150+ rows per cell
-        if cell_matrix.shape[0] >= 30:
+        if label_columns:
+            index_mapping = {
+                name: index for index, name in label_columns.items()
+            }
+        elif cell_matrix.shape[0] >= 30:
+            # PhysiCell 1.10+ typically has 150+ rows per cell
             index_mapping = CELL_DATA_INDICES_V2
         else:
             index_mapping = CELL_DATA_INDICES_LEGACY
@@ -423,6 +563,14 @@ def parse_cells_mat(
         # Fall back to inferring from phase code
         dead_flags = np.array([1 if p >= 100 else 0 for p in phases])
 
+    # Every labelled column, including ones with no entry in the index tables
+    # (is_motile, migration_speed) and model-specific custom variables.
+    columns = {
+        name: cell_matrix[index, :]
+        for index, name in sorted(label_columns.items())
+        if index < cell_matrix.shape[0]
+    }
+
     return {
         'positions': positions,
         'cell_types': cell_types,
@@ -432,7 +580,9 @@ def parse_cells_mat(
         'phases': phases,
         'dead_flags': dead_flags,
         'ids': ids,
-        'raw_data': cell_matrix,
+        'columns': columns,
+        'raw_data': raw_data,
+        'orientation': orientation,
     }
 
 
@@ -479,7 +629,23 @@ def parse_microenvironment_mat(
     if me_matrix is None:
         raise ValueError(f"Could not find microenvironment data in {mat_path}")
 
-    # Structure: rows 0-2 are x,y,z; row 3 is volume; rows 4+ are substrates
+    raw_data = me_matrix
+
+    # Structure: rows 0-2 are x,y,z; row 3 is volume; rows 4+ are substrates.
+    # PhysiCell always writes (4 + n_substrates, n_voxels), but validate it
+    # rather than assume -- an unchecked transpose here is silent corruption.
+    if substrate_names is not None:
+        n_expected = 4 + len(substrate_names)
+        if me_matrix.shape[0] != n_expected:
+            if me_matrix.shape[1] == n_expected:
+                me_matrix = me_matrix.T
+            else:
+                raise ValueError(
+                    f"microenvironment matrix {me_matrix.shape} matches neither "
+                    f"orientation for {len(substrate_names)} declared substrates "
+                    f"in {mat_path}"
+                )
+
     voxel_positions = me_matrix[:3, :].T  # (n_voxels, 3)
 
     # Extract substrate concentrations
@@ -495,7 +661,7 @@ def parse_microenvironment_mat(
     return {
         'voxel_positions': voxel_positions,
         'concentrations': concentrations,
-        'raw_data': me_matrix,
+        'raw_data': raw_data,
     }
 
 
diff --git a/spatialtissuepy/synthetic/physicell/reader.py b/spatialtissuepy/synthetic/physicell/reader.py
index 86c52c3..6ca8132 100644
--- a/spatialtissuepy/synthetic/physicell/reader.py
+++ b/spatialtissuepy/synthetic/physicell/reader.py
@@ -8,6 +8,7 @@
 from __future__ import annotations
 
 import re
+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
@@ -74,9 +75,21 @@ class PhysiCellTimeStep(ABMTimeStep):
     def _load_cell_data(self) -> Dict[str, np.ndarray]:
         """Load cell data from MAT file (cached)."""
         if self._cell_data is None:
+            # Resolve columns and matrix orientation from this frame's own
+            # <labels> block, so models with differing variable counts each
+            # parse against their own declaration.
+            labels = None
+            try:
+                labels = self._load_metadata().extra.get('custom_labels') or None
+            except (OSError, ET.ParseError):
+                # No readable XML: parse_cells_mat falls back to its row-count
+                # heuristic and warns.
+                pass
+
             self._cell_data = parse_cells_mat(
                 self.cells_mat_path,
-                self.cell_type_mapping
+                self.cell_type_mapping,
+                labels=labels,
             )
         return self._cell_data
 
@@ -180,35 +193,72 @@ def to_spatial_data(self) -> SpatialTissueData:
             metadata=extra_metadata,
         )
 
-    def to_dataframe(self) -> pd.DataFrame:
+    def to_dataframe(
+        self,
+        include_dead_cells: Optional[bool] = None,
+        extra_columns: bool = False,
+    ) -> pd.DataFrame:
         """
         Convert cell data to a pandas DataFrame.
 
+        Parameters
+        ----------
+        include_dead_cells : bool, optional
+            Whether to include dead cells. Defaults to the instance attribute,
+            so rows stay aligned with ``positions`` and ``to_spatial_data()``.
+            Pass ``True`` to recover the unfiltered frame.
+        extra_columns : bool, default False
+            Append every labelled column from the frame's XML, including
+            standard fields absent from the index tables (``is_motile``,
+            ``migration_speed``) and model-specific custom variables. Off by
+            default because a full PhysiCell frame carries 150+ variables.
+
         Returns
         -------
         pd.DataFrame
-            DataFrame with cell properties.
+            DataFrame with cell properties, in the same row order as
+            ``positions``.
         """
         data = self._load_cell_data()
 
+        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
+
         df = pd.DataFrame({
-            'cell_id': data['ids'],
-            'x': data['positions'][:, 0],
-            'y': data['positions'][:, 1],
-            'z': data['positions'][:, 2],
-            'cell_type': data['cell_types'],
-            'cell_type_id': data['cell_type_ids'],
-            'volume': data['volumes'],
-            'radius': data['radii'],
-            'phase': data['phases'],
-            'is_dead': data['dead_flags'].astype(bool),
-            'is_alive': ~data['dead_flags'].astype(bool),
+            'cell_id': data['ids'][mask],
+            'x': data['positions'][mask, 0],
+            'y': data['positions'][mask, 1],
+            'z': data['positions'][mask, 2],
+            'cell_type': data['cell_types'][mask],
+            'cell_type_id': data['cell_type_ids'][mask],
+            'volume': data['volumes'][mask],
+            'radius': data['radii'][mask],
+            'phase': data['phases'][mask],
+            'is_dead': data['dead_flags'][mask].astype(bool),
+            'is_alive': ~data['dead_flags'][mask].astype(bool),
         })
 
-        df['time'] = self.time
-        df['time_index'] = self.time_index
+        trailing = pd.DataFrame({
+            'time': np.full(len(df), self.time),
+            'time_index': np.full(len(df), self.time_index),
+        })
+
+        if extra_columns:
+            # Build in one concat; a full PhysiCell frame adds 150+ columns and
+            # inserting them one at a time fragments the frame badly.
+            extra = pd.DataFrame({
+                name: values[mask]
+                for name, values in data['columns'].items()
+                if name not in df.columns
+            })
+            return pd.concat([df, extra, trailing], axis=1)
 
-        return df
+        return pd.concat([df, trailing], axis=1)
 
     def cell_counts_by_type(self) -> Dict[str, int]:
         """Get cell counts by type."""
diff --git a/tests/test_physicell_labels.py b/tests/test_physicell_labels.py
new file mode 100644
index 0000000..a47a58a
--- /dev/null
+++ b/tests/test_physicell_labels.py
@@ -0,0 +1,330 @@
+"""
+Regression tests for PhysiCell reader defects fixed in v0.2.1.
+
+These cover the label-driven column extraction and orientation handling
+reported by the llm-abm-consistency harness (OHSU ChangLab, 2026-07-21).
+
+The bundled example simulation cannot exercise the orientation bug on its own:
+it writes 154 variables per frame and its smallest frame holds 906 cells, so
+``n_cells > n_variables`` holds everywhere. These tests build small matrices
+where that assumption fails.
+"""
+
+import importlib.metadata
+import warnings
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import pytest
+from scipy.io import loadmat, savemat
+
+import spatialtissuepy
+from spatialtissuepy.synthetic.physicell import (
+    discover_physicell_timesteps,
+    read_physicell_timestep,
+)
+from spatialtissuepy.synthetic.physicell.parser import (
+    declared_variable_count,
+    expand_cell_labels,
+    parse_cells_mat,
+    parse_microenvironment_mat,
+    parse_physicell_xml,
+)
+
+N_VARS = 87
+
+
+@pytest.fixture
+def example_physicell_dir():
+    """Path to example PhysiCell simulation output folder."""
+    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 labels():
+    """A minimal <labels> block in the shape parse_physicell_xml produces."""
+    return {
+        0: ('ID', 1),
+        1: ('position', 3),
+        4: ('total_volume', 1),
+        5: ('cell_type', 1),
+        7: ('current_phase', 1),
+        26: ('dead', 1),
+        37: ('radius', 1),
+        50: ('is_motile', 1),
+        52: ('migration_speed', 1),
+        86: ('phenotype_state', 1),
+    }
+
+
+@pytest.fixture
+def cell_matrix():
+    """A (87, 40) matrix with each cell's values keyed off its column index."""
+    matrix = np.zeros((N_VARS, 40), dtype=float)
+    for cell in range(40):
+        matrix[:, cell] = np.arange(N_VARS) * 1000.0 + cell
+    return matrix
+
+
+def _write(tmp_path, name, matrix):
+    path = tmp_path / name
+    savemat(str(path), {'cells': matrix})
+    return path
+
+
+class TestLabelExpansion:
+    """Expansion of the MultiCellDS <labels> block into flat columns."""
+
+    def test_scalar_labels_keep_their_name(self):
+        assert expand_cell_labels({0: ('ID', 1)}) == {0: 'ID'}
+
+    def test_size_three_expands_to_xyz(self):
+        assert expand_cell_labels({1: ('position', 3)}) == {
+            1: 'position_x', 2: 'position_y', 3: 'position_z'
+        }
+
+    def test_other_vector_sizes_expand_positionally(self):
+        assert expand_cell_labels({5: ('state', 2)}) == {
+            5: 'state_0', 6: 'state_1'
+        }
+
+    def test_variable_count_accounts_for_trailing_vector_width(self):
+        # max(labels) + 1 would undercount this by two.
+        assert declared_variable_count({0: ('ID', 1), 1: ('position', 3)}) == 4
+
+
+class TestOrientation:
+    """Issue #2: orientation must come from the label count, not magnitude."""
+
+    def test_fewer_cells_than_variables(self, tmp_path, labels, cell_matrix):
+        """An 87x40 matrix with 87 declared labels parses as 40 cells."""
+        path = _write(tmp_path, 'small.mat', cell_matrix)
+        result = parse_cells_mat(path, labels=labels)
+
+        assert result['positions'].shape[0] == 40
+        np.testing.assert_array_equal(
+            result['positions'][:, 0], cell_matrix[1, :]
+        )
+
+    def test_genuine_on_disk_transpose(self, tmp_path, labels, cell_matrix):
+        """A 40x87 matrix also parses as 40 cells."""
+        path = _write(tmp_path, 'transposed.mat', cell_matrix.T)
+        result = parse_cells_mat(path, labels=labels)
+
+        assert result['positions'].shape[0] == 40
+        np.testing.assert_array_equal(
+            result['positions'][:, 0], cell_matrix[1, :]
+        )
+        assert result['orientation'] == 'cells_x_variables'
+
+    def test_square_matrix_prefers_variables_on_axis_zero(self, tmp_path, labels):
+        square = np.zeros((N_VARS, N_VARS), dtype=float)
+        for cell in range(N_VARS):
+            square[:, cell] = np.arange(N_VARS) * 1000.0 + cell
+
+        result = parse_cells_mat(_write(tmp_path, 'sq.mat', square), labels=labels)
+
+        assert result['positions'].shape[0] == N_VARS
+        np.testing.assert_array_equal(result['positions'][:, 0], square[1, :])
+        assert result['orientation'] == 'variables_x_cells'
+
+    def test_neither_orientation_raises(self, tmp_path, labels):
+        path = _write(tmp_path, 'bad.mat', np.zeros((99, 33)))
+
+        with pytest.raises(ValueError, match='matches neither orientation'):
+            parse_cells_mat(path, labels=labels)
+
+    def test_raw_data_round_trips(self, tmp_path, labels, cell_matrix):
+        """raw_data is what loadmat returned, not the reoriented view."""
+        path = _write(tmp_path, 'raw.mat', cell_matrix)
+        result = parse_cells_mat(path, labels=labels)
+
+        np.testing.assert_array_equal(result['raw_data'], loadmat(str(path))['cells'])
+        assert result['orientation'] == 'variables_x_cells'
+
+    def test_warns_when_orientation_must_be_guessed(self, tmp_path, cell_matrix):
+        path = _write(tmp_path, 'noguess.mat', cell_matrix)
+
+        with pytest.warns(UserWarning, match='Guessing cell matrix orientation'):
+            parse_cells_mat(path)
+
+
+class TestLabelledColumns:
+    """Issue #1: every labelled column must be reachable."""
+
+    def test_non_indexed_standard_fields(self, tmp_path, labels, cell_matrix):
+        result = parse_cells_mat(
+            _write(tmp_path, 'cols.mat', cell_matrix), labels=labels
+        )
+
+        np.testing.assert_array_equal(
+            result['columns']['is_motile'], cell_matrix[50, :]
+        )
+        np.testing.assert_array_equal(
+            result['columns']['migration_speed'], cell_matrix[52, :]
+        )
+
+    def test_model_specific_custom_column(self, tmp_path, labels, cell_matrix):
+        result = parse_cells_mat(
+            _write(tmp_path, 'custom.mat', cell_matrix), labels=labels
+        )
+
+        np.testing.assert_array_equal(
+            result['columns']['phenotype_state'], cell_matrix[86, :]
+        )
+
+    def test_differing_label_counts_do_not_leak(self, tmp_path, labels, cell_matrix):
+        """Two models parsed in one session each resolve from their own XML."""
+        model_b_labels = dict(labels)
+        model_b_labels[87] = ('hif', 1)
+        model_b_labels[88] = ('hypoxia_timer', 1)
+        model_b_matrix = np.vstack([
+            cell_matrix,
+            np.full((2, cell_matrix.shape[1]), 7.0),
+        ])
+
+        result_a = parse_cells_mat(
+            _write(tmp_path, 'a.mat', cell_matrix), labels=labels
+        )
+        result_b = parse_cells_mat(
+            _write(tmp_path, 'b.mat', model_b_matrix), labels=model_b_labels
+        )
+
+        assert 'hif' not in result_a['columns']
+        assert 'hif' in result_b['columns']
+        # Index 86 carries a different meaning in each model.
+        assert result_a['positions'].shape[0] == result_b['positions'].shape[0]
+
+    def test_no_labels_preserves_legacy_output(self, tmp_path):
+        """With labels=None, output matches v0.2.0 for a default-layout frame."""
+        wide = np.zeros((N_VARS, 200), dtype=float)
+        for cell in range(200):
+            wide[:, cell] = np.arange(N_VARS) * 1000.0 + cell
+        path = _write(tmp_path, 'wide.mat', wide)
+
+        result = parse_cells_mat(path)
+
+        assert result['columns'] == {}
+        np.testing.assert_array_equal(result['positions'][:, 0], wide[1, :])
+        np.testing.assert_array_equal(result['dead_flags'], wide[26, :].astype(int))
+
+
+class TestMicroenvironmentGuard:
+    """Issue #6: validate the microenvironment matrix orientation."""
+
+    def test_correct_orientation_parses(self, tmp_path):
+        me = np.arange(6 * 50, dtype=float).reshape(6, 50)
+        path = tmp_path / 'me.mat'
+        savemat(str(path), {'multiscale_microenvironment': me})
+
+        result = parse_microenvironment_mat(path, ['oxygen', 'glucose'])
+
+        assert result['voxel_positions'].shape == (50, 3)
+        np.testing.assert_array_equal(result['concentrations']['oxygen'], me[4, :])
+
+    def test_transposed_input_is_corrected(self, tmp_path):
+        me = np.arange(6 * 50, dtype=float).reshape(6, 50)
+        path = tmp_path / 'me_t.mat'
+        savemat(str(path), {'multiscale_microenvironment': me.T})
+
+        result = parse_microenvironment_mat(path, ['oxygen', 'glucose'])
+
+        np.testing.assert_array_equal(result['concentrations']['oxygen'], me[4, :])
+
+    def test_mismatched_substrate_count_raises(self, tmp_path):
+        path = tmp_path / 'me_bad.mat'
+        savemat(str(path), {'multiscale_microenvironment': np.zeros((9, 11))})
+
+        with pytest.raises(ValueError, match='matches neither orientation'):
+            parse_microenvironment_mat(path, ['oxygen', 'glucose'])
+
+
+class TestToDataFrameDeadCells:
+    """Issue #4: to_dataframe must honour include_dead_cells."""
+
+    @pytest.fixture
+    def timestep_with_dead(self, example_physicell_dir):
+        for _, xml_path, _ in discover_physicell_timesteps(example_physicell_dir):
+            timestep = read_physicell_timestep(xml_path)
+            if timestep.n_dead_cells > 0:
+                return timestep
+        pytest.skip("No frame with dead cells in the example simulation")
+
+    def test_excludes_dead_by_default(self, timestep_with_dead):
+        df = timestep_with_dead.to_dataframe()
+
+        assert len(df) == timestep_with_dead.n_cells
+        assert not df['is_dead'].any()
+
+    def test_includes_dead_when_requested(self, timestep_with_dead):
+        df = timestep_with_dead.to_dataframe(include_dead_cells=True)
+
+        assert len(df) == timestep_with_dead.n_cells_total
+
+    def test_row_order_matches_positions(self, timestep_with_dead):
+        df = timestep_with_dead.to_dataframe()
+
+        np.testing.assert_allclose(
+            df[['x', 'y', 'z']].to_numpy(), timestep_with_dead.positions
+        )
+
+    def test_extra_columns_are_opt_in(self, timestep_with_dead):
+        base = timestep_with_dead.to_dataframe()
+        extended = timestep_with_dead.to_dataframe(extra_columns=True)
+
+        assert extended.shape[1] > base.shape[1]
+        assert 'is_motile' in extended.columns
+        assert len(extended) == len(base)
+
+    def test_no_fragmentation_warning(self, timestep_with_dead):
+        with warnings.catch_warnings(record=True) as record:
+            warnings.simplefilter("always")
+            timestep_with_dead.to_dataframe(extra_columns=True)
+
+        assert not [
+            w for w in record
+            if issubclass(w.category, pd.errors.PerformanceWarning)
+        ]
+
+
+class TestReaderUsesFrameLabels:
+    """The reader resolves columns from each frame's own XML."""
+
+    def test_custom_column_reachable_through_reader(self, example_physicell_dir):
+        _, xml_path, _ = discover_physicell_timesteps(example_physicell_dir)[0]
+        timestep = read_physicell_timestep(xml_path)
+
+        labels = parse_physicell_xml(xml_path).extra['custom_labels']
+        expected = expand_cell_labels(labels)
+
+        data = timestep._load_cell_data()
+        assert data['columns'], "reader should pass labels through to the parser"
+        assert set(data['columns']) == set(expected.values())
+
+    def test_reader_emits_no_orientation_warning(self, example_physicell_dir):
+        _, xml_path, _ = discover_physicell_timesteps(example_physicell_dir)[0]
+
+        with warnings.catch_warnings(record=True) as record:
+            warnings.simplefilter("always")
+            read_physicell_timestep(xml_path)._load_cell_data()
+
+        assert not [
+            w for w in record
+            if 'Guessing cell matrix orientation' in str(w.message)
+        ]
+
+
+class TestVersionSingleSource:
+    """Issue #5: __version__ must not drift from distribution metadata."""
+
+    def test_version_matches_package_metadata(self):
+        assert spatialtissuepy.__version__ == importlib.metadata.version(
+            "spatialtissuepy"
+        )