diff --git a/.ed_planning/visreg/temporal-regularization/PLAN.md b/.ed_planning/visreg/temporal-regularization/PLAN.md
new file mode 100644
index 000000000..bc59e9c74
--- /dev/null
+++ b/.ed_planning/visreg/temporal-regularization/PLAN.md
@@ -0,0 +1,59 @@
+# Temporal regularization of DynaCLR embeddings — plan
+
+**Status (2026-07-27):** implementation, PR-review hardening, and local integration validation complete; remote CI remains.
+
+## Objective
+
+Keep NT-Xent for state separation, then add two losses on encoder features:
+
+- a shared next-state predictor, `L_pred`, to reward transition structure shared across cells;
+- a three-frame curvature loss, `L_curv`, to make distance along a track better reflect progression.
+
+The losses create temporal structure; they are not gated on straightness already being present in a frozen contrastive embedding.
+
+## Locked design
+
+- Emit fixed `K`-frame, `tau`-spaced sequences as row-major `B*K` tensors.
+- A valid sequence uses one exact `global_track_id` and marker. It stops at divisions and never switches siblings.
+- Apply temporal losses directly to encoder features; do not add a separate temporal head.
+- Stop-gradient the predictor target at `t+1`.
+- Keep temporal weights independently schedulable.
+- Default to track-consistent stochastic augmentation. Retain `independent` and `none` as explicit ablations.
+- Mask incomplete sequences while keeping every configured module in the distributed backward graph.
+
+## Milestones
+
+- [x] Sequence sampling, rectangular batching, and validity masks.
+- [x] Curvature loss and shared predictor.
+- [x] Lightning integration, schedules, and configuration.
+- [x] Local end-to-end smoke coverage.
+- [x] PR review hardening and regressions.
+- [x] Run the final combined local suite.
+- [ ] Confirm remote CI after push.
+- [ ] Compare biology and representation metrics against the NT-Xent baseline.
+
+## PR review resolutions
+
+| Priority | Risk | Resolution |
+| --- | --- | --- |
+| P1 | A lineage stencil could jump between siblings. | Match every frame by exact `global_track_id` and invalidate at division boundaries (`061d6e31`). |
+| P1 | Flattened `B*K` frames received independent random transforms. | Reuse one transform realization per track; expose `consistent`, `independent`, and `none` modes (`061d6e31`). |
+| P1 | An all-invalid rank left predictor parameters unused in DDP. | Forward the empty masked tensor through the predictor and backpropagate a graph-connected zero (`651d693c`). |
+| P2 | `positive_cell_source=self` skipped the lookup needed by sequence emission. | Build the lineage/timepoint lookup whenever sequences are enabled (`061d6e31`). |
+| P2 | Half-open HPI bins omitted a maximum on an exact boundary. | Share an edge builder that always adds a terminal bin (`f9426fab`). |
+| P3 | Zero, negative, or non-finite HPI widths were accepted. | Reject invalid widths during Pydantic validation (`f9426fab`). |
+
+## Acceptance checks
+
+- Exact-track and division-boundary sampler regressions pass.
+- Self-positive sequence emission passes.
+- Track-consistent flip regression passes.
+- Predictor parameters receive non-`None`, zero gradients for all-invalid batches.
+- HPI terminal-boundary and invalid-width regressions pass.
+- Dataset/datamodule, engine, witness-GMM, and final integration suites pass.
+
+Local combined validation passed on 2026-07-27. The two skipped inference-reproducibility tests require external HPC data and CUDA; there were no failures.
+
+## Follow-up experiment
+
+Train matched seeds for NT-Xent-only versus NT-Xent + `L_pred` + `L_curv`. Compare contrastive retrieval, collapse indicators, temporal prediction, per-track progression, and held-out biological separation before promoting the temporal objective to a default recipe.
diff --git a/.ed_planning/visreg/temporal-regularization/straightening-on-dynaclr.md b/.ed_planning/visreg/temporal-regularization/straightening-on-dynaclr.md
new file mode 100644
index 000000000..545136142
--- /dev/null
+++ b/.ed_planning/visreg/temporal-regularization/straightening-on-dynaclr.md
@@ -0,0 +1,63 @@
+# Straightening on DynaCLR — exploration
+
+**Status (2026-07-27):** promoted to implementation; review findings incorporated.
+
+## Question
+
+Can DynaCLR learn a common biological progression without erasing state organization?
+
+## Working hypothesis
+
+Straightening alone is a per-track geometric prior; it does not directly reward a transition shared across cells. The selected objective therefore combines:
+
+`L = L_NT-Xent + lambda_pred * L_pred + lambda_curv * L_curv`
+
+The shared predictor supplies the common-transition pressure, curvature regularization makes local motion easier to interpret, and NT-Xent resists constant-state collapse.
+
+## Implementation findings
+
+### Identity is stricter than lineage
+
+A lineage may contain a parent and multiple daughters. Sampling any same-lineage row at each timepoint can create a synthetic path that changes physical cells. Temporal stencils must match the anchor's exact `global_track_id`; a division makes the stencil invalid.
+
+### Augmentation is part of the temporal model
+
+Applying random transforms to flattened `B*K` frames independently injects artificial motion. The default must reuse one random realization across the `K` frames of each track while allowing different realizations across tracks. Independent and augmentation-free modes remain useful ablations.
+
+### Masked loss still has distributed semantics
+
+On a rank with no valid sequences, returning a zero connected only to encoder output leaves predictor parameters unused. Passing the empty input through the predictor produces zero gradients for its parameters and keeps DDP iteration state consistent.
+
+### Evaluation bins need explicit boundary semantics
+
+HPI loops use half-open intervals `[lo, hi)`. `np.arange(start, max + width, width)` still omits `max` when it is exactly a boundary. The edge builder must create one additional terminal edge, and widths must be finite and positive before arithmetic.
+
+## Decisions retained
+
+- Operate on encoder features, matching the representation being regularized.
+- Use three or more points for curvature; two points provide smoothing only.
+- Keep fixed-frame spacing initially for a small, auditable implementation.
+- Treat velocity-based phenotype splitting as a possible discovery signal, not automatically as representation damage.
+- Evaluate biology and collapse jointly; straightness by itself is not a success criterion.
+
+## Rejected shortcuts
+
+- Same-lineage sampling without exact track identity.
+- Independent stochastic transforms as the default temporal input.
+- Skipping configured modules on empty masked batches under DDP.
+- Using frozen-embedding straightness as a gate for whether training may induce temporal structure.
+
+## Evidence added by review hardening
+
+- Exact-track selection chooses one sibling consistently and invalidates parent-to-daughter stencils.
+- `self` positives can emit sequences without a missing lookup.
+- Repeated frames remain identical after track-consistent random flips.
+- Every predictor parameter receives a zero, non-`None` gradient on all-invalid batches.
+- Boundary-aligned maximum HPI values fall inside a bin; invalid widths fail at config construction.
+
+## Open empirical questions
+
+- Does the predictor learn biology rather than acquisition-time drift?
+- Which temporal weight schedule preserves contrastive retrieval best?
+- Is `K=3` sufficient, or do longer stencils improve robustness enough to justify their sampling cost?
+- Do consistent spatial transforms improve temporal metrics without weakening useful augmentation diversity?
diff --git a/applications/dynaclr/configs/training/debug/temporal-straightening-smoke.yml b/applications/dynaclr/configs/training/debug/temporal-straightening-smoke.yml
new file mode 100644
index 000000000..7f848eb26
--- /dev/null
+++ b/applications/dynaclr/configs/training/debug/temporal-straightening-smoke.yml
@@ -0,0 +1,82 @@
+# Smoke test: predictor + straightening temporal terms wired into DynaCLR.
+# Fast-dev run on the 2-FOV test parquet to confirm the sequence sampler emits
+# (B*K, C, Z, Y, X), the engine encodes it in one forward, and both temporal
+# losses (loss/curv, loss/pred) are computed without NaN/shape/DDP errors.
+#
+# Run:
+#   uv run --no-sync python -m dynaclr.cli fit -c \
+#     applications/dynaclr/configs/training/debug/temporal-straightening-smoke.yml
+
+seed_everything: 42
+
+trainer:
+  accelerator: gpu
+  devices: 1
+  precision: bf16-mixed
+  fast_dev_run: 4
+  logger: false
+  enable_checkpointing: false
+  enable_model_summary: false
+  use_distributed_sampler: false
+
+model:
+  class_path: dynaclr.engine.ContrastiveModule
+  init_args:
+    encoder:
+      class_path: viscy_models.contrastive.ContrastiveEncoder
+      init_args:
+        backbone: convnext_tiny
+        in_channels: 1
+        embedding_dim: 768
+        in_stack_depth: 1
+        stem_kernel_size: [1, 4, 4]
+        stem_stride: [1, 4, 4]
+        projection_dim: 32
+    loss_function:
+      class_path: viscy_models.contrastive.loss.NTXentLoss
+      init_args:
+        temperature: 0.2
+    lr: 0.00002
+    example_input_array_shape: [1, 1, 1, 160, 160]
+    straightening_loss:
+      class_path: viscy_models.contrastive.loss.TemporalStraighteningLoss
+    predictor:
+      class_path: viscy_models.contrastive.predictor.Predictor
+      init_args:
+        dim: 768
+    lambda_curv: 0.01
+    lambda_pred: 0.1
+
+data:
+  class_path: dynaclr.data.datamodule.MultiExperimentDataModule
+  init_args:
+    cell_index_path: /hpc/projects/organelle_phenotyping/models/collections/DynaCLR-2D-MIP-BagOfChannels-v3.parquet
+    z_window: 1
+    yx_patch_size: [192, 192]
+    final_yx_patch_size: [160, 160]
+    channels_per_sample: 1
+    positive_cell_source: lookup
+    positive_match_columns: [lineage_id]
+    tau_range: [0.5, 2.0]
+    tau_decay_rate: 2.0
+    emit_sequence: true
+    sequence_length: 3
+    sequence_tau_frames: 1
+    stratify_by: [perturbation, marker]
+    split_ratio: 0.8
+    batch_size: 16
+    num_workers: 1
+    seed: 42
+    normalizations:
+      - class_path: viscy_transforms.NormalizeSampled
+        init_args:
+          keys: [channel_0]
+          level: timepoint_statistics
+          subtrahend: mean
+          divisor: std
+    augmentations:
+      - class_path: viscy_transforms.BatchedRandFlipd
+        init_args:
+          keys: [channel_0]
+          spatial_axes: [1, 2]
+          prob: 0.5
diff --git a/applications/dynaclr/src/dynaclr/data/datamodule.py b/applications/dynaclr/src/dynaclr/data/datamodule.py
index 1c49aaf82..3cd05a273 100644
--- a/applications/dynaclr/src/dynaclr/data/datamodule.py
+++ b/applications/dynaclr/src/dynaclr/data/datamodule.py
@@ -11,7 +11,9 @@
 
 from __future__ import annotations
 
+import copy
 import logging
+from typing import Literal
 
 import numpy as np
 import pandas as pd
@@ -145,6 +147,12 @@ class MultiExperimentDataModule(LightningDataModule):
         Mapping from ``batch_key`` (used by classification heads) to
         dataframe column name.  E.g. ``{"gene_label": "condition"}``.
         Default: ``None``.
+    sequence_augment : {"consistent", "independent", "none"}
+        How training transforms are applied to emitted temporal sequences.
+        ``"consistent"`` (default) reuses one random transform realization
+        across all frames of each track. ``"independent"`` transforms the
+        flattened frames independently, and ``"none"`` applies normalization
+        and the final crop without stochastic augmentation.
     """
 
     def __init__(
@@ -190,6 +198,10 @@ def __init__(
         positive_match_columns: list[str] | None = None,
         positive_channel_source: str = "same",
         label_columns: dict[str, str] | None = None,
+        emit_sequence: bool = False,
+        sequence_length: int = 3,
+        sequence_tau_frames: int = 1,
+        sequence_augment: Literal["consistent", "independent", "none"] = "consistent",
         max_border_shift: int = -1,
         shuffle_val: bool = False,
         pin_memory: bool = True,
@@ -252,6 +264,14 @@ def __init__(
         self.positive_match_columns = positive_match_columns
         self.positive_channel_source = positive_channel_source
         self.label_columns = label_columns
+        self.emit_sequence = emit_sequence
+        self.sequence_length = sequence_length
+        self.sequence_tau_frames = sequence_tau_frames
+        if sequence_augment not in {"consistent", "independent", "none"}:
+            raise ValueError(
+                f"sequence_augment must be one of {{'consistent', 'independent', 'none'}}, got {sequence_augment!r}."
+            )
+        self.sequence_augment = sequence_augment
         self.max_border_shift = max_border_shift
         self.shuffle_val = shuffle_val
         self.pin_memory = pin_memory
@@ -318,6 +338,12 @@ def setup(self, stage: str | None = None) -> None:
             self._augmentation_transform = Compose(
                 self.normalizations + self.augmentations + [self._train_final_crop()]
             )
+            sequence_transforms = self.normalizations + [self._train_final_crop()]
+            if self.sequence_augment != "none":
+                sequence_transforms = self.normalizations + self.augmentations + [self._train_final_crop()]
+            # This transform has independent random state so resetting it for
+            # temporal consistency cannot perturb anchor/positive augmentation.
+            self._sequence_transform = Compose(copy.deepcopy(sequence_transforms))
 
             _logger.info(
                 "MultiExperimentDataModule setup: %d train anchors, %d val anchors",
@@ -419,6 +445,9 @@ def _setup_experiment_split(self, registry: ExperimentRegistry, cell_index_df: p
             positive_match_columns=self.positive_match_columns,
             positive_channel_source=self.positive_channel_source,
             label_columns=self.label_columns,
+            emit_sequence=self.emit_sequence,
+            sequence_length=self.sequence_length,
+            sequence_tau_frames=self.sequence_tau_frames,
         )
 
         if val_names:
@@ -445,6 +474,9 @@ def _setup_experiment_split(self, registry: ExperimentRegistry, cell_index_df: p
                 positive_match_columns=self.positive_match_columns,
                 positive_channel_source=self.positive_channel_source,
                 label_columns=self.label_columns,
+                emit_sequence=self.emit_sequence,
+                sequence_length=self.sequence_length,
+                sequence_tau_frames=self.sequence_tau_frames,
             )
 
     def _setup_fov_split(self, registry: ExperimentRegistry, cell_index_df: pd.DataFrame) -> None:
@@ -583,6 +615,9 @@ def _materialize_strings(df: pd.DataFrame) -> pd.DataFrame:
             positive_match_columns=self.positive_match_columns,
             positive_channel_source=self.positive_channel_source,
             label_columns=self.label_columns,
+            emit_sequence=self.emit_sequence,
+            sequence_length=self.sequence_length,
+            sequence_tau_frames=self.sequence_tau_frames,
         )
 
         if not val_tracks.empty:
@@ -602,6 +637,9 @@ def _materialize_strings(df: pd.DataFrame) -> pd.DataFrame:
                 positive_match_columns=self.positive_match_columns,
                 positive_channel_source=self.positive_channel_source,
                 label_columns=self.label_columns,
+                emit_sequence=self.emit_sequence,
+                sequence_length=self.sequence_length,
+                sequence_tau_frames=self.sequence_tau_frames,
             )
 
     # ------------------------------------------------------------------
@@ -729,6 +767,54 @@ def _train_final_crop(self) -> BatchedCenterSpatialCropd:
             ),
         )
 
+    def _transform_sequence_consistently(
+        self,
+        patch: Tensor,
+        norm_meta: list | None,
+        extra: dict[str, Tensor] | None,
+    ) -> Tensor:
+        """Apply one stochastic transform realization per temporal track.
+
+        The dataset flattens sequences row-major as ``B * K``. Transforming
+        each temporal slot as a batch of ``B`` and restoring the same random
+        state before every slot gives sample ``i`` identical random parameters
+        at all ``K`` timepoints while retaining independent choices across
+        samples.
+        """
+        k = self.sequence_length
+        if patch.shape[0] % k:
+            raise ValueError(f"Sequence batch size {patch.shape[0]} is not divisible by sequence_length={k}.")
+
+        seed = int(torch.randint(0, torch.iinfo(torch.int32).max, ()).item())
+        slot_outputs: list[Tensor] = []
+        for slot in range(k):
+            slot_indices = torch.arange(slot, patch.shape[0], k, device=patch.device)
+            slot_norm_meta = norm_meta[slot::k] if norm_meta is not None else None
+            slot_extra = None
+            if extra is not None:
+                slot_extra = {
+                    name: value.index_select(0, slot_indices)
+                    if isinstance(value, Tensor) and value.shape[0] == patch.shape[0]
+                    else value
+                    for name, value in extra.items()
+                }
+
+            devices = [patch.device] if patch.is_cuda else []
+            with torch.random.fork_rng(devices=devices):
+                torch.manual_seed(seed)
+                self._sequence_transform.set_random_state(seed=seed)
+                slot_outputs.append(
+                    _transform_channel_wise(
+                        transform=self._sequence_transform,
+                        channel_names=self._channel_names,
+                        patch=patch.index_select(0, slot_indices),
+                        norm_meta=slot_norm_meta,
+                        extra=slot_extra,
+                    )
+                )
+
+        return torch.stack(slot_outputs, dim=1).reshape(-1, *slot_outputs[0].shape[1:])
+
     def on_after_batch_transfer(self, batch, dataloader_idx: int):
         """Apply normalizations, augmentations, final crop, and ChannelDropout.
 
@@ -778,9 +864,7 @@ def on_after_batch_transfer(self, batch, dataloader_idx: int):
             batch.pop("anchor_meta", None)
             return batch
 
-        transform = self._augmentation_transform
-
-        for key in ["anchor", "positive", "negative"]:
+        for key in ["anchor", "positive", "negative", "sequence"]:
             if key in batch:
                 norm_meta_key = f"{key}_norm_meta"
                 norm_meta = batch.get(norm_meta_key)
@@ -805,13 +889,21 @@ def on_after_batch_transfer(self, batch, dataloader_idx: int):
                                 device=batch[key].device,
                             )
                         }
-                transformed = _transform_channel_wise(
-                    transform=transform,
-                    channel_names=self._channel_names,
-                    patch=batch[key],
-                    norm_meta=norm_meta,
-                    extra=extra,
-                )
+                if key == "sequence" and self.sequence_augment == "consistent":
+                    transformed = self._transform_sequence_consistently(
+                        patch=batch[key],
+                        norm_meta=norm_meta,
+                        extra=extra,
+                    )
+                else:
+                    transform = self._sequence_transform if key == "sequence" else self._augmentation_transform
+                    transformed = _transform_channel_wise(
+                        transform=transform,
+                        channel_names=self._channel_names,
+                        patch=batch[key],
+                        norm_meta=norm_meta,
+                        extra=extra,
+                    )
                 batch[key] = transformed
                 if norm_meta_key in batch:
                     del batch[norm_meta_key]
diff --git a/applications/dynaclr/src/dynaclr/data/dataset.py b/applications/dynaclr/src/dynaclr/data/dataset.py
index a682b744e..080f3ce97 100644
--- a/applications/dynaclr/src/dynaclr/data/dataset.py
+++ b/applications/dynaclr/src/dynaclr/data/dataset.py
@@ -199,6 +199,9 @@ def __init__(
         positive_match_columns: list[str] | None = None,
         positive_channel_source: str = "same",
         label_columns: dict[str, str] | None = None,
+        emit_sequence: bool = False,
+        sequence_length: int = 3,
+        sequence_tau_frames: int = 1,
     ) -> None:
         if ts is None:
             raise ImportError(
@@ -241,6 +244,28 @@ def __init__(
         self.positive_match_columns = positive_match_columns if positive_match_columns is not None else ["lineage_id"]
         self.positive_channel_source = positive_channel_source
 
+        # Temporal sequence emission (for straightening / predictor losses):
+        # emit K consecutive same-track same-marker frames at fixed frame stride.
+        self.emit_sequence = emit_sequence
+        self.sequence_length = sequence_length
+        self.sequence_tau_frames = sequence_tau_frames
+        if emit_sequence:
+            if sequence_length < 3:
+                raise ValueError(f"sequence_length must be >= 3 for curvature, got {sequence_length}")
+            if sequence_tau_frames < 1:
+                raise ValueError(f"sequence_tau_frames must be >= 1, got {sequence_tau_frames}")
+            if self._channel_mode != "from_index":
+                raise ValueError(
+                    "emit_sequence requires bag-of-channels mode (channels_per_sample=1); "
+                    f"got channel_mode={self._channel_mode!r}. Per-marker sequences need one "
+                    "channel per sample for a well-defined marker filter."
+                )
+            if "lineage_id" not in self.positive_match_columns:
+                raise ValueError(
+                    "emit_sequence requires 'lineage_id' in positive_match_columns "
+                    "so the lineage-timepoint lookup is built."
+                )
+
         self._label_encoders: dict[str, tuple[str, dict[str, int]]] = {}
         if label_columns:
             for batch_key, col in label_columns.items():
@@ -265,7 +290,8 @@ def __init__(
     def _build_match_lookup(self) -> None:
         """Build lookup structures for O(1) positive candidate lookup.
 
-        For ``positive_cell_source="self"``, no lookup is needed.
+        For ``positive_cell_source="self"``, no positive lookup is needed, but
+        sequence emission still needs the lineage/timepoint lookup.
 
         For temporal mode (``"lineage_id"`` in ``positive_match_columns``),
         builds ``_lineage_timepoints``:
@@ -274,7 +300,7 @@ def _build_match_lookup(self) -> None:
         For generic column-match mode, builds ``_match_lookup``:
         ``{match_key_tuple: [row_indices_in_tracks]}``.
         """
-        if self.positive_cell_source == "self":
+        if self.positive_cell_source == "self" and not self.emit_sequence:
             return
 
         tracks = self.index.tracks
@@ -349,6 +375,7 @@ def _cache_columns(df: pd.DataFrame, columns: list[str]) -> dict:
         hot_cols: set[str] = {
             "channel_name",
             "experiment",
+            "global_track_id",
             "lineage_id",
             "t",
             "marker",
@@ -456,6 +483,20 @@ def __getitems__(self, indices: list[int]) -> dict:
                 sample["positive"] = positive_patches
                 sample["positive_norm_meta"] = positive_norms
                 sample["positive_meta"] = self._extract_meta(positive_rows)
+
+            if self.emit_sequence:
+                seq_track_indices, seq_valid = self._sample_sequence_indices(anchor_positions=indices)
+                tr_chan_arr = self._tr_arrays["channel_name"]
+                seq_forced_channel_names = [[tr_chan_arr[i]] for i in seq_track_indices]
+                seq_patches, seq_norms = self._slice_patches(
+                    self._tr_arrays, seq_track_indices, seq_forced_channel_names
+                )
+                seq_rows = self.index.tracks.iloc[seq_track_indices].reset_index(drop=True)
+                sample["sequence"] = seq_patches  # (B*K, C, Z, Y, X), row-major over (B, K)
+                sample["sequence_norm_meta"] = seq_norms
+                sample["sequence_meta"] = self._extract_meta(seq_rows)
+                sample["sequence_valid"] = torch.from_numpy(seq_valid)
+                sample["sequence_length"] = self.sequence_length
         else:
             # Build per-sample index dicts via NumPy column arrays (no .iterrows).
             all_cols = list(ULTRACK_INDEX_COLUMNS) + [
@@ -633,6 +674,87 @@ def _sample_positive_indices_temporal(self, anchor_positions: list[int]) -> np.n
 
         return pos_track_indices
 
+    def _sample_sequence_indices(
+        self,
+        anchor_positions: list[int],
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Sample K consecutive same-lineage same-marker frames per anchor.
+
+        Builds a fixed stencil ``[t0, t0+tau, ..., t0+(K-1)*tau]`` (frames) for
+        the anchor's exact ``global_track_id`` and marker. Sequences never cross
+        a division boundary or switch between sibling tracks. When every
+        stencil timepoint has an exact-track candidate, the sequence is valid;
+        otherwise it is marked invalid and filled with the first available
+        exact-track row K times (placeholder, keeping the tensor rectangular;
+        the loss drops it via ``valid``).
+
+        Frames are grouped row-major per sample: ``[s0_t0, s0_t1, ..., s0_t(K-1),
+        s1_t0, ...]`` so the engine can ``view(B, K, D)`` to recover grouping.
+
+        Parameters
+        ----------
+        anchor_positions : list[int]
+            Positional indices into ``valid_anchors`` for the batch.
+
+        Returns
+        -------
+        tuple[np.ndarray, np.ndarray]
+            ``seq_track_indices`` of shape ``(B * K,)`` (positional indices into
+            ``self.index.tracks``) and ``seq_valid`` of shape ``(B,)`` (bool).
+        """
+        rng = self._rng
+        exp_arr = self._va_arrays["experiment"]
+        lid_arr = self._va_arrays["lineage_id"]
+        t_arr = self._va_arrays["t"]
+        anchor_marker_arr = self._va_arrays["marker"]
+        anchor_track_arr = self._va_arrays["global_track_id"]
+        tr_marker_arr = self._tr_arrays["marker"]
+        tr_track_arr = self._tr_arrays["global_track_id"]
+        lt_map = self._lineage_timepoints
+        k = self.sequence_length
+        tau = self.sequence_tau_frames
+
+        b = len(anchor_positions)
+        seq_track_indices = np.empty(b * k, dtype=np.int64)
+        seq_valid = np.zeros(b, dtype=bool)
+
+        for i, ai in enumerate(anchor_positions):
+            exp_name = str(exp_arr[ai])
+            lineage_id = str(lid_arr[ai])
+            anchor_t = int(t_arr[ai])
+            anchor_marker = anchor_marker_arr[ai]
+            anchor_track = anchor_track_arr[ai]
+            timepoints = lt_map.get((exp_name, lineage_id))
+
+            rows: list[int] | None = None
+            picked: list[int] = []
+            if timepoints is not None:
+                wanted = [anchor_t + j * tau for j in range(k)]
+                for wt in wanted:
+                    cands = timepoints.get(wt)
+                    if not cands:
+                        break
+                    idx_arr = np.asarray(cands, dtype=np.int64)
+                    mask = (tr_marker_arr[idx_arr] == anchor_marker) & (tr_track_arr[idx_arr] == anchor_track)
+                    filtered = idx_arr[mask]
+                    if len(filtered) == 0:
+                        break
+                    picked.append(int(filtered[rng.integers(len(filtered))]))
+                if len(picked) == k:
+                    rows = picked
+
+            if rows is None:
+                # Reuse the anchor-time exact-track row when it was found;
+                # otherwise row 0 is only a shape-preserving fallback. The
+                # invalid sample is masked from both temporal losses.
+                placeholder = picked[0] if picked else 0
+                rows = [placeholder] * k
+            else:
+                seq_valid[i] = True
+            seq_track_indices[i * k : (i + 1) * k] = rows
+
+        return seq_track_indices, seq_valid
+
     # ------------------------------------------------------------------
     # Patch extraction (tensorstore I/O)
     # ------------------------------------------------------------------
diff --git a/applications/dynaclr/src/dynaclr/engine.py b/applications/dynaclr/src/dynaclr/engine.py
index 924ee9b5c..ba9f70d6b 100644
--- a/applications/dynaclr/src/dynaclr/engine.py
+++ b/applications/dynaclr/src/dynaclr/engine.py
@@ -51,12 +51,32 @@ def __init__(
         freeze_backbone: bool = False,
         projection: nn.Module | None = None,
         auxiliary_heads: dict[str, BaseHead] | None = None,
+        straightening_loss: nn.Module | None = None,
+        predictor: nn.Module | None = None,
+        lambda_curv: float = 0.0,
+        lambda_pred: float = 0.0,
+        curv_schedule: Literal["cosine", "constant"] = "constant",
+        pred_schedule: Literal["cosine", "constant"] = "constant",
+        curv_weight_start: float = 0.0,
+        pred_weight_start: float = 0.0,
+        temporal_warmup_epochs: int = 50,
     ) -> None:
         super().__init__()
         self.model = encoder
         if projection is not None:
             self.model.projection = projection
         self.loss_function = loss_function
+        self.straightening_loss = straightening_loss
+        self.predictor = predictor
+        self.lambda_curv = lambda_curv
+        self.lambda_pred = lambda_pred
+        self.curv_schedule = curv_schedule
+        self.pred_schedule = pred_schedule
+        self.curv_weight_start = curv_weight_start
+        self.pred_weight_start = pred_weight_start
+        self.temporal_warmup_epochs = temporal_warmup_epochs
+        self._curv_weight = curv_weight_start if curv_schedule == "cosine" else lambda_curv
+        self._pred_weight = pred_weight_start if pred_schedule == "cosine" else lambda_pred
         self.lr = lr
         self.schedule = schedule
         self.log_batches_per_epoch = log_batches_per_epoch
@@ -93,6 +113,58 @@ def on_train_epoch_start(self) -> None:  # noqa: D102
         for head in self.auxiliary_heads.values():
             head.step(self.current_epoch)
             self.log(f"hparams/loss_weight/{head.head_name}", head.get_weight())
+        if self.straightening_loss is not None:
+            if self.curv_schedule == "cosine":
+                self._curv_weight = cosine_anneal(
+                    self.curv_weight_start, self.lambda_curv, self.current_epoch, self.temporal_warmup_epochs
+                )
+            self.log("hparams/loss_weight/curv", self._curv_weight)
+        if self.predictor is not None:
+            if self.pred_schedule == "cosine":
+                self._pred_weight = cosine_anneal(
+                    self.pred_weight_start, self.lambda_pred, self.current_epoch, self.temporal_warmup_epochs
+                )
+            self.log("hparams/loss_weight/pred", self._pred_weight)
+
+    def _temporal_losses(self, batch: TripletSample, stage: Literal["train", "val"]) -> Tensor:
+        """Compute weighted straightening + predictor losses on the sequence batch.
+
+        Encodes the ``(B*K, C, Z, Y, X)`` sequence in ONE forward, reshapes to
+        ``(B, K, 768)``, and applies both temporal terms directly to the encoder
+        features ``z`` (no head). Returns the weighted sum (zero tensor when no
+        sequence / no temporal modules are configured).
+        """
+        if "sequence" not in batch or (self.straightening_loss is None and self.predictor is None):
+            return torch.zeros((), device=self.device)
+        seq = batch["sequence"]  # (B*K, C, Z, Y, X)
+        k = int(batch["sequence_length"])
+        b = seq.size(0) // k
+        z, _ = self(seq)  # ONE forward -> (B*K, 768)
+        z = z.view(b, k, -1)  # (B, K, 768)
+        valid = batch.get("sequence_valid")
+        on_step = stage == "train"
+        total = torch.zeros((), device=self.device)
+        if self.straightening_loss is not None:
+            curv = self.straightening_loss(z, valid=valid)
+            total = total + self._curv_weight * curv
+            self.log(f"loss/curv/{stage}", curv, on_step=on_step, on_epoch=not on_step, sync_dist=True, batch_size=b)
+        if self.predictor is not None:
+            src = z[:, :-1, :].reshape(-1, z.size(-1))  # (B*(K-1), 768) z_t
+            tgt = z[:, 1:, :].reshape(-1, z.size(-1)).detach()  # sg(z_{t+1})
+            if valid is not None:
+                mask = valid.view(b, 1).expand(b, k - 1).reshape(-1)
+                src, tgt = src[mask], tgt[mask]
+            pred = self.predictor(src)
+            if pred.numel() == 0:
+                # Keep predictor parameters in the graph on ranks whose batch
+                # contains no valid sequence, avoiding DDP unused-parameter
+                # failures on the following iteration.
+                l_pred = pred.sum() * 0.0
+            else:
+                l_pred = F.mse_loss(pred, tgt)
+            total = total + self._pred_weight * l_pred
+            self.log(f"loss/pred/{stage}", l_pred, on_step=on_step, on_epoch=not on_step, sync_dist=True, batch_size=b)
+        return total
 
     def on_fit_start(self) -> None:  # noqa: D102
         if self.freeze_backbone:
@@ -284,6 +356,7 @@ def training_step(self, batch: TripletSample, batch_idx: int) -> Tensor:  # noqa
             stage="train",
         )
         loss = loss + self._run_auxiliary_heads(anchor_features, batch, "train")
+        loss = loss + self._temporal_losses(batch, "train")
         return loss
 
     def on_train_epoch_end(self) -> None:  # noqa: D102
@@ -316,6 +389,7 @@ def validation_step(self, batch: TripletSample, batch_idx: int) -> Tensor:  # no
             stage="val",
         )
         loss = loss + self._run_auxiliary_heads(anchor_features, batch, "val")
+        loss = loss + self._temporal_losses(batch, "val")
         n = self.log_embeddings_every_n_epochs
         if n is not None and self.current_epoch % n == 0 and not self.trainer.sanity_checking:
             self._embedding_outputs.append((anchor_features.detach().cpu(), batch.get("anchor_meta", [])))
diff --git a/applications/dynaclr/src/dynaclr/evaluation/evaluate_config.py b/applications/dynaclr/src/dynaclr/evaluation/evaluate_config.py
index 619ad0d99..4676d5aa8 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/evaluate_config.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/evaluate_config.py
@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import math
 from typing import Literal
 
 from pydantic import BaseModel, model_validator
@@ -322,6 +323,24 @@ class WitnessGmmLabelsConfig(BaseModel):
         the witness (bounds kernel cost). None = use all. Default: 5000.
     random_seed : int
         Seed for reference subsampling and the GMM. Default: 42.
+    mmd_hpi_bin_hours : float or None
+        If set, additionally compute MMD²(control, condition) **per time bin** of
+        this width (hours post perturbation) and write an MMD-vs-HPI diagnostic
+        plot — the population divergence kinetics, plus a control-vs-control null
+        band. Requires an ``hours_post_perturbation`` obs column. None disables
+        the diagnostic. Default: None.
+    witness_time_bin_hours : float or None
+        If set, score the witness with **time-matched references**: each cell is
+        scored against control/perturbed reference cells drawn from its own
+        ``hours_post_perturbation`` bin of this width, rather than against a single
+        pooled all-timepoint reference. The DynaCLR embedding carries a strong
+        time/culture axis (uninfected cells drift over a long timelapse), so a
+        pooled reference leaks that axis into the witness score and the labels —
+        worst for weak channels. Time-matching cancels the shared time component
+        so the witness axis reflects perturbation, not culture-time. A single
+        global bandwidth (median heuristic on the pooled reference) is kept so
+        scores stay comparable across bins. None = pooled reference (original
+        behavior). Requires an ``hours_post_perturbation`` obs column. Default: None.
     """
 
     experiments: list[WitnessGmmExperiment]
@@ -336,7 +355,11 @@ class WitnessGmmLabelsConfig(BaseModel):
     mmd_n_permutations: int = 1000
     bandwidth: float | None = None
     max_reference_cells: int | None = 5000
+    witness_time_bin_hours: float | None = None
+    gate: str = "gmm"
+    control_fp_target: float = 0.05
     random_seed: int = 42
+    mmd_hpi_bin_hours: float | None = None
 
     @model_validator(mode="after")
     def _validate(self) -> "WitnessGmmLabelsConfig":
@@ -347,10 +370,18 @@ def _validate(self) -> "WitnessGmmLabelsConfig":
             raise ValueError(f"class_map must define {sorted(missing)} (got keys {sorted(self.class_map)})")
         if not 0.0 < self.gmm_pos_threshold <= 1.0:
             raise ValueError(f"gmm_pos_threshold must be in (0, 1], got {self.gmm_pos_threshold}")
+        if self.gate not in ("gmm", "control_anchored"):
+            raise ValueError(f"gate must be 'gmm' or 'control_anchored', got {self.gate!r}")
+        if not 0.0 < self.control_fp_target < 1.0:
+            raise ValueError(f"control_fp_target must be in (0, 1), got {self.control_fp_target}")
         if not 0.0 < self.mmd_pvalue_threshold <= 1.0:
             raise ValueError(f"mmd_pvalue_threshold must be in (0, 1], got {self.mmd_pvalue_threshold}")
         if self.annotation_format not in ("csv", "parquet"):
             raise ValueError(f"annotation_format must be 'csv' or 'parquet', got {self.annotation_format!r}")
+        for field in ("witness_time_bin_hours", "mmd_hpi_bin_hours"):
+            width = getattr(self, field)
+            if width is not None and (not math.isfinite(width) or width <= 0):
+                raise ValueError(f"{field} must be a finite positive number, got {width}")
         return self
 
 
diff --git a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels.py b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels.py
index 6e05b0559..92f348c90 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels.py
@@ -39,12 +39,18 @@ class vocabulary — a file indistinguishable from a hand annotation.
 
 from dynaclr.evaluation.linear_classifiers.witness_gmm_plots import (
     plot_mmd_null,
+    plot_mmd_vs_hpi,
     plot_remodeling_vs_time,
     plot_witness_gmm,
 )
 from viscy_utils.cli_utils import load_config
 from viscy_utils.evaluation.mmd import median_heuristic, mmd_permutation_test, subsample, witness_function
-from viscy_utils.evaluation.witness_gmm import fit_gmm_labels
+from viscy_utils.evaluation.witness_gmm import (
+    ControlAnchoredResult,
+    _gaussian_pdf,
+    fit_control_anchored_labels,
+    fit_gmm_labels,
+)
 
 if TYPE_CHECKING:
     from dynaclr.evaluation.evaluate_config import WitnessGmmExperiment, WitnessGmmLabelsConfig
@@ -192,6 +198,9 @@ class _MarkerScores:
     scores: np.ndarray  # witness score per cell
     conditions: np.ndarray  # obs[condition_column] as array
     cond_mmd: dict  # condition -> _MmdResult (raw MMD² + p-value + null)
+    # Optional MMD-vs-HPI kinetics: condition -> list of (hpi_bin_center, mmd2, p).
+    # Key "__control_null__" holds the control-vs-control baseline (should be ~0).
+    hpi_mmd: dict | None = None
 
     @property
     def cond_pvalues(self) -> dict:
@@ -254,8 +263,16 @@ def compute_marker_scores(
 
     X_ref = subsample(X_all[control_mask], config.max_reference_cells, rng)
     Y_ref = subsample(X_all[perturbed_mask], config.max_reference_cells, rng)
+    # Global bandwidth from the pooled reference — kept even in time-matched mode
+    # so witness scores are on one comparable scale across timepoints.
     bandwidth = config.bandwidth if config.bandwidth is not None else median_heuristic(X_ref, Y_ref)
-    scores = witness_function(X_all, X_ref, Y_ref, bandwidth=bandwidth)
+
+    if config.witness_time_bin_hours is None:
+        scores = witness_function(X_all, X_ref, Y_ref, bandwidth=bandwidth)
+    else:
+        scores = _time_matched_witness_scores(
+            X_all, obs, control_mask, perturbed_mask, X_ref, Y_ref, bandwidth, config, rng
+        )
 
     # Per-condition MMD significance: is this condition's cloud distinct from the
     # control reference? Raw p-values here; FDR-corrected run-wide by the caller.
@@ -274,7 +291,145 @@ def compute_marker_scores(
         )
         cond_mmd[cond] = _MmdResult(mmd2=float(mmd2), p_value=float(p_value), null=np.asarray(null))
 
-    return _MarkerScores(obs, control_mask, perturbed_mask, scores, conditions, cond_mmd)
+    hpi_mmd = None
+    if config.mmd_hpi_bin_hours is not None:
+        hpi_mmd = _compute_hpi_mmd(X_all, obs, control_mask, perturbed_mask, conditions, bandwidth, config, rng)
+
+    return _MarkerScores(obs, control_mask, perturbed_mask, scores, conditions, cond_mmd, hpi_mmd)
+
+
+def _time_matched_witness_scores(
+    X_all: np.ndarray,
+    obs: pd.DataFrame,
+    control_mask: np.ndarray,
+    perturbed_mask: np.ndarray,
+    X_ref_pooled: np.ndarray,
+    Y_ref_pooled: np.ndarray,
+    bandwidth: float,
+    config: WitnessGmmLabelsConfig,
+    rng: np.random.Generator,
+) -> np.ndarray:
+    """Witness score per cell using per-HPI-bin (time-matched) references.
+
+    For each ``witness_time_bin_hours``-wide window, cells in that window are
+    scored against control/perturbed reference cells **from the same window**, so
+    the witness axis reflects perturbation rather than culture-time (the embedding
+    has a strong time axis; see :func:`_compute_hpi_mmd`). Bins lacking enough
+    reference cells on either side fall back to the pooled reference, so every
+    cell is scored. The bandwidth is the shared global value for a comparable
+    scale across bins.
+    """
+    if "hours_post_perturbation" not in obs.columns:
+        _logger.warning("witness_time_bin_hours set but no hours_post_perturbation column; using pooled reference.")
+        return witness_function(X_all, X_ref_pooled, Y_ref_pooled, bandwidth=bandwidth)
+
+    hpi = obs["hours_post_perturbation"].to_numpy(dtype=float)
+    width = config.witness_time_bin_hours
+    scores = np.full(len(obs), np.nan, dtype=np.float64)
+    finite = np.isfinite(hpi)
+    edges = _hpi_bin_edges(hpi, width)
+    min_ref = 5
+    n_fallback = 0
+    for lo, hi in zip(edges[:-1], edges[1:]):
+        in_bin = finite & (hpi >= lo) & (hpi < hi)
+        if not in_bin.any():
+            continue
+        ctrl_bin = X_all[control_mask & in_bin]
+        pert_bin = X_all[perturbed_mask & in_bin]
+        if len(ctrl_bin) >= min_ref and len(pert_bin) >= min_ref:
+            xr = subsample(ctrl_bin, config.max_reference_cells, rng)
+            yr = subsample(pert_bin, config.max_reference_cells, rng)
+        else:
+            xr, yr = X_ref_pooled, Y_ref_pooled  # sparse bin → pooled fallback
+            n_fallback += int(in_bin.sum())
+        scores[in_bin] = witness_function(X_all[in_bin], xr, yr, bandwidth=bandwidth)
+    # Cells with non-finite HPI (never assigned) get the pooled score.
+    missing = np.isnan(scores)
+    if missing.any():
+        scores[missing] = witness_function(X_all[missing], X_ref_pooled, Y_ref_pooled, bandwidth=bandwidth)
+    if n_fallback:
+        _logger.info("Time-matched witness: %d cells in sparse bins used the pooled reference.", n_fallback)
+    return scores
+
+
+def _hpi_bin_edges(hpi: np.ndarray, width: float) -> np.ndarray:
+    """Return half-open bin edges that include the maximum finite HPI.
+
+    An extra terminal edge is required when the maximum lies exactly on a bin
+    boundary; otherwise ``[lo, hi)`` loops silently omit those cells.
+    """
+    finite = np.asarray(hpi, dtype=float)
+    finite = finite[np.isfinite(finite)]
+    if finite.size == 0:
+        return np.empty(0, dtype=float)
+    start = np.floor(finite.min() / width) * width
+    n_bins = int(np.floor((finite.max() - start) / width)) + 1
+    return start + np.arange(n_bins + 1, dtype=float) * width
+
+
+def _compute_hpi_mmd(
+    X_all: np.ndarray,
+    obs: pd.DataFrame,
+    control_mask: np.ndarray,
+    perturbed_mask: np.ndarray,
+    conditions: np.ndarray,
+    bandwidth: float,
+    config: WitnessGmmLabelsConfig,
+    rng: np.random.Generator,
+) -> dict:
+    """Time-matched MMD²(control, condition) per HPI bin — infection kinetics.
+
+    For each ``mmd_hpi_bin_hours``-wide window of ``hours_post_perturbation``,
+    tests each condition's cells against the **control cells in the SAME window**
+    (not the pooled all-timepoint control reference). This time-matching is
+    essential: the DynaCLR embedding carries a strong time/culture axis (uninfected
+    cells drift substantially over a 36 h timelapse), so comparing a bin's
+    perturbed cells to a time-pooled control conflates infection with that time
+    axis. Matching control and perturbed within the bin cancels the shared time
+    component and isolates the infection difference. Also emits a
+    control-vs-control null (random half-split of the bin's control cells) under
+    ``"__control_null__"`` as the no-difference floor. Bins with too few cells on
+    either side are skipped. Returns ``{key: [(hpi_center, mmd2, p), ...]}``.
+    """
+    if "hours_post_perturbation" not in obs.columns:
+        _logger.warning("mmd_hpi_bin_hours set but no hours_post_perturbation column; skipping HPI-MMD.")
+        return {}
+    hpi = obs["hours_post_perturbation"].to_numpy(dtype=float)
+    width = config.mmd_hpi_bin_hours
+    edges = _hpi_bin_edges(hpi, width)
+    if edges.size == 0:
+        return {}
+
+    def _mmd(a: np.ndarray, b: np.ndarray) -> tuple[float, float]:
+        mmd2, p, _null = mmd_permutation_test(
+            subsample(a, config.max_reference_cells, rng),
+            subsample(b, config.max_reference_cells, rng),
+            n_permutations=config.mmd_n_permutations,
+            bandwidth=bandwidth,
+            seed=config.random_seed,
+        )
+        return float(mmd2), float(p)
+
+    out: dict = {}
+    for lo, hi in zip(edges[:-1], edges[1:]):
+        in_bin = (hpi >= lo) & (hpi < hi)
+        center = float(lo + width / 2)
+        ctrl_bin = X_all[control_mask & in_bin]
+        if len(ctrl_bin) < 5:
+            continue  # need a time-matched control reference for this bin
+        for cond in pd.unique(conditions[perturbed_mask]):
+            cells = X_all[perturbed_mask & (conditions == cond) & in_bin]
+            if len(cells) < 5:
+                continue
+            mmd2, p = _mmd(ctrl_bin, cells)
+            out.setdefault(str(cond), []).append((center, mmd2, p))
+        # Control-vs-control null: split this bin's control cells in half.
+        if len(ctrl_bin) >= 10:
+            perm = rng.permutation(len(ctrl_bin))
+            half = len(ctrl_bin) // 2
+            mmd2, p = _mmd(ctrl_bin[perm[:half]], ctrl_bin[perm[half:]])
+            out.setdefault("__control_null__", []).append((center, mmd2, p))
+    return out
 
 
 def label_marker(
@@ -314,6 +469,7 @@ def label_marker(
     time_col = "hours_post_perturbation" if "hours_post_perturbation" in obs.columns else "t"
     time_values = obs[time_col].to_numpy() if time_col in obs.columns else None
 
+    control_scores_all = scores[control_mask]
     cond_gmm: dict = {}
     cond_scores: dict = {}
     cond_time: dict = {}
@@ -325,16 +481,27 @@ def label_marker(
         if cond not in significant_conditions:
             _logger.warning("MMD not significant (FDR) for condition %r; no positives labeled.", cond)
             continue
-        res = fit_gmm_labels(scores[cond_mask], pos_threshold=config.gmm_pos_threshold, random_state=config.random_seed)
-        cond_gmm[cond] = res
         cond_scores[cond] = scores[cond_mask]
         if time_values is not None:
             cond_time[cond] = time_values[cond_mask]
-        if not res.separated:
-            _logger.warning("GMM unimodal for condition %r; no positives labeled.", cond)
-            continue
-        any_separated = True
         cond_idx = np.flatnonzero(cond_mask)
+        if config.gate == "control_anchored":
+            # Baseline frozen to control; label the excess over baseline. Never
+            # abstains — right for graded, subtle shifts (weak channels).
+            res = fit_control_anchored_labels(
+                scores[cond_mask], control_scores_all, control_fp_target=config.control_fp_target
+            )
+            cond_gmm[cond] = res
+            any_separated = True
+        else:
+            res = fit_gmm_labels(
+                scores[cond_mask], pos_threshold=config.gmm_pos_threshold, random_state=config.random_seed
+            )
+            cond_gmm[cond] = res
+            if not res.separated:
+                _logger.warning("GMM unimodal for condition %r; no positives labeled.", cond)
+                continue
+            any_separated = True
         idx = cond_idx[res.hard_label == 1]
         labels[idx] = pos_label
         gmm_posterior[idx] = res.posterior[res.hard_label == 1]
@@ -359,13 +526,19 @@ def label_marker(
     # drives the remodeling-vs-time control line as an empirical FALSE-POSITIVE
     # rate, not a hardcoded zero. NaN posterior for controls when no GMM separated.
     control_pos = np.zeros(control_mask.sum(), dtype=bool)
-    control_post = np.full(control_mask.sum(), np.nan, dtype=float)
-    ctrl_X = scores[control_mask].reshape(-1, 1)
+    ctrl_s = scores[control_mask]
     for res in cond_gmm.values():
+        if isinstance(res, ControlAnchoredResult):
+            # Control-anchored: score controls under the fitted mixture, threshold
+            # at the calibrated cut (FP ≈ control_fp_target by construction).
+            base = res.pi_baseline * _gaussian_pdf(ctrl_s, res.mu_c, res.sigma_c)
+            rem = (1 - res.pi_baseline) * _gaussian_pdf(ctrl_s, res.mu_r, res.sigma_r)
+            post = rem / (base + rem + 1e-300)
+            control_pos |= post >= res.threshold
+            continue
         if not res.separated:
             continue
-        post = res.gmm.predict_proba(ctrl_X)[:, res.remod_component]
-        control_post = post if np.isnan(control_post).all() else np.maximum(control_post, post)
+        post = res.gmm.predict_proba(ctrl_s.reshape(-1, 1))[:, res.remod_component]
         control_pos |= post >= config.gmm_pos_threshold
     return _MarkerLabels(
         frame.reset_index(drop=True),
@@ -378,7 +551,7 @@ def label_marker(
     )
 
 
-def generate_witness_gmm_annotation(config: WitnessGmmLabelsConfig) -> Path:
+def generate_witness_gmm_annotation(config: WitnessGmmLabelsConfig) -> Path | None:
     """Run Stage A end to end and write the annotation file plus diagnostics.
 
     Loads each experiment's embeddings zarr, pools per marker, labels via
@@ -473,6 +646,14 @@ class vs time per condition, for every labeled marker.
                 condition=str(cond),
                 output_path=plots_dir / f"mmd_null_{marker}_{cond}.png",
             )
+        # MMD-vs-HPI kinetics (population divergence over time) when enabled.
+        if ms.hpi_mmd:
+            plot_mmd_vs_hpi(
+                ms.hpi_mmd,
+                config.mmd_pvalue_threshold,
+                marker=str(marker),
+                output_path=plots_dir / f"mmd_vs_hpi_{marker}.png",
+            )
 
     # Pass 2: GMM-label each marker using the FDR-significant conditions.
     frames: list[pd.DataFrame] = []
@@ -483,26 +664,40 @@ class vs time per condition, for every labeled marker.
         # Provenance row per tested (marker, condition): MMD gate + GMM summary.
         for cond, mmd in ms.cond_mmd.items():
             res = cond_gmm.get(cond)
-            mmd_rows.append(
-                {
-                    "marker": str(marker),
-                    "condition": str(cond),
-                    "n_perturbed": int((ms.perturbed_mask & (ms.conditions == cond)).sum()),
-                    "mmd2": mmd.mmd2,
-                    "p_raw": mmd.p_value,
-                    "p_adjusted": p_adj_by_key[(marker, cond)],
-                    "mmd_significant": cond in significant.get(marker, set()),
-                    "gmm_separated": bool(res.separated) if res is not None else False,
-                    "gmm_remodel_weight": float(res.gmm.weights_[res.remod_component]) if res is not None else np.nan,
-                    "n_confident_positive": int((res.hard_label == 1).sum()) if res is not None else 0,
-                    "gmm_bic": res.bic if res is not None else np.nan,
-                    "gmm_aic": res.aic if res is not None else np.nan,
-                    "gmm_bic_1comp": res.bic_1comp if res is not None else np.nan,
-                    "gmm_delta_bic": (res.bic_1comp - res.bic) if res is not None else np.nan,
-                }
-            )
-        # Witness-GMM diagnostic for every fitted condition (labeled or unimodal-skipped).
+            row = {
+                "marker": str(marker),
+                "condition": str(cond),
+                "n_perturbed": int((ms.perturbed_mask & (ms.conditions == cond)).sum()),
+                "mmd2": mmd.mmd2,
+                "p_raw": mmd.p_value,
+                "p_adjusted": p_adj_by_key[(marker, cond)],
+                "mmd_significant": cond in significant.get(marker, set()),
+                "gate": config.gate,
+                "n_confident_positive": int((res.hard_label == 1).sum()) if res is not None else 0,
+            }
+            if isinstance(res, ControlAnchoredResult):
+                row.update(
+                    {
+                        "remodel_fraction": 1.0 - res.pi_baseline,
+                        "control_fp": res.control_fp,
+                        "posterior_threshold": res.threshold,
+                    }
+                )
+            elif res is not None:
+                row.update(
+                    {
+                        "gmm_separated": bool(res.separated),
+                        "gmm_remodel_weight": float(res.gmm.weights_[res.remod_component]),
+                        "gmm_bic": res.bic,
+                        "gmm_aic": res.aic,
+                        "gmm_delta_bic": res.bic_1comp - res.bic,
+                    }
+                )
+            mmd_rows.append(row)
+        # Witness-GMM diagnostic (GMM gate only; control-anchored uses a different model).
         for cond, res in cond_gmm.items():
+            if isinstance(res, ControlAnchoredResult):
+                continue
             plot_witness_gmm(
                 result.cond_scores[cond],
                 result.control_scores,
@@ -532,28 +727,38 @@ class vs time per condition, for every labeled marker.
         )
         frames.append(frame)
 
-    if not frames:
-        raise RuntimeError("No markers produced labels — check references, thresholds, and condition_column.")
-
-    out = pd.concat(frames, ignore_index=True)
     # Disambiguate by marker: sibling single-marker configs often share a
     # label_column (e.g. three organelle markers → organelle_remodeling_state),
     # which would clobber a bare <label_column>.<fmt>. Prefix with the filtered
     # marker(s) so each config writes its own file.
     stem = f"{'_'.join(config.marker_filters)}_{config.label_column}" if config.marker_filters else config.label_column
-    output_path = labels_dir / f"{stem}.{config.annotation_format}"
     labels_dir.mkdir(parents=True, exist_ok=True)
+
+    # Population-level provenance sidecar always written — even when a marker
+    # abstains (no bimodal + significant condition), so the MMD/GMM evidence for
+    # the abstain decision is auditable.
+    mmd_path = labels_dir / f"{stem}_mmd.csv"
+    pd.DataFrame(mmd_rows).to_csv(mmd_path, index=False)
+    _logger.info("Wrote MMD/GMM provenance to %s", mmd_path)
+
+    if not frames:
+        # Legitimate abstain (e.g. time-matched witness leaves a weak channel
+        # unimodal): write no annotation, but do not crash — the diagnostics +
+        # sidecar above record why.
+        _logger.warning(
+            "No markers produced confident labels (no significant + bimodal condition). "
+            "Wrote diagnostics + %s but no annotation file.",
+            mmd_path.name,
+        )
+        return None
+
+    out = pd.concat(frames, ignore_index=True)
+    output_path = labels_dir / f"{stem}.{config.annotation_format}"
     if output_path.suffix == ".parquet":
         out.to_parquet(output_path, index=False)
     else:
         out.to_csv(output_path, index=False)
     _logger.info("Wrote %d annotations (%s) to %s", len(out), config.label_column, output_path)
-
-    # Population-level provenance sidecar: one row per (marker, condition) with the
-    # MMD gate (mmd2, raw/BY-adjusted p, significance) and the GMM summary.
-    mmd_path = labels_dir / f"{stem}_mmd.csv"
-    pd.DataFrame(mmd_rows).to_csv(mmd_path, index=False)
-    _logger.info("Wrote MMD/GMM provenance to %s", mmd_path)
     return output_path
 
 
@@ -574,7 +779,10 @@ def main(config_path: Path) -> None:
     raw = load_config(config_path)
     config = WitnessGmmLabelsConfig(**raw["witness_gmm_labels"])
     out = generate_witness_gmm_annotation(config)
-    click.echo(f"Wrote witness-GMM annotation to {out}")
+    if out is None:
+        click.echo("No confident labels produced (marker abstained); diagnostics written, no annotation file.")
+    else:
+        click.echo(f"Wrote witness-GMM annotation to {out}")
 
 
 if __name__ == "__main__":
diff --git a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels_test.py b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels_test.py
index 06abacd5b..38b5669e3 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels_test.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_labels_test.py
@@ -3,9 +3,11 @@
 import anndata as ad
 import numpy as np
 import pandas as pd
+import pytest
 
 from dynaclr.evaluation.evaluate_config import WitnessGmmExperiment, WitnessGmmLabelsConfig
 from dynaclr.evaluation.linear_classifiers.witness_gmm_labels import (
+    _hpi_bin_edges,
     _well_prefix_mask,
     compute_marker_scores,
     generate_witness_gmm_annotation,
@@ -75,7 +77,7 @@ def _make_separable_embeddings(
     return ad.AnnData(X=X, obs=obs, var=var)
 
 
-def _config(output_dir, experiment="exp_A", embeddings_zarr="unused.zarr", annotation_format="csv"):
+def _config(output_dir, experiment="exp_A", embeddings_zarr="unused.zarr", annotation_format="csv", **overrides):
     return WitnessGmmLabelsConfig(
         experiments=[
             WitnessGmmExperiment(
@@ -91,9 +93,26 @@ def _config(output_dir, experiment="exp_A", embeddings_zarr="unused.zarr", annot
         condition_column="perturbation",
         output_dir=str(output_dir),
         annotation_format=annotation_format,
+        **overrides,
     )
 
 
+def test_hpi_bin_edges_include_maximum_on_boundary():
+    hpi = np.asarray([0.0, 2.0, 4.0])
+
+    edges = _hpi_bin_edges(hpi, width=2.0)
+
+    assert edges.tolist() == [0.0, 2.0, 4.0, 6.0]
+    assert any(lo <= hpi.max() < hi for lo, hi in zip(edges[:-1], edges[1:]))
+
+
+@pytest.mark.parametrize("field", ["witness_time_bin_hours", "mmd_hpi_bin_hours"])
+@pytest.mark.parametrize("width", [0.0, -1.0, np.inf, np.nan])
+def test_hpi_bin_width_must_be_finite_and_positive(tmp_path, field, width):
+    with pytest.raises(ValueError, match=field):
+        _config(tmp_path, **{field: width})
+
+
 def test_well_prefix_mask_no_spurious_prefix_match():
     """'A/1' must not match 'A/10/...' — matching is on path components."""
     fov = pd.Series(["A/1/000000", "A/10/000000", "B/2/000000"])
diff --git a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_plots.py b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_plots.py
index b685433c0..b2319fe64 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_plots.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/linear_classifiers/witness_gmm_plots.py
@@ -13,6 +13,9 @@
 - :func:`plot_remodeling_vs_time` — per marker: the fraction of cells in the
   perturbed (remodeled/infected) class vs time, per condition — the biological
   kinetics the labels imply.
+- :func:`plot_mmd_vs_hpi` — per marker: MMD²(control, condition) per HPI bin —
+  the population-divergence kinetics (label-free), with a control-vs-control null
+  band. Answers "when, and how much, does the population diverge from control?"
 """
 
 from __future__ import annotations
@@ -322,3 +325,69 @@ def plot_remodeling_vs_time(
     ax.legend(frameon=True, fontsize=9)
     fig.savefig(output_path, dpi=150, bbox_inches="tight")
     plt.close(fig)
+
+
+def plot_mmd_vs_hpi(
+    hpi_mmd: dict,
+    pvalue_threshold: float,
+    marker: str,
+    output_path: Path,
+) -> None:
+    """Plot MMD²(control, condition) per HPI bin — population-divergence kinetics.
+
+    Each condition curve traces how far its cell cloud sits from the control
+    reference at each time window (label-free — no GMM, no per-cell labels). The
+    ``__control_null__`` series (control split against itself) is drawn as a grey
+    band: the no-difference floor the condition curves must exceed to be real.
+    A condition curve that **grows** with HPI is the signature of a genuine,
+    progressive perturbation; one that sits flat near the null (or as a constant
+    offset that never grows) is weak/ambiguous — possibly a fixed well effect
+    rather than a time-developing response. Points failing the per-bin
+    significance test (raw p > ``pvalue_threshold``) are drawn hollow.
+
+    Parameters
+    ----------
+    hpi_mmd : dict
+        ``{condition: [(hpi_center, mmd2, p), ...]}`` from ``_compute_hpi_mmd``;
+        key ``"__control_null__"`` is the control-vs-control baseline.
+    pvalue_threshold : float
+        Per-bin raw-p threshold for the hollow/filled marker distinction.
+    marker : str
+        Marker name (title).
+    output_path : Path
+        Output file path.
+    """
+    if not hpi_mmd:
+        return
+
+    fig, ax = plt.subplots(figsize=(9, 5))
+    colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
+
+    null = hpi_mmd.get("__control_null__")
+    if null:
+        arr = np.array(sorted(null))
+        ax.plot(arr[:, 0], arr[:, 1], color="0.5", ls="--", lw=1.2, label="control vs control (null)")
+        # Shade up to the max null MMD² as the no-difference band.
+        ax.axhspan(0, float(arr[:, 1].max()), color="0.85", alpha=0.5, zorder=0)
+
+    ci = 0
+    for cond, series in hpi_mmd.items():
+        if cond == "__control_null__":
+            continue
+        arr = np.array(sorted(series))
+        centers, mmd2, pvals = arr[:, 0], arr[:, 1], arr[:, 2]
+        color = colors[ci % len(colors)]
+        ci += 1
+        ax.plot(centers, mmd2, color=color, lw=1.8, label=f"control vs {cond}")
+        sig = pvals <= pvalue_threshold
+        ax.scatter(centers[sig], mmd2[sig], color=color, s=30, zorder=3)
+        ax.scatter(centers[~sig], mmd2[~sig], facecolors="none", edgecolors=color, s=30, zorder=3, label="_nolegend_")
+
+    ax.set_ylim(bottom=0)
+    ax.set_title(f"MMD² vs time — {marker}\n(population divergence from control; hollow = n.s.)", fontsize=11)
+    ax.set_xlabel("hours post perturbation", fontsize=11)
+    ax.set_ylabel("MMD²  (control vs condition)", fontsize=11)
+    ax.grid(True, alpha=0.3)
+    ax.legend(frameon=True, fontsize=9)
+    fig.savefig(output_path, dpi=150, bbox_inches="tight")
+    plt.close(fig)
diff --git a/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction.py b/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction.py
index d318ed610..9334037f8 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction.py
@@ -221,11 +221,10 @@ def apply_lot_correction(
     """Apply a fitted LOT pipeline to an embedding zarr.
 
     Transforms all cells through StandardScaler → (optional PCA) → LOT and writes
-    an AnnData zarr whose ``.X`` contains the corrected embeddings. ``.obs`` and
-    the input ``.uns`` are preserved (plus a ``uns["lot_correction"]`` provenance
-    entry). ``obsm`` (e.g. ``X_backbone``, ``X_umap``, ``X_phate``, ``X_pca``),
-    ``varm``, ``obsp``, and ``layers`` are intentionally dropped: they were
-    computed in the *uncorrected* space and would contradict the corrected ``.X``.
+    an AnnData zarr whose ``.X`` contains the corrected embeddings. The matching
+    pre-LOT coordinates are stored in ``obsm["X_pre_lot"]`` for correction QC.
+    Input metadata are preserved; other arrays derived from the uncorrected space
+    are dropped.
 
     Parameters
     ----------
@@ -237,14 +236,17 @@ def apply_lot_correction(
         Path to write the corrected AnnData zarr.
     overwrite : bool, optional
         If ``False`` (default) and *output_zarr* already exists, raise.
+        Existing output is replaced only after the new store is fully written.
     """
     import shutil
+    import tempfile
 
+    input_zarr = Path(input_zarr)
     output_zarr = Path(output_zarr)
-    if output_zarr.exists():
-        if not overwrite:
-            raise FileExistsError(f"Output path already exists: {output_zarr}. Set overwrite=true to overwrite.")
-        shutil.rmtree(output_zarr)
+    if input_zarr.resolve() == output_zarr.resolve():
+        raise ValueError("input_zarr and output_zarr must be different paths.")
+    if output_zarr.exists() and not overwrite:
+        raise FileExistsError(f"Output path already exists: {output_zarr}. Set overwrite=true to overwrite.")
 
     _logger.info("Loading input zarr: %s", input_zarr)
     adata_in = ad.read_zarr(input_zarr)
@@ -275,6 +277,7 @@ def apply_lot_correction(
         pass
 
     adata_out = ad.AnnData(X=Z_corrected.astype(np.float32), obs=obs, uns=dict(adata_in.uns))
+    adata_out.obsm["X_pre_lot"] = np.asarray(Z, dtype=np.float32)
     adata_out.var.index = adata_out.var.index.astype(object)
     adata_out.uns["lot_correction"] = {
         "source_zarr": str(input_zarr),
@@ -283,9 +286,31 @@ def apply_lot_correction(
         "pca_variance_explained": pipeline.get("pca_variance_explained"),
     }
 
-    _logger.info("Writing corrected zarr: %s", output_zarr)
-    adata_out.write_zarr(output_zarr, convert_strings_to_categoricals=False)
-    _logger.info("Done.")
+    output_zarr.parent.mkdir(parents=True, exist_ok=True)
+    temp_root = Path(tempfile.mkdtemp(prefix=f".{output_zarr.name}.", dir=output_zarr.parent))
+    temp_output = temp_root / "new.zarr"
+    backup_output = temp_root / "previous.zarr"
+    keep_temp = False
+    try:
+        _logger.info("Writing corrected zarr: %s", temp_output)
+        adata_out.write_zarr(temp_output, convert_strings_to_categoricals=False)
+        if output_zarr.exists():
+            output_zarr.rename(backup_output)
+        try:
+            temp_output.rename(output_zarr)
+        except Exception:
+            if backup_output.exists():
+                try:
+                    backup_output.rename(output_zarr)
+                except Exception:
+                    keep_temp = True
+                    _logger.exception("Could not restore previous output; backup retained at %s", backup_output)
+            raise
+    finally:
+        if not keep_temp:
+            shutil.rmtree(temp_root, ignore_errors=True)
+
+    _logger.info("Done: %s", output_zarr)
 
 
 def save_lot_pipeline(pipeline: dict, path: Union[str, Path]) -> None:
diff --git a/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction_test.py b/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction_test.py
index 28ca92307..7a363e471 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction_test.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/lot_correction/lot_correction_test.py
@@ -30,6 +30,12 @@ def _make_adata(n: int, d: int, seed: int, shift: float = 0.0) -> ad.AnnData:
     return ad.AnnData(X=X, obs=obs)
 
 
+def _write_zarr(adata: ad.AnnData, path: Path) -> None:
+    ad.settings.allow_write_nullable_strings = True
+    adata.var.index = adata.var.index.astype(object)
+    adata.write_zarr(path, convert_strings_to_categoricals=False)
+
+
 def test_pool_embeddings_concatenates_rows():
     a = _make_adata(10, 8, seed=0)
     b = _make_adata(15, 8, seed=1)
@@ -118,9 +124,65 @@ def test_fit_save_load_then_apply_to_many(tmp_path, n_pca):
         apply_lot_correction(input_zarr, loaded, output_zarr)
         out = ad.read_zarr(output_zarr)
         assert out.shape == (20, expected_dim)
+        assert out.obsm["X_pre_lot"].shape == (20, expected_dim)
+        expected_pre = loaded["scaler"].transform(input_adata.X)
+        if loaded["pca"] is not None:
+            expected_pre = loaded["pca"].transform(expected_pre)
+        np.testing.assert_allclose(out.obsm["X_pre_lot"], expected_pre, rtol=1e-5)
         assert out.uns["lot_correction"]["channel"] == "Phase3D"
 
 
+def test_apply_rejects_same_input_and_output(tmp_path):
+    source = [_make_adata(40, 16, seed=1)]
+    target = [_make_adata(50, 16, seed=2)]
+    pipeline = fit_lot_correction(source, target, n_pca=5, ns_lot=50, random_seed=0)
+    input_zarr = tmp_path / "input.zarr"
+    _write_zarr(source[0], input_zarr)
+
+    with pytest.raises(ValueError, match="must be different"):
+        apply_lot_correction(input_zarr, pipeline, input_zarr, overwrite=True)
+
+    assert ad.read_zarr(input_zarr).shape == source[0].shape
+
+
+def test_apply_preserves_existing_output_on_failure(tmp_path):
+    source = [_make_adata(40, 16, seed=1)]
+    target = [_make_adata(50, 16, seed=2)]
+    pipeline = fit_lot_correction(source, target, n_pca=5, ns_lot=50, random_seed=0)
+
+    bad_input = _make_adata(20, 15, seed=3)
+    input_zarr = tmp_path / "bad_input.zarr"
+    _write_zarr(bad_input, input_zarr)
+
+    previous = _make_adata(7, 3, seed=4)
+    output_zarr = tmp_path / "output.zarr"
+    _write_zarr(previous, output_zarr)
+    previous_x = previous.X.copy()
+
+    with pytest.raises(ValueError):
+        apply_lot_correction(input_zarr, pipeline, output_zarr, overwrite=True)
+
+    preserved = ad.read_zarr(output_zarr)
+    assert preserved.shape == previous.shape
+    np.testing.assert_array_equal(preserved.X, previous_x)
+
+
+def test_apply_replaces_existing_output_after_success(tmp_path):
+    source = [_make_adata(40, 16, seed=1)]
+    target = [_make_adata(50, 16, seed=2)]
+    pipeline = fit_lot_correction(source, target, n_pca=5, ns_lot=50, random_seed=0)
+
+    input_zarr = tmp_path / "input.zarr"
+    _write_zarr(source[0], input_zarr)
+    output_zarr = tmp_path / "output.zarr"
+    _write_zarr(_make_adata(7, 3, seed=4), output_zarr)
+
+    apply_lot_correction(input_zarr, pipeline, output_zarr, overwrite=True)
+
+    assert ad.read_zarr(output_zarr).shape == (40, 5)
+    assert not list(tmp_path.glob(f".{output_zarr.name}.*"))
+
+
 def test_coerce_obs_for_zarr_handles_categorical_string():
     """Categorical-over-string obs must coerce to a zarr-writable object dtype.
 
diff --git a/applications/dynaclr/src/dynaclr/evaluation/mmd/compute_mmd.py b/applications/dynaclr/src/dynaclr/evaluation/mmd/compute_mmd.py
index 5c049bd8b..cd8772cc2 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/mmd/compute_mmd.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/mmd/compute_mmd.py
@@ -46,6 +46,37 @@ def _extract_embeddings(adata: ad.AnnData, embedding_key: str | None) -> np.ndar
     return np.asarray(X)
 
 
+def _load_adatas_by_experiment(input_paths: list[str]) -> dict[str, ad.AnnData]:
+    """Load and concatenate every store belonging to the same experiment.
+
+    Prediction writes one store per experiment and marker. Grouping only by
+    experiment with a dict comprehension would silently retain the final marker.
+    """
+    grouped: dict[str, list[ad.AnnData]] = {}
+    for input_path in input_paths:
+        adata = ad.read_zarr(input_path)
+        if "experiment" not in adata.obs.columns:
+            raise KeyError(f"obs column 'experiment' not found in {input_path}.")
+        experiments = adata.obs["experiment"].dropna().unique()
+        if len(experiments) != 1:
+            raise ValueError(f"Expected exactly one experiment in {input_path}, found {list(experiments)}.")
+        grouped.setdefault(str(experiments[0]), []).append(adata)
+
+    return {
+        experiment: (parts[0] if len(parts) == 1 else ad.concat(parts, join="inner", merge="same", index_unique="-"))
+        for experiment, parts in grouped.items()
+    }
+
+
+def _experiment_marker_pairs(input_paths: list[str]) -> set[tuple[str, str]]:
+    pairs = set()
+    for experiment, adata in _load_adatas_by_experiment(input_paths).items():
+        if "marker" not in adata.obs:
+            raise KeyError(f"obs column 'marker' not found for {experiment}.")
+        pairs.update((experiment, str(marker)) for marker in adata.obs["marker"].unique())
+    return pairs
+
+
 def _run_one_comparison(
     emb_a: np.ndarray,
     emb_b: np.ndarray,
@@ -414,7 +445,7 @@ def run_mmd_combined(config: MMDCombinedConfig) -> pd.DataFrame:
     """
     from itertools import combinations
 
-    adatas = {ad.read_zarr(p).obs["experiment"].iloc[0]: ad.read_zarr(p) for p in config.input_paths}
+    adatas = _load_adatas_by_experiment(config.input_paths)
 
     if config.obs_filter:
         filtered = {}
@@ -521,32 +552,48 @@ def run_mmd_combined(config: MMDCombinedConfig) -> pd.DataFrame:
 
 
 def run_mmd_over_time(config: MMDOverTimeConfig) -> pd.DataFrame:
-    """Run combined cross-experiment MMD on pre- and post-correction embeddings.
-
-    Runs :func:`run_mmd_combined` twice — once on ``input_paths`` (uncorrected,
-    ``correction="pre"``) and once on ``corrected_paths`` (LOT-corrected,
-    ``correction="post"``) — using identical settings, then concatenates the
-    results with a ``correction`` column so the batch effect before and after
-    correction can be compared over time in a single output.
-
-    Parameters
-    ----------
-    config : MMDOverTimeConfig
-        Over-time analysis configuration (combined config + ``corrected_paths``).
+    """Compare pre- and post-LOT MMD in the same scaler/PCA space.
 
-    Returns
-    -------
-    pd.DataFrame
-        Same columns as :func:`run_mmd_combined` plus a ``correction`` column
-        (``"pre"`` / ``"post"``) and a ``pair_kind`` column (``"cross"`` for
-        source↔target pairs, ``"within"`` for source↔source pairs).
+    Corrected stores keep pre-LOT coordinates in ``obsm["X_pre_lot"]`` and
+    post-LOT coordinates in ``.X``.
     """
-    base = config.model_dump(exclude={"corrected_paths", "target_experiments"})
-
-    pre = run_mmd_combined(MMDCombinedConfig(**base))
+    if config.embedding_key is not None:
+        raise ValueError("embedding_key must be unset for over-time LOT MMD.")
+
+    raw_pairs = _experiment_marker_pairs(config.input_paths)
+    corrected_pairs = _experiment_marker_pairs(config.corrected_paths)
+    if raw_pairs != corrected_pairs:
+        raise ValueError("input_paths and corrected_paths contain different experiment-marker populations.")
+
+    pre_lot_key = "X_pre_lot"
+    for path in config.corrected_paths:
+        if pre_lot_key not in ad.read_zarr(path).obsm:
+            raise ValueError(f"{path} has no obsm['{pre_lot_key}']; re-run LOT correction first.")
+
+    base = config.model_dump(
+        exclude={
+            "input_paths",
+            "corrected_paths",
+            "target_experiments",
+            "embedding_key",
+        }
+    )
+    pre = run_mmd_combined(
+        MMDCombinedConfig(
+            **base,
+            input_paths=config.corrected_paths,
+            embedding_key=pre_lot_key,
+        )
+    )
     pre["correction"] = "pre"
 
-    post = run_mmd_combined(MMDCombinedConfig(**{**base, "input_paths": config.corrected_paths}))
+    post = run_mmd_combined(
+        MMDCombinedConfig(
+            **base,
+            input_paths=config.corrected_paths,
+            embedding_key=None,
+        )
+    )
     post["correction"] = "post"
 
     df = pd.concat([pre, post], ignore_index=True)
diff --git a/applications/dynaclr/src/dynaclr/evaluation/mmd/config.py b/applications/dynaclr/src/dynaclr/evaluation/mmd/config.py
index b2f439a26..60bb46798 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/mmd/config.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/mmd/config.py
@@ -202,12 +202,10 @@ class MMDCombinedConfig(_MMDBaseConfig):
 class MMDOverTimeConfig(MMDCombinedConfig):
     """Pre/post batch-effect MMD over time in a single run.
 
-    Runs the combined cross-experiment MMD (per marker × condition × time bin)
-    on both the uncorrected embeddings (``input_paths``) and their LOT-corrected
-    counterparts (``corrected_paths``), tags each result set with a
-    ``correction`` column (``"pre"`` / ``"post"``), and returns one combined
-    DataFrame — so the batch effect before and after correction can be plotted
-    as two series over time instead of living in two separate output folders.
+    Runs combined cross-experiment MMD on the pre-LOT coordinates stored in
+    ``corrected_paths[*].obsm["X_pre_lot"]`` and post-LOT coordinates in
+    ``corrected_paths[*].X``. Both therefore use the same scaler/PCA space.
+    ``input_paths`` identify the expected experiment-marker populations.
 
     Parameters
     ----------
diff --git a/applications/dynaclr/src/dynaclr/evaluation/orchestration/matrix.py b/applications/dynaclr/src/dynaclr/evaluation/orchestration/matrix.py
index 41fe19aa7..8cfcb2f41 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/orchestration/matrix.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/orchestration/matrix.py
@@ -197,10 +197,15 @@ def build_train_cmd(model: dict) -> list[str]:
 def build_predict_cmd(model: dict, ckpt_name: str, checkpoint: str) -> list[str]:
     """The predict sbatch command for one checkpoint of a model.
 
-    ``markers`` (optional list in the matrix row) is passed as a comma-joined 8th
-    positional; empty = all channels.
+    ``markers`` and ``predict_flags`` are forwarded as positional values to
+    the prediction wrapper.
     """
     markers = ",".join(model["markers"]) if model.get("markers") else ""
+    predict_flags = model.get("predict_flags", {})
+    z_range = predict_flags.get("z_range")
+    if z_range is not None and len(z_range) != 2:
+        raise ValueError("predict_flags.z_range must contain exactly two values.")
+    z_start, z_end = ("", "") if z_range is None else map(str, z_range)
     return [
         "sbatch",
         str(_PREDICT_SBATCH),
@@ -211,12 +216,25 @@ def build_predict_cmd(model: dict, ckpt_name: str, checkpoint: str) -> list[str]
         str(model["datasets_root"]),
         checkpoint,  # empty = derive from run dir
         markers,  # empty = all channels
+        z_start,
+        z_end,
+        str(predict_flags.get("z_reduction") or ""),
+        str(predict_flags.get("reference_pixel_size") or ""),
+        str(predict_flags.get("batch_size") or ""),
     ]
 
 
 def build_eval_cmd(model: dict, ckpt_name: str) -> list[str]:
     """The eval sbatch command for one checkpoint of a model."""
-    return ["sbatch", str(_EVAL_SBATCH), model["eval_config"], model["family"], model["run"], ckpt_name]
+    return [
+        "sbatch",
+        str(_EVAL_SBATCH),
+        model["eval_config"],
+        model["family"],
+        model["run"],
+        ckpt_name,
+        str(model["datasets_root"]),
+    ]
 
 
 def build_stage_cmds(model: dict, stages: tuple[str, ...]) -> list[tuple[str, list[str]]]:
diff --git a/applications/dynaclr/src/dynaclr/evaluation/orchestration/predict_batch.py b/applications/dynaclr/src/dynaclr/evaluation/orchestration/predict_batch.py
index 9b0150f3e..bbd79b8f3 100644
--- a/applications/dynaclr/src/dynaclr/evaluation/orchestration/predict_batch.py
+++ b/applications/dynaclr/src/dynaclr/evaluation/orchestration/predict_batch.py
@@ -253,6 +253,20 @@ def build_predict_cmd(
     show_default=True,
     help="Base under which datasets live.",
 )
+@click.option("--z-range", nargs=2, type=int, default=None, help="Z window forwarded to predict-triplet.")
+@click.option(
+    "--z-reduction",
+    type=click.Choice(["mip", "center"]),
+    default=None,
+    help="Z reduction forwarded to predict-triplet.",
+)
+@click.option(
+    "--reference-pixel-size",
+    type=float,
+    default=None,
+    help="Reference pixel size forwarded to predict-triplet.",
+)
+@click.option("--batch-size", type=int, default=None, help="Batch size forwarded to predict-triplet.")
 @click.option("--markers", multiple=True, default=None, help="Marker subset (default: all channels).")
 @click.option("--no-labelfree", is_flag=True, help="Skip label-free (phase/brightfield) channels.")
 @click.option("--num-workers", type=int, default=0, help="Predict dataloader workers (must be 0).")
@@ -281,6 +295,10 @@ def main(
     checkpoint: Path | None,
     models_root: Path,
     datasets_root: Path,
+    z_range: tuple[int, int] | None,
+    z_reduction: str | None,
+    reference_pixel_size: float | None,
+    batch_size: int | None,
     markers: tuple[str, ...],
     no_labelfree: bool,
     num_workers: int,
@@ -298,6 +316,16 @@ def main(
         preflight(collection, workspace_dir, auto_normalize=auto_normalize)
 
     resolved_checkpoint = checkpoint or checkpoint_path(model_family, run, ckpt_name, models_root=models_root)
+    predict_flags = {}
+    if z_range is not None:
+        predict_flags["z_range"] = list(z_range)
+    if z_reduction is not None:
+        predict_flags["z_reduction"] = z_reduction
+    if reference_pixel_size is not None:
+        predict_flags["reference_pixel_size"] = reference_pixel_size
+    if batch_size is not None:
+        predict_flags["batch_size"] = batch_size
+
     cmd = build_predict_cmd(
         collection,
         resolved_checkpoint,
@@ -305,6 +333,7 @@ def main(
         run,
         ckpt_name,
         datasets_root,
+        predict_flags=predict_flags or None,
         markers=list(markers) if markers else None,
         no_labelfree=no_labelfree,
         num_workers=num_workers,
diff --git a/applications/dynaclr/tests/test_datamodule.py b/applications/dynaclr/tests/test_datamodule.py
index 25cac7e52..9c08981a0 100644
--- a/applications/dynaclr/tests/test_datamodule.py
+++ b/applications/dynaclr/tests/test_datamodule.py
@@ -5,6 +5,7 @@
 from __future__ import annotations
 
 import pytest
+import torch
 
 from viscy_data.cell_index import build_timelapse_cell_index
 
@@ -147,6 +148,43 @@ def test_init_exposes_all_hyperparameters(self, two_experiments):
         assert dm.cache_pool_bytes == 1024
         assert dm.seed == 42
 
+    def test_rejects_unknown_sequence_augmentation_mode(self, two_experiments):
+        from dynaclr.data.datamodule import MultiExperimentDataModule
+
+        parquet_path, _ = two_experiments
+        with pytest.raises(ValueError, match="sequence_augment"):
+            MultiExperimentDataModule(
+                cell_index_path=str(parquet_path),
+                z_window=1,
+                sequence_augment="per_frame",
+            )
+
+
+class TestSequenceAugmentation:
+    def test_consistent_mode_reuses_transform_per_track(self, two_experiments):
+        from monai.transforms import Compose
+
+        from dynaclr.data.datamodule import MultiExperimentDataModule
+        from viscy_transforms import BatchedRandFlipd
+
+        parquet_path, _ = two_experiments
+        dm = MultiExperimentDataModule(
+            cell_index_path=str(parquet_path),
+            z_window=1,
+            sequence_length=3,
+            sequence_augment="consistent",
+        )
+        dm._channel_names = ["Phase"]
+        dm._sequence_transform = Compose([BatchedRandFlipd(keys=["Phase"], spatial_axes=[2], prob=0.5)])
+
+        tracks = torch.arange(8 * 15, dtype=torch.float32).reshape(8, 1, 1, 3, 5)
+        sequence = tracks.repeat_interleave(dm.sequence_length, dim=0)
+        transformed = dm._transform_sequence_consistently(sequence, None, None)
+        transformed = transformed.reshape(8, dm.sequence_length, 1, 1, 3, 5)
+
+        assert torch.equal(transformed[:, 0], transformed[:, 1])
+        assert torch.equal(transformed[:, 1], transformed[:, 2])
+
 
 class TestTrainValSplitByExperiment:
     """DATA-04: Train/val split is by whole experiments, not individual FOVs."""
diff --git a/applications/dynaclr/tests/test_dataset.py b/applications/dynaclr/tests/test_dataset.py
index ab058d369..3da21439a 100644
--- a/applications/dynaclr/tests/test_dataset.py
+++ b/applications/dynaclr/tests/test_dataset.py
@@ -538,6 +538,72 @@ def test_self_positive_pixel_values_identical(self, single_experiment_index):
             "Self-positive: anchor and positive tensors must be identical before augmentation"
         )
 
+    def test_self_positive_can_emit_sequence(self, single_experiment_index):
+        """Self positives still build the temporal lookup used by sequences."""
+        from dynaclr.data.dataset import MultiExperimentTripletDataset
+
+        ds = MultiExperimentTripletDataset(
+            index=single_experiment_index,
+            fit=True,
+            positive_cell_source="self",
+            channels_per_sample=1,
+            emit_sequence=True,
+            sequence_length=3,
+        )
+
+        assert ds._lineage_timepoints
+        batch = ds.__getitems__([0])
+        assert batch["sequence"].shape[0] == 3
+        assert batch["sequence_valid"].shape == (1,)
+
+
+class TestSequenceSampling:
+    """Temporal stencils must preserve exact track identity."""
+
+    @staticmethod
+    def _sampler_dataset(track_ids, timepoints):
+        from dynaclr.data.dataset import MultiExperimentTripletDataset
+
+        ds = object.__new__(MultiExperimentTripletDataset)
+        ds._va_arrays = {
+            "experiment": np.asarray(["exp"]),
+            "lineage_id": np.asarray(["lineage"]),
+            "t": np.asarray([0]),
+            "marker": np.asarray(["GFP"]),
+            "global_track_id": np.asarray([track_ids[0]]),
+        }
+        ds._tr_arrays = {
+            "marker": np.asarray(["GFP"] * len(track_ids)),
+            "global_track_id": np.asarray(track_ids),
+        }
+        ds._lineage_timepoints = {("exp", "lineage"): timepoints}
+        ds.sequence_length = 3
+        ds.sequence_tau_frames = 1
+        ds._rng = np.random.default_rng(0)
+        return ds
+
+    def test_selects_same_track_when_siblings_share_lineage(self):
+        ds = self._sampler_dataset(
+            ["track-a", "track-b", "track-a", "track-b", "track-a", "track-b"],
+            {0: [0, 1], 1: [2, 3], 2: [4, 5]},
+        )
+
+        rows, valid = ds._sample_sequence_indices(np.asarray([0]))
+
+        assert valid.tolist() == [True]
+        assert rows.tolist() == [0, 2, 4]
+
+    def test_invalidates_stencil_at_division_boundary(self):
+        ds = self._sampler_dataset(
+            ["parent", "daughter-a", "daughter-b", "daughter-a", "daughter-b"],
+            {0: [0], 1: [1, 2], 2: [3, 4]},
+        )
+
+        rows, valid = ds._sample_sequence_indices(np.asarray([0]))
+
+        assert valid.tolist() == [False]
+        assert rows.tolist() == [0, 0, 0]
+
 
 class TestTimepointStatisticsResolution:
     """Verify that timepoint_statistics norm_meta resolves the correct timepoint."""
diff --git a/applications/dynaclr/tests/test_engine.py b/applications/dynaclr/tests/test_engine.py
index 5c04fb4d0..555fb10ab 100644
--- a/applications/dynaclr/tests/test_engine.py
+++ b/applications/dynaclr/tests/test_engine.py
@@ -6,6 +6,7 @@
 
 from dynaclr.engine import ContrastiveModule
 from viscy_models.components.heads import ClassificationHead
+from viscy_models.contrastive.predictor import Predictor
 
 
 def test_contrastive_module_init(_SimpleEncoder, synth_dims):
@@ -48,6 +49,42 @@ def test_contrastive_module_forward(_SimpleEncoder, synth_dims):
     assert projections.shape == (2, 32)
 
 
+def test_all_invalid_sequences_keep_predictor_in_backward_graph(_SimpleEncoder, synth_dims):
+    """Empty masked inputs produce zero, non-None predictor gradients for DDP."""
+    predictor = Predictor(dim=64)
+    module = ContrastiveModule(
+        encoder=_SimpleEncoder(),
+        predictor=predictor,
+        lambda_pred=1.0,
+        example_input_array_shape=(
+            1,
+            synth_dims["c"],
+            synth_dims["d"],
+            synth_dims["h"],
+            synth_dims["w"],
+        ),
+    )
+    module.log = lambda *args, **kwargs: None
+    batch = {
+        "sequence": torch.randn(
+            6,
+            synth_dims["c"],
+            synth_dims["d"],
+            synth_dims["h"],
+            synth_dims["w"],
+        ),
+        "sequence_length": 3,
+        "sequence_valid": torch.zeros(2, dtype=torch.bool),
+    }
+
+    loss = module._temporal_losses(batch, "train")
+    loss.backward()
+
+    assert loss.item() == 0.0
+    assert all(parameter.grad is not None for parameter in predictor.parameters())
+    assert all(torch.count_nonzero(parameter.grad) == 0 for parameter in predictor.parameters())
+
+
 def test_embedding_pca_logged_every_n_epochs(_SimpleEncoder, _SyntheticTripletDataModule, synth_dims):
     """PCA logging is triggered at epochs 0, n, 2n, ... and not in between."""
     seed_everything(0)
diff --git a/applications/dynaclr/tests/test_mmd.py b/applications/dynaclr/tests/test_mmd.py
index 08637a8d8..e3831d53b 100644
--- a/applications/dynaclr/tests/test_mmd.py
+++ b/applications/dynaclr/tests/test_mmd.py
@@ -7,8 +7,20 @@
 import pandas as pd
 import pytest
 
-from dynaclr.evaluation.mmd.compute_mmd import run_mmd_analysis, run_mmd_pooled
-from dynaclr.evaluation.mmd.config import ComparisonSpec, MMDEvalConfig, MMDPooledConfig, MMDSettings
+from dynaclr.evaluation.mmd.compute_mmd import (
+    run_mmd_analysis,
+    run_mmd_combined,
+    run_mmd_over_time,
+    run_mmd_pooled,
+)
+from dynaclr.evaluation.mmd.config import (
+    ComparisonSpec,
+    MMDCombinedConfig,
+    MMDEvalConfig,
+    MMDOverTimeConfig,
+    MMDPooledConfig,
+    MMDSettings,
+)
 from viscy_utils.evaluation.mmd import (
     compute_mmd_unbiased,
     median_heuristic,
@@ -467,6 +479,80 @@ def _save_adata_zarr(adata: ad.AnnData, path: str) -> None:
     adata.write_zarr(path)
 
 
+def test_run_mmd_combined_keeps_all_marker_stores(tmp_path):
+    """Separate marker stores for one experiment must not overwrite each other."""
+    paths = []
+    for exp_index, experiment in enumerate(["exp_a", "exp_b"]):
+        for marker_index, marker in enumerate(["SEC61B", "TOMM20"]):
+            rng = np.random.default_rng(exp_index * 10 + marker_index)
+            obs = pd.DataFrame(
+                {
+                    "experiment": [experiment] * 30,
+                    "marker": [marker] * 30,
+                    "perturbation": ["uninfected"] * 30,
+                }
+            )
+            adata = ad.AnnData(X=rng.normal(size=(30, 8)).astype(np.float32), obs=obs)
+            store = str(tmp_path / f"{experiment}_{marker}.zarr")
+            _save_adata_zarr(adata, store)
+            paths.append(store)
+
+    result = run_mmd_combined(
+        MMDCombinedConfig(
+            input_paths=paths,
+            output_dir=str(tmp_path / "out"),
+            mmd=MMDSettings(n_permutations=10, min_cells=10),
+        )
+    )
+
+    assert set(result["marker"]) == {"SEC61B", "TOMM20"}
+    assert len(result) == 2
+
+
+def test_run_mmd_over_time_uses_pre_lot_coordinates(tmp_path):
+    """Pre/post results must use matching scaler/PCA coordinates."""
+    raw_paths = []
+    corrected_paths = []
+    for index, experiment in enumerate(["exp_a", "exp_b"]):
+        rng = np.random.default_rng(index)
+        obs = pd.DataFrame(
+            {
+                "experiment": [experiment] * 30,
+                "marker": ["SEC61B"] * 30,
+                "perturbation": ["uninfected"] * 30,
+            }
+        )
+
+        raw_path = str(tmp_path / f"{experiment}_raw.zarr")
+        _save_adata_zarr(
+            ad.AnnData(X=rng.normal(size=(30, 7)).astype(np.float32), obs=obs.copy()),
+            raw_path,
+        )
+        raw_paths.append(raw_path)
+
+        corrected = ad.AnnData(
+            X=rng.normal(size=(30, 3)).astype(np.float32),
+            obs=obs.copy(),
+        )
+        corrected.obsm["X_pre_lot"] = rng.normal(size=(30, 3)).astype(np.float32)
+        corrected_path = str(tmp_path / f"{experiment}_corrected.zarr")
+        _save_adata_zarr(corrected, corrected_path)
+        corrected_paths.append(corrected_path)
+
+    result = run_mmd_over_time(
+        MMDOverTimeConfig(
+            input_paths=raw_paths,
+            corrected_paths=corrected_paths,
+            output_dir=str(tmp_path / "out"),
+            mmd=MMDSettings(n_permutations=10, min_cells=10),
+        )
+    )
+
+    assert set(result["correction"]) == {"pre", "post"}
+    assert set(result.loc[result["correction"] == "pre", "embedding_key"]) == {"X_pre_lot"}
+    assert set(result.loc[result["correction"] == "post", "embedding_key"]) == {"X"}
+
+
 def test_run_mmd_pooled_columns(tmp_path):
     """run_mmd_pooled returns expected columns including activity_zscore and q_value."""
     adata1 = _make_adata(n_cells=200, seed=0)
diff --git a/applications/dynaclr/tests/test_submit_matrix.py b/applications/dynaclr/tests/test_submit_matrix.py
index 3b6c8ac59..ebefe9774 100644
--- a/applications/dynaclr/tests/test_submit_matrix.py
+++ b/applications/dynaclr/tests/test_submit_matrix.py
@@ -1,5 +1,7 @@
 """Unit tests for submit_matrix: .sh parsing, defaults merge, chained commands."""
 
+import os
+import subprocess
 from pathlib import Path
 
 import pytest
@@ -114,9 +116,8 @@ def test_predict_empty_checkpoint_when_derived(tmp_path):
         {"ckpt_name": "last", "collection": "c.yml", "eval_config": "e.yaml", "datasets_root": "/d"},
     )
     predict_cmd = dict(submit_matrix.build_stage_cmds(model, ("predict",)))["predict"]
-    # positionals end with: … checkpoint(""), markers("")
-    assert predict_cmd[-2] == ""  # empty checkpoint positional
-    assert predict_cmd[-1] == ""  # empty markers positional
+    assert predict_cmd[7] == ""  # empty checkpoint positional
+    assert predict_cmd[8] == ""  # empty markers positional
 
 
 def test_predict_passes_markers(tmp_path):
@@ -127,7 +128,103 @@ def test_predict_passes_markers(tmp_path):
         {"ckpt_name": "last", "collection": "c.yml", "eval_config": "e.yaml", "datasets_root": "/d"},
     )
     predict_cmd = dict(submit_matrix.build_stage_cmds(model, ("predict",)))["predict"]
-    assert predict_cmd[-1] == "SEC61B,Phase3D"
+    assert predict_cmd[8] == "SEC61B,Phase3D"
+
+
+def test_matrix_forwards_predict_flags_and_eval_root():
+    model = {
+        "collection": "c.yml",
+        "family": "family",
+        "run": "run",
+        "datasets_root": "/datasets",
+        "eval_config": "e.yml",
+        "predict_flags": {
+            "z_range": [15, 45],
+            "z_reduction": "mip",
+            "reference_pixel_size": 0.1494,
+            "batch_size": 64,
+        },
+    }
+
+    predict_cmd = submit_matrix.build_predict_cmd(model, "last", "")
+    assert predict_cmd[-5:] == ["15", "45", "mip", "0.1494", "64"]
+
+    eval_cmd = submit_matrix.build_eval_cmd(model, "last")
+    assert eval_cmd[-1] == "/datasets"
+
+
+def test_predict_sbatch_repeats_marker_option(tmp_path):
+    """Each marker must have its own Click ``--markers`` option."""
+    fake_srun = tmp_path / "srun"
+    fake_srun.write_text('#!/bin/bash\nprintf "%s\\n" "$@"\n')
+    fake_srun.chmod(0o755)
+
+    workspace = Path(__file__).parents[3]
+    env = os.environ.copy()
+    env["PATH"] = f"{tmp_path}{os.pathsep}{env['PATH']}"
+    env["WORKSPACE_DIR"] = str(workspace)
+    result = subprocess.run(
+        [
+            "bash",
+            str(workspace / "applications/dynaclr/tools/predict.sbatch"),
+            "collection.yml",
+            "family",
+            "run",
+            "last",
+            "/datasets",
+            "",
+            "SEC61B,TOMM20",
+            "15",
+            "45",
+            "mip",
+            "0.1494",
+            "64",
+        ],
+        check=True,
+        capture_output=True,
+        text=True,
+        env=env,
+    )
+
+    args = result.stdout.splitlines()
+    marker_positions = [i for i, arg in enumerate(args) if arg == "--markers"]
+    assert [args[i + 1] for i in marker_positions] == ["SEC61B", "TOMM20"]
+    z_range_index = args.index("--z-range")
+    assert args[z_range_index + 1 : z_range_index + 3] == ["15", "45"]
+    assert args[args.index("--z-reduction") + 1] == "mip"
+    assert args[args.index("--reference-pixel-size") + 1] == "0.1494"
+    assert args[args.index("--batch-size") + 1] == "64"
+
+
+def test_eval_sbatch_forwards_datasets_root(tmp_path):
+    fake_module = tmp_path / "module"
+    fake_module.write_text("#!/bin/bash\nexit 0\n")
+    fake_module.chmod(0o755)
+    fake_uv = tmp_path / "uv"
+    fake_uv.write_text('#!/bin/bash\nprintf "%s\\n" "$@"\n')
+    fake_uv.chmod(0o755)
+
+    workspace = Path(__file__).parents[3]
+    env = os.environ.copy()
+    env["PATH"] = f"{tmp_path}{os.pathsep}{env['PATH']}"
+    env["WORKSPACE_DIR"] = str(workspace)
+    model = {
+        "eval_config": "e.yml",
+        "family": "family",
+        "run": "run",
+        "datasets_root": "/datasets",
+    }
+    cmd = submit_matrix.build_eval_cmd(model, "last")
+    result = subprocess.run(
+        ["bash", *cmd[1:]],
+        check=True,
+        capture_output=True,
+        text=True,
+        env=env,
+    )
+
+    args = result.stdout.splitlines()
+    assert args[args.index("--datasets-root") + 1] == "/datasets"
 
 
 # --- matrix-level preflight -------------------------------------------------
diff --git a/applications/dynaclr/tests/test_submit_predict.py b/applications/dynaclr/tests/test_submit_predict.py
index 6c604ed69..1b6cc804d 100644
--- a/applications/dynaclr/tests/test_submit_predict.py
+++ b/applications/dynaclr/tests/test_submit_predict.py
@@ -2,6 +2,8 @@
 
 from pathlib import Path
 
+from click.testing import CliRunner
+
 from dynaclr.evaluation.orchestration import predict_batch as submit_predict
 
 MF = "DynaCLR-2D-MIP-BagOfChannels"
@@ -57,3 +59,42 @@ def test_predict_cmd_no_optional_flags():
     assert "--markers" not in cmd
     assert "--no-labelfree" not in cmd
     assert "--z-range" not in cmd
+
+
+def test_predict_batch_cli_forwards_prediction_flags():
+    result = CliRunner().invoke(
+        submit_predict.main,
+        [
+            "-c",
+            "c.yml",
+            "--model-family",
+            MF,
+            "--run",
+            RUN,
+            "--ckpt-name",
+            "last",
+            "--checkpoint",
+            "/c.ckpt",
+            "--datasets-root",
+            "/data",
+            "--z-range",
+            "15",
+            "45",
+            "--z-reduction",
+            "mip",
+            "--reference-pixel-size",
+            "0.1494",
+            "--batch-size",
+            "64",
+            "--skip-preflight",
+            "--print-cmd",
+        ],
+    )
+
+    assert result.exit_code == 0, result.output
+    args = result.output.splitlines()
+    z_range_index = args.index("--z-range")
+    assert args[z_range_index + 1 : z_range_index + 3] == ["15", "45"]
+    assert args[args.index("--z-reduction") + 1] == "mip"
+    assert args[args.index("--reference-pixel-size") + 1] == "0.1494"
+    assert args[args.index("--batch-size") + 1] == "64"
diff --git a/applications/dynaclr/tools/eval.sbatch b/applications/dynaclr/tools/eval.sbatch
index 1dfd1e484..c27f75bd6 100644
--- a/applications/dynaclr/tools/eval.sbatch
+++ b/applications/dynaclr/tools/eval.sbatch
@@ -10,6 +10,7 @@
 #   $2 model_family    embedding-tree key
 #   $3 run             embedding-tree key
 #   $4 ckpt_name       embedding-tree key
+#   $5 datasets_root   base under which datasets live
 #
 #SBATCH --job-name=dynaclr_eval
 #SBATCH --nodes=1
@@ -31,9 +32,11 @@ EVAL_CONFIG="$1"
 MODEL_FAMILY="$2"
 RUN="$3"
 CKPT_NAME="$4"
+DATASETS_ROOT="$5"
 
 uv run --project "$WORKSPACE_DIR" dynaclr eval \
   --eval-config "$EVAL_CONFIG" \
   --model-family "$MODEL_FAMILY" \
   --run "$RUN" \
-  --ckpt-name "$CKPT_NAME"
+  --ckpt-name "$CKPT_NAME" \
+  --datasets-root "$DATASETS_ROOT"
diff --git a/applications/dynaclr/tools/predict.sbatch b/applications/dynaclr/tools/predict.sbatch
index fc35aad4a..03cb28cd7 100644
--- a/applications/dynaclr/tools/predict.sbatch
+++ b/applications/dynaclr/tools/predict.sbatch
@@ -10,6 +10,12 @@
 #   $4 ckpt_name       checkpoint label (last | epochN-stepM)
 #   $5 datasets_root   base under which datasets live
 #   $6 checkpoint      (optional) explicit checkpoint path; empty = derive from run dir
+#   $7 markers         (optional) comma-separated marker subset
+#   $8 z_start         (optional) first --z-range value
+#   $9 z_end           (optional) second --z-range value
+#   $10 z_reduction    (optional) mip | center
+#   $11 reference_pixel_size (optional) microns per pixel
+#   $12 batch_size     (optional) prediction batch size
 #
 #SBATCH --job-name=dynaclr_predict
 #SBATCH --nodes=1
@@ -33,6 +39,11 @@ CKPT_NAME="$4"
 DATASETS_ROOT="$5"
 CHECKPOINT="${6:-}"
 MARKERS="${7:-}"
+Z_START="${8:-}"
+Z_END="${9:-}"
+Z_REDUCTION="${10:-}"
+REFERENCE_PIXEL_SIZE="${11:-}"
+BATCH_SIZE="${12:-}"
 
 CKPT_FLAG=()
 if [ -n "$CHECKPOINT" ]; then
@@ -41,9 +52,28 @@ fi
 
 MARKER_FLAG=()
 if [ -n "$MARKERS" ]; then
-  MARKER_FLAG=(--markers)
   IFS=',' read -ra _M <<< "$MARKERS"
-  MARKER_FLAG+=("${_M[@]}")
+  for _m in "${_M[@]}"; do
+    MARKER_FLAG+=(--markers "$_m")
+  done
+fi
+
+PREDICT_FLAG=()
+if [ -n "$Z_START" ] || [ -n "$Z_END" ]; then
+  if [ -z "$Z_START" ] || [ -z "$Z_END" ]; then
+    echo "Both z_start and z_end are required." >&2
+    exit 2
+  fi
+  PREDICT_FLAG+=(--z-range "$Z_START" "$Z_END")
+fi
+if [ -n "$Z_REDUCTION" ]; then
+  PREDICT_FLAG+=(--z-reduction "$Z_REDUCTION")
+fi
+if [ -n "$REFERENCE_PIXEL_SIZE" ]; then
+  PREDICT_FLAG+=(--reference-pixel-size "$REFERENCE_PIXEL_SIZE")
+fi
+if [ -n "$BATCH_SIZE" ]; then
+  PREDICT_FLAG+=(--batch-size "$BATCH_SIZE")
 fi
 
 srun uv run --project "$WORKSPACE_DIR" dynaclr predict-batch \
@@ -54,4 +84,5 @@ srun uv run --project "$WORKSPACE_DIR" dynaclr predict-batch \
   --datasets-root "$DATASETS_ROOT" \
   "${CKPT_FLAG[@]}" \
   "${MARKER_FLAG[@]}" \
+  "${PREDICT_FLAG[@]}" \
   --num-workers 0
diff --git a/packages/viscy-data/src/viscy_data/triplet.py b/packages/viscy-data/src/viscy_data/triplet.py
index b58e35391..e17bbfbb2 100644
--- a/packages/viscy-data/src/viscy_data/triplet.py
+++ b/packages/viscy-data/src/viscy_data/triplet.py
@@ -602,14 +602,11 @@ def __init__(
                 f"Extracting {self.initial_yx_patch_size} px patches "
                 f"and resizing to {final_yx_patch_size} px."
             )
-            scale_yx = (
-                final_yx_patch_size[0] / self.initial_yx_patch_size[0],
-                final_yx_patch_size[1] / self.initial_yx_patch_size[1],
-            )
             extra_transforms.append(
                 BatchedZoomd(
                     keys=list(self.source_channel),
-                    scale_factor=(1.0, *scale_yx),
+                    scale_factor=None,
+                    size=(extraction_width, *final_yx_patch_size),
                     mode="nearest-exact",
                 )
             )
diff --git a/packages/viscy-data/tests/test_triplet.py b/packages/viscy-data/tests/test_triplet.py
index 2f2de9866..736252a93 100644
--- a/packages/viscy-data/tests/test_triplet.py
+++ b/packages/viscy-data/tests/test_triplet.py
@@ -346,13 +346,9 @@ def test_timepoint_statistics_resolved_in_triplet_dataset(
 def test_reference_pixel_size_rescale_output_shape(preprocessed_hcs_dataset, tracks_hcs_dataset, reference_pixel_size):
     """reference_pixel_size rescale must land exactly on final_yx_patch_size.
 
-    The fixture pixel size is 1.0 µm/px, so a reference of 1.08 makes the naive
-    ``initial = round(final * scale)`` land on an odd number (35). The dataset
-    extracts a centered window of width ``2 * (initial // 2)`` = 34, one pixel
-    short, so a scale-factor resize would undershoot to 31 rather than 32 and the
-    datamodule's spatial-shape check would raise. Rounding the extraction size to
-    an even number keeps the resize exact. 1.3186 mirrors the SEC61B_DENV run
-    (0.1494 / 0.1133).
+    The extraction size is kept even for centered crops, while interpolation
+    receives ``final_yx_patch_size`` explicitly to avoid scale-factor flooring.
+    1.3186 mirrors the SEC61B_DENV run (0.1494 / 0.1133).
     """
     z_range = (4, 9)
     final_yx = (32, 32)
diff --git a/packages/viscy-models/src/viscy_models/contrastive/__init__.py b/packages/viscy-models/src/viscy_models/contrastive/__init__.py
index 1d3cc12e7..9d63a504a 100644
--- a/packages/viscy-models/src/viscy_models/contrastive/__init__.py
+++ b/packages/viscy-models/src/viscy_models/contrastive/__init__.py
@@ -1,7 +1,16 @@
 """Contrastive learning architectures."""
 
 from viscy_models.contrastive.encoder import ContrastiveEncoder, projection_mlp
-from viscy_models.contrastive.loss import NTXentHCL
+from viscy_models.contrastive.loss import NTXentHCL, NTXentLoss, TemporalStraighteningLoss
+from viscy_models.contrastive.predictor import Predictor
 from viscy_models.contrastive.resnet3d import ResNet3dEncoder
 
-__all__ = ["ContrastiveEncoder", "NTXentHCL", "ResNet3dEncoder", "projection_mlp"]
+__all__ = [
+    "ContrastiveEncoder",
+    "NTXentHCL",
+    "NTXentLoss",
+    "Predictor",
+    "ResNet3dEncoder",
+    "TemporalStraighteningLoss",
+    "projection_mlp",
+]
diff --git a/packages/viscy-models/src/viscy_models/contrastive/loss.py b/packages/viscy-models/src/viscy_models/contrastive/loss.py
index 1fa46cef6..d8203140e 100644
--- a/packages/viscy-models/src/viscy_models/contrastive/loss.py
+++ b/packages/viscy-models/src/viscy_models/contrastive/loss.py
@@ -3,6 +3,7 @@
 Provides:
 - ``NTXentLoss``: re-exported with a ``step()`` method for temperature scheduling.
 - ``NTXentHCL``: adds hard-negative concentration on top.
+- ``TemporalStraighteningLoss``: curvature loss on consecutive track embeddings.
 """
 
 from __future__ import annotations
@@ -10,9 +11,10 @@
 from typing import Literal
 
 import torch
+import torch.nn.functional as F
 from pytorch_metric_learning.losses import NTXentLoss as _NTXentLossBase
 from pytorch_metric_learning.utils import common_functions as c_f
-from torch import Tensor
+from torch import Tensor, nn
 
 from viscy_models.schedule import cosine_anneal
 
@@ -184,3 +186,62 @@ def _compute_loss(
                 }
             }
         return self.zero_losses()
+
+
+class TemporalStraighteningLoss(nn.Module):
+    """Curvature loss straightening latent trajectories (Wang et al. 2026).
+
+    Penalizes turning between consecutive latent velocity vectors along a track,
+    ``L_curv = mean_t (1 - cos(v_t, v_{t+1}))`` with ``v_t = z_{t+1} - z_t``.
+    Minimizing it makes successive velocities parallel, so the trajectory
+    approaches a straight line and Euclidean distance tracks progression.
+
+    Acts directly on the encoder embedding ``z`` (no projection head): cosine is
+    scale-invariant, so no normalization is applied inside the loss. Composes
+    additively with any anti-collapse objective (e.g. NT-Xent).
+
+    Parameters
+    ----------
+    eps : float
+        Numerical floor for the cosine similarity denominator. Default: 1e-8.
+    """
+
+    def __init__(self, eps: float = 1e-8) -> None:
+        super().__init__()
+        self.eps = eps
+
+    def step(self, epoch: int) -> None:
+        """No-op; the loss weight schedule lives on the LightningModule."""
+
+    def forward(self, z_seq: Tensor, valid: Tensor | None = None) -> Tensor:
+        """Compute mean curvature over a batch of consecutive-frame sequences.
+
+        Parameters
+        ----------
+        z_seq : Tensor
+            Embeddings of shape ``(B, K, D)`` with ``K >= 3`` consecutive frames
+            of the same track in time order.
+        valid : Tensor or None
+            Boolean mask of shape ``(B,)`` selecting samples with a real
+            (non-placeholder) sequence. Invalid samples are dropped. When every
+            sample is invalid, a graph-connected zero (``(z_seq * 0).sum()``) is
+            returned so DDP all-reduce stays balanced. Default: None (all used).
+
+        Returns
+        -------
+        Tensor
+            Scalar curvature loss in ``[0, 2]``.
+        """
+        if z_seq.size(1) < 3:
+            raise ValueError(f"z_seq must have K >= 3 frames, got K={z_seq.size(1)}")
+        v = z_seq[:, 1:, :] - z_seq[:, :-1, :]  # (B, K-1, D) velocities
+        v1, v2 = v[:, :-1, :], v[:, 1:, :]  # (B, K-2, D) consecutive pairs
+        cos = F.cosine_similarity(v1, v2, dim=-1, eps=self.eps)  # (B, K-2)
+        per_sample = (1.0 - cos).mean(dim=1)  # (B,)
+        if valid is not None:
+            per_sample = per_sample[valid]
+            if per_sample.numel() == 0:
+                # Graph-connected zero: keeps every DDP rank producing a
+                # gradient-tracked scalar so all-reduce stays balanced.
+                return (z_seq * 0.0).sum()
+        return per_sample.mean()
diff --git a/packages/viscy-models/src/viscy_models/contrastive/predictor.py b/packages/viscy-models/src/viscy_models/contrastive/predictor.py
new file mode 100644
index 000000000..f3a25b328
--- /dev/null
+++ b/packages/viscy-models/src/viscy_models/contrastive/predictor.py
@@ -0,0 +1,53 @@
+"""Latent predictor for temporal-dynamics learning (Wang et al. 2026 ``f_theta``).
+
+A single MLP shared across all cells that maps a latent state to its next state,
+``z_hat_{t+1} = predictor(z_t)``. Trained with a stop-gradient target
+(``L_pred = ||predictor(z_t) - sg(z_{t+1})||^2``). One shared predictor working
+across many cells is what rewards a *common* transition structure; the
+stop-gradient on the target prevents the trivial collapse where the encoder maps
+all frames to a constant.
+"""
+
+from __future__ import annotations
+
+import torch.nn as nn
+from torch import Tensor
+
+__all__ = ["Predictor"]
+
+
+class Predictor(nn.Module):
+    """Shared latent next-state predictor: an MLP ``dim -> hidden -> dim``.
+
+    Parameters
+    ----------
+    dim : int
+        Latent dimension (matches the encoder embedding, e.g. 768). Input and
+        output dimension.
+    hidden_dim : int or None
+        Hidden layer width. Defaults to ``dim`` when None.
+    """
+
+    def __init__(self, dim: int, hidden_dim: int | None = None) -> None:
+        super().__init__()
+        hidden_dim = hidden_dim or dim
+        self.net = nn.Sequential(
+            nn.Linear(dim, hidden_dim),
+            nn.GELU(),
+            nn.Linear(hidden_dim, dim),
+        )
+
+    def forward(self, z: Tensor) -> Tensor:
+        """Predict the next latent state.
+
+        Parameters
+        ----------
+        z : Tensor
+            Current latent states of shape ``(N, dim)``.
+
+        Returns
+        -------
+        Tensor
+            Predicted next latent states of shape ``(N, dim)``.
+        """
+        return self.net(z)
diff --git a/packages/viscy-models/tests/test_contrastive/test_predictor.py b/packages/viscy-models/tests/test_contrastive/test_predictor.py
new file mode 100644
index 000000000..b5b7a22fd
--- /dev/null
+++ b/packages/viscy-models/tests/test_contrastive/test_predictor.py
@@ -0,0 +1,71 @@
+"""Tests for the shared latent Predictor (next-state forecaster)."""
+
+import torch
+import torch.nn.functional as F
+
+from viscy_models.contrastive.predictor import Predictor
+
+
+def test_output_shape():
+    """Predictor maps (N, dim) -> (N, dim)."""
+    pred = Predictor(dim=32)
+    z = torch.randn(10, 32)
+    assert pred(z).shape == (10, 32)
+
+
+def test_custom_hidden_dim():
+    """Hidden width is configurable; output dim still equals input dim."""
+    pred = Predictor(dim=16, hidden_dim=64)
+    assert pred(torch.randn(4, 16)).shape == (4, 16)
+
+
+def test_learns_a_predictable_transition():
+    """A shared predictor can fit a common linear transition z_{t+1} = z_t + v."""
+    torch.manual_seed(0)
+    dim = 8
+    v = torch.randn(dim)  # one shared step common to all samples
+    pred = Predictor(dim=dim)
+    opt = torch.optim.Adam(pred.parameters(), lr=1e-2)
+    for _ in range(300):
+        z_t = torch.randn(64, dim)
+        z_next = z_t + v  # common transition
+        opt.zero_grad()
+        loss = F.mse_loss(pred(z_t), z_next.detach())
+        loss.backward()
+        opt.step()
+    # after fitting, prediction error on the common transition is small
+    z_t = torch.randn(256, dim)
+    final = F.mse_loss(pred(z_t), z_t + v)
+    assert final.item() < 0.05
+
+
+def test_random_transition_stays_high():
+    """A non-learnable (per-sample random) target cannot be fit -> loss stays high."""
+    torch.manual_seed(0)
+    dim = 8
+    pred = Predictor(dim=dim)
+    opt = torch.optim.Adam(pred.parameters(), lr=1e-2)
+    for _ in range(300):
+        z_t = torch.randn(64, dim)
+        z_next = torch.randn(64, dim)  # unrelated to z_t: no common rule
+        opt.zero_grad()
+        loss = F.mse_loss(pred(z_t), z_next.detach())
+        loss.backward()
+        opt.step()
+    z_t = torch.randn(256, dim)
+    z_next = torch.randn(256, dim)
+    final = F.mse_loss(pred(z_t), z_next)
+    # best a predictor can do for a zero-mean unit-var random target is ~var=1
+    assert final.item() > 0.5
+
+
+def test_stop_grad_target_no_grad_to_target():
+    """With a detached target, no gradient flows into the target tensor."""
+    dim = 8
+    pred = Predictor(dim=dim)
+    z_t = torch.randn(4, dim, requires_grad=True)
+    z_next = torch.randn(4, dim, requires_grad=True)
+    loss = F.mse_loss(pred(z_t), z_next.detach())
+    loss.backward()
+    assert z_next.grad is None  # target is stop-grad
+    assert z_t.grad is not None  # input side still receives gradient
diff --git a/packages/viscy-models/tests/test_contrastive/test_straightening.py b/packages/viscy-models/tests/test_contrastive/test_straightening.py
new file mode 100644
index 000000000..845a78bb6
--- /dev/null
+++ b/packages/viscy-models/tests/test_contrastive/test_straightening.py
@@ -0,0 +1,94 @@
+"""Tests for TemporalStraighteningLoss (curvature loss on latent tracks)."""
+
+import pytest
+import torch
+
+from viscy_models.contrastive.loss import TemporalStraighteningLoss
+
+
+def test_straight_trajectory_zero_curvature():
+    """A perfectly straight (constant-velocity) track has curvature ~0."""
+    # z_t = t * direction: consecutive velocities are identical -> cos=1 -> loss=0
+    t = torch.arange(5, dtype=torch.float32).view(1, 5, 1)
+    direction = torch.randn(1, 1, 8)
+    z_seq = t * direction  # (1, 5, 8), constant velocity = direction
+    loss = TemporalStraighteningLoss()(z_seq)
+    assert loss.item() == pytest.approx(0.0, abs=1e-5)
+
+
+def test_reversing_trajectory_max_curvature():
+    """A track that reverses direction each step has curvature ~2 (cos=-1)."""
+    # frames alternate between two points -> v_t and v_{t+1} are antiparallel
+    a = torch.zeros(1, 8)
+    b = torch.ones(1, 8)
+    z_seq = torch.stack([a, b, a, b, a], dim=1)  # (1, 5, 8)
+    loss = TemporalStraighteningLoss()(z_seq)
+    assert loss.item() == pytest.approx(2.0, abs=1e-5)
+
+
+def test_random_iid_points_curvature_near_1p5():
+    """I.i.d. random *points* give curvature ~1.5, not 1.0.
+
+    Consecutive velocities v_t = z_{t+1}-z_t and v_{t+1} = z_{t+2}-z_{t+1} share
+    the term z_{t+1} with opposite sign, so E[cos] = -0.5 -> E[1-cos] = 1.5.
+    (A value < 1.5 on real data indicates residual straightness — useful as a
+    reference: the DynaCLR probe measured ~1.29.)
+    """
+    torch.manual_seed(0)
+    z_seq = torch.randn(512, 3, 128)
+    loss = TemporalStraighteningLoss()(z_seq)
+    assert loss.item() == pytest.approx(1.5, abs=0.1)
+
+
+def test_valid_mask_drops_samples():
+    """The valid mask selects which samples contribute to the mean."""
+    straight = (torch.arange(3, dtype=torch.float32).view(1, 3, 1) * torch.ones(1, 1, 4)).repeat(2, 1, 1)
+    reversing = torch.stack([torch.zeros(2, 4), torch.ones(2, 4), torch.zeros(2, 4)], dim=1)
+    z_seq = torch.cat([straight, reversing], dim=0)  # (4, 3, 4): 2 straight, 2 reversing
+    loss_fn = TemporalStraighteningLoss()
+    # only the two straight samples counted -> ~0
+    valid = torch.tensor([True, True, False, False])
+    assert loss_fn(z_seq, valid=valid).item() == pytest.approx(0.0, abs=1e-5)
+    # only the two reversing samples counted -> ~2
+    valid = torch.tensor([False, False, True, True])
+    assert loss_fn(z_seq, valid=valid).item() == pytest.approx(2.0, abs=1e-5)
+
+
+def test_all_invalid_returns_graph_connected_zero():
+    """When no sample is valid, return a zero that still carries grad (DDP-safe)."""
+    z_seq = torch.randn(4, 3, 8, requires_grad=True)
+    valid = torch.zeros(4, dtype=torch.bool)
+    loss = TemporalStraighteningLoss()(z_seq, valid=valid)
+    assert loss.item() == 0.0
+    loss.backward()  # must not raise; grad graph is connected via new_zeros
+    assert z_seq.grad is not None
+
+
+def test_gradients_flow():
+    """Loss is differentiable w.r.t. the embeddings."""
+    z_seq = torch.randn(8, 3, 16, requires_grad=True)
+    loss = TemporalStraighteningLoss()(z_seq)
+    loss.backward()
+    assert z_seq.grad is not None
+    assert torch.isfinite(z_seq.grad).all()
+
+
+def test_k_less_than_three_raises():
+    """Curvature is undefined with fewer than 3 frames."""
+    z_seq = torch.randn(4, 2, 8)
+    with pytest.raises(ValueError, match="K >= 3"):
+        TemporalStraighteningLoss()(z_seq)
+
+
+def test_k_greater_than_three_sliding_window():
+    """K>3 averages curvature over the sliding stencil; straight stays ~0."""
+    z_seq = torch.arange(6, dtype=torch.float32).view(1, 6, 1) * torch.ones(1, 1, 4)
+    loss = TemporalStraighteningLoss()(z_seq)
+    assert loss.item() == pytest.approx(0.0, abs=1e-5)
+
+
+def test_step_is_noop():
+    """step() exists (for the on_train_epoch_start hook) and does nothing."""
+    loss_fn = TemporalStraighteningLoss()
+    loss_fn.step(0)
+    loss_fn.step(100)  # must not raise
diff --git a/packages/viscy-transforms/src/viscy_transforms/_zoom.py b/packages/viscy-transforms/src/viscy_transforms/_zoom.py
index afee93be4..062c8b710 100644
--- a/packages/viscy-transforms/src/viscy_transforms/_zoom.py
+++ b/packages/viscy-transforms/src/viscy_transforms/_zoom.py
@@ -15,16 +15,15 @@
 
 
 class BatchedZoom(Transform):
-    """Zoom (resize) a batched tensor by a scale factor.
+    """Resize a batched tensor by a scale factor or explicit size.
 
     Uses ``torch.nn.functional.interpolate`` for GPU-efficient resizing
     of batched 3D data. Supports various interpolation modes.
 
     Parameters
     ----------
-    scale_factor : float | tuple[float, float, float]
-        Multiplier for spatial size. If float, same factor is used for
-        all dimensions. If tuple, specifies (depth, height, width) factors.
+    scale_factor : float | tuple[float, float, float] | None
+        Multiplier for spatial size. Must be ``None`` when ``size`` is set.
     mode : str
         Interpolation algorithm. Options:
         - "nearest": Nearest neighbor interpolation
@@ -44,6 +43,9 @@ class BatchedZoom(Transform):
     antialias : bool
         If True, applies anti-aliasing when downsampling.
         Only effective for bilinear and bicubic modes. Default: False.
+    size : tuple[int, int, int] | None
+        Exact output ``(depth, height, width)``. Exactly one of ``size`` and
+        ``scale_factor`` must be provided.
 
     Returns
     -------
@@ -61,7 +63,7 @@ class BatchedZoom(Transform):
 
     def __init__(
         self,
-        scale_factor: float | tuple[float, float, float],
+        scale_factor: float | tuple[float, float, float] | None,
         mode: Literal[
             "nearest",
             "nearest-exact",
@@ -74,8 +76,12 @@ def __init__(
         align_corners: bool | None = None,
         recompute_scale_factor: bool | None = None,
         antialias: bool = False,
+        size: tuple[int, int, int] | None = None,
     ) -> None:
+        if (scale_factor is None) == (size is None):
+            raise ValueError("Provide exactly one of scale_factor or size.")
         self.scale_factor = scale_factor
+        self.size = size
         self.mode = mode
         self.align_corners = align_corners
         self.recompute_scale_factor = recompute_scale_factor
@@ -96,6 +102,7 @@ def __call__(self, sample: Tensor) -> Tensor:
         """
         return torch.nn.functional.interpolate(
             sample,
+            size=self.size,
             scale_factor=self.scale_factor,
             mode=self.mode,
             align_corners=self.align_corners,
@@ -113,9 +120,8 @@ class BatchedZoomd(MapTransform):
     ----------
     keys : Sequence[str]
         Keys of the data dictionary to apply zoom to.
-    scale_factor : float | tuple[float, float, float]
-        Multiplier for spatial size. If float, same factor is used for
-        all dimensions. If tuple, specifies (depth, height, width) factors.
+    scale_factor : float | tuple[float, float, float] | None
+        Multiplier for spatial size. Must be ``None`` when ``size`` is set.
     mode : str
         Interpolation algorithm. See :class:`BatchedZoom` for options.
     align_corners : bool | None
@@ -124,6 +130,9 @@ class BatchedZoomd(MapTransform):
         If True, recomputes scale_factor for interpolation. Default: None.
     antialias : bool
         If True, applies anti-aliasing when downsampling. Default: False.
+    size : tuple[int, int, int] | None
+        Exact output ``(depth, height, width)``. Exactly one of ``size`` and
+        ``scale_factor`` must be provided.
 
     Returns
     -------
@@ -148,7 +157,7 @@ class BatchedZoomd(MapTransform):
     def __init__(
         self,
         keys: Sequence[str],
-        scale_factor: float | tuple[float, float, float],
+        scale_factor: float | tuple[float, float, float] | None,
         mode: Literal[
             "nearest",
             "nearest-exact",
@@ -161,6 +170,7 @@ def __init__(
         align_corners: bool | None = None,
         recompute_scale_factor: bool | None = None,
         antialias: bool = False,
+        size: tuple[int, int, int] | None = None,
     ) -> None:
         super().__init__(keys)
         self.transform = BatchedZoom(
@@ -169,6 +179,7 @@ def __init__(
             align_corners=align_corners,
             recompute_scale_factor=recompute_scale_factor,
             antialias=antialias,
+            size=size,
         )
 
     def __call__(self, data: dict[str, Tensor]) -> dict[str, Tensor]:
diff --git a/packages/viscy-transforms/tests/test_zoom.py b/packages/viscy-transforms/tests/test_zoom.py
index 467202633..0d232a2b6 100644
--- a/packages/viscy-transforms/tests/test_zoom.py
+++ b/packages/viscy-transforms/tests/test_zoom.py
@@ -36,6 +36,21 @@ def test_batched_zoomd():
     assert result["label"].shape == expected_shape
 
 
+def test_batched_zoomd_exact_size_avoids_scale_rounding():
+    """Explicit size must avoid PyTorch scale-factor flooring (154 -> 159)."""
+    data = {"image": torch.rand(1, 1, 1, 154, 154)}
+    transform = BatchedZoomd(
+        keys=["image"],
+        scale_factor=None,
+        size=(1, 160, 160),
+        mode="nearest-exact",
+    )
+
+    result = transform(data)
+
+    assert result["image"].shape == (1, 1, 1, 160, 160)
+
+
 def test_batched_zoom_roundtrip():
     """Test roundtrip zoom (2x then 0.5x) returns close to original."""
     batch_size = 4
diff --git a/packages/viscy-utils/src/viscy_utils/evaluation/witness_gmm.py b/packages/viscy-utils/src/viscy_utils/evaluation/witness_gmm.py
index ae2b8c8e4..a61f8a0c0 100644
--- a/packages/viscy-utils/src/viscy_utils/evaluation/witness_gmm.py
+++ b/packages/viscy-utils/src/viscy_utils/evaluation/witness_gmm.py
@@ -148,3 +148,158 @@ def fit_gmm_labels(
         aic=float(gmm.aic(X)),
         bic_1comp=float(gmm_1.bic(X)),
     )
+
+
+def _gaussian_pdf(x: NDArray, mu: float, sigma: float) -> NDArray:
+    """1-D Gaussian density (sigma floored to avoid divide-by-zero)."""
+    sigma = max(float(sigma), 1e-9)
+    return np.exp(-0.5 * ((x - mu) / sigma) ** 2) / (sigma * np.sqrt(2 * np.pi))
+
+
+@dataclass
+class ControlAnchoredResult:
+    """Result of gating perturbed witness scores against a control-anchored baseline.
+
+    Models the perturbed scores as a two-component mixture
+    ``perturbed = pi * N(mu_c, sigma_c) + (1 - pi) * N(mu_r, sigma_r)`` where the
+    baseline component is **frozen to the control distribution** and only the
+    remodel component + mixing weight are fit. "Remodel" therefore means *excess
+    over the control baseline*, not one of two modes discovered within the
+    perturbed cells — the right frame when infection shifts the proportions of a
+    shared state-space rather than creating a distinct new state.
+
+    Parameters
+    ----------
+    mu_c, sigma_c : float
+        Frozen baseline Gaussian, fit to the control witness scores.
+    mu_r, sigma_r : float
+        Fitted remodel Gaussian (the perturbed-enriched excess).
+    pi_baseline : float
+        Estimated fraction of perturbed cells still explained by the baseline
+        component; ``1 - pi_baseline`` is the remodeled fraction. Weakly
+        identified when the shift is small — treat as approximate.
+    posterior : NDArray
+        Per perturbed cell: P(remodel | score), shape ``(n_perturbed,)``.
+    hard_label : NDArray
+        Per perturbed cell, int8: ``1`` where ``posterior >= threshold`` (remodel),
+        else ``-1``.
+    threshold : float
+        Posterior threshold, calibrated so the control false-positive rate equals
+        ``control_fp_target``.
+    control_fp : float
+        Realized control false-positive rate at ``threshold`` (≈ the target).
+    converged : bool
+        Whether the constrained EM converged.
+    """
+
+    mu_c: float
+    sigma_c: float
+    mu_r: float
+    sigma_r: float
+    pi_baseline: float
+    posterior: NDArray
+    hard_label: NDArray
+    threshold: float
+    control_fp: float
+    converged: bool
+
+
+def fit_control_anchored_labels(
+    perturbed_scores: NDArray,
+    control_scores: NDArray,
+    control_fp_target: float = 0.05,
+    max_iter: int = 200,
+    tol: float = 1e-6,
+) -> ControlAnchoredResult:
+    """Gate perturbed witness scores against a control-anchored baseline.
+
+    The baseline Gaussian is frozen to the control scores' mean/std; the perturbed
+    scores are then modeled as ``pi * baseline + (1 - pi) * remodel`` and only the
+    remodel Gaussian + mixing weight ``pi`` are fit by a constrained EM (the
+    baseline component's parameters never move). A perturbed cell's remodel
+    posterior is thresholded, with the threshold **calibrated so the control
+    false-positive rate equals ``control_fp_target``** — a controlled error rate
+    rather than a fixed posterior cut. Unlike :func:`fit_gmm_labels`, this never
+    abstains on a unimodal perturbed distribution: it measures the *excess* over
+    baseline, so a subtle proportion shift still yields labels (with a small
+    remodeled fraction).
+
+    Parameters
+    ----------
+    perturbed_scores : NDArray
+        Witness scores of the condition's perturbed cells, shape ``(n,)``.
+    control_scores : NDArray
+        Witness scores of the (time-matched) control-reference cells, shape ``(m,)``.
+    control_fp_target : float, optional
+        Target control false-positive rate; the posterior threshold is set so this
+        fraction of control cells is called remodel. By default 0.05.
+    max_iter : int, optional
+        Max constrained-EM iterations. By default 200.
+    tol : float, optional
+        Log-likelihood convergence tolerance. By default 1e-6.
+
+    Returns
+    -------
+    ControlAnchoredResult
+    """
+    p = np.asarray(perturbed_scores, dtype=np.float64).ravel()
+    c = np.asarray(control_scores, dtype=np.float64).ravel()
+
+    # Frozen baseline from control.
+    mu_c = float(c.mean())
+    sigma_c = float(c.std())
+
+    # Initialize the remodel component on the side of the perturbed cloud away from
+    # control (witness is negative = perturbed-leaning, so init below the mean).
+    mu_r = float(p.mean() - p.std())
+    sigma_r = float(p.std())
+    pi = 0.5  # baseline weight
+
+    prev_ll = -np.inf
+    converged = False
+    for _ in range(max_iter):
+        base = pi * _gaussian_pdf(p, mu_c, sigma_c)
+        rem = (1 - pi) * _gaussian_pdf(p, mu_r, sigma_r)
+        total = base + rem + 1e-300
+        resp_r = rem / total  # responsibility of the remodel component
+        # M-step — update ONLY the remodel component and the mixing weight.
+        nr = resp_r.sum()
+        if nr > 1e-6:
+            mu_r = float((resp_r * p).sum() / nr)
+            sigma_r = float(np.sqrt((resp_r * (p - mu_r) ** 2).sum() / nr))
+            sigma_r = max(sigma_r, 1e-6)
+        pi = float(1.0 - nr / len(p))
+        pi = min(max(pi, 1e-6), 1 - 1e-6)
+        ll = float(np.log(total).sum())
+        if abs(ll - prev_ll) < tol:
+            converged = True
+            break
+        prev_ll = ll
+
+    def _posterior(x: NDArray) -> NDArray:
+        base = pi * _gaussian_pdf(x, mu_c, sigma_c)
+        rem = (1 - pi) * _gaussian_pdf(x, mu_r, sigma_r)
+        return rem / (base + rem + 1e-300)
+
+    post_p = _posterior(p)
+    post_c = _posterior(c)
+    # Calibrate the threshold so the control FP rate matches the target: the
+    # (1 - target) quantile of the control posteriors.
+    threshold = float(np.quantile(post_c, 1.0 - control_fp_target))
+    control_fp = float((post_c >= threshold).mean())
+
+    hard_label = np.full(len(p), -1, dtype=np.int8)
+    hard_label[post_p >= threshold] = 1
+
+    return ControlAnchoredResult(
+        mu_c=mu_c,
+        sigma_c=sigma_c,
+        mu_r=mu_r,
+        sigma_r=sigma_r,
+        pi_baseline=pi,
+        posterior=post_p,
+        hard_label=hard_label,
+        threshold=threshold,
+        control_fp=control_fp,
+        converged=converged,
+    )
diff --git a/packages/viscy-utils/tests/test_witness_gmm.py b/packages/viscy-utils/tests/test_witness_gmm.py
index a8802cd3b..4547a74d6 100644
--- a/packages/viscy-utils/tests/test_witness_gmm.py
+++ b/packages/viscy-utils/tests/test_witness_gmm.py
@@ -52,3 +52,35 @@ def test_fit_gmm_labels_threshold_strictness():
     strict = fit_gmm_labels(scores, pos_threshold=0.99)
     # A stricter posterior bar yields no more confident positives than a lax one.
     assert (strict.hard_label == 1).sum() <= (lax.hard_label == 1).sum()
+
+
+def test_control_anchored_labels_calibrated_fp():
+    """Control-anchored gate: labels the excess over baseline and holds the
+    control false-positive rate near the target, without needing bimodality."""
+    from viscy_utils.evaluation.witness_gmm import fit_control_anchored_labels
+
+    rng = np.random.default_rng(0)
+    control = rng.normal(0.0, 1.0, 4000)
+    # Perturbed = mostly baseline + a shifted (remodel) minority (subtle, overlapping).
+    perturbed = np.concatenate([rng.normal(0.0, 1.0, 3000), rng.normal(-2.0, 1.0, 1000)])
+    res = fit_control_anchored_labels(perturbed, control, control_fp_target=0.05)
+    # FP rate is calibrated to the target.
+    assert abs(res.control_fp - 0.05) < 0.02
+    # Some perturbed cells are labeled remodel; remodel fraction is positive and < 1.
+    assert 0.0 < (1.0 - res.pi_baseline) < 1.0
+    assert (res.hard_label == 1).any()
+    # Remodel component sits below (more negative than) the control baseline.
+    assert res.mu_r < res.mu_c
+
+
+def test_control_anchored_no_shift_low_positives():
+    """When perturbed == control (no real shift), few cells clear the calibrated
+    threshold — the FP-calibrated cut keeps the positive rate near the target."""
+    from viscy_utils.evaluation.witness_gmm import fit_control_anchored_labels
+
+    rng = np.random.default_rng(1)
+    control = rng.normal(0.0, 1.0, 4000)
+    perturbed = rng.normal(0.0, 1.0, 4000)  # identical distribution
+    res = fit_control_anchored_labels(perturbed, control, control_fp_target=0.05)
+    # No real excess → perturbed positive rate stays near the target FP (~5%).
+    assert (res.hard_label == 1).mean() < 0.15