diff --git a/docs/tutorials/README.md b/docs/tutorials/README.md
index 3f7d0a6..f4072e2 100644
--- a/docs/tutorials/README.md
+++ b/docs/tutorials/README.md
@@ -11,6 +11,8 @@ These tutorials introduce the MicroSimulator modeling interface through runnable
 5. [Plasmid segregation, contacts, and conjugation](discrete-state-and-contacts.md)
 6. [Checkpoints, contact graphs, and quantitative analysis](analysis.md)
 7. [SimBOL circuit examples](simbol.md)
+8. [Microfluidic devices: walls, flow, and washout](microfluidics.md)
+9. [Solved flow: a pillar channel, Brinkman feedback, and the benchmarks](flow-solvers.md)
 
 The examples use `uv`, the `microsimulator` command, data-only checkpoints, and the standalone viewer. Each model selects its backend explicitly and can be run headlessly for batch experiments.
 
@@ -21,3 +23,5 @@ Teaching models are under [`examples/tutorials`](../../examples/tutorials). Scen
 The tutorials state numerical assumptions where they affect interpretation, including the meaning of cell length and volume, concentration dilution, time-step-dependent probabilities, signal units, and boundary conditions. For quantitative studies, follow the convergence and comparison guidance in each lesson rather than relying on viewer appearance alone.
 
 Readers comparing these models with the CellModeller wiki, legacy examples, or SimBOL sources can consult [tutorial sources and model translations](../compatibility/tutorial-source-provenance.md).
+
+The [nutrient validation study](nutrient-validation.md) supplies quantitative spatial-growth, conservation, and refinement evidence for the microfluidic stack.
diff --git a/docs/tutorials/flow-solvers.md b/docs/tutorials/flow-solvers.md
new file mode 100644
index 0000000..27ec792
--- /dev/null
+++ b/docs/tutorials/flow-solvers.md
@@ -0,0 +1,60 @@
+# Solved flow in a pillar channel
+
+The [pillar-channel model](../../examples/tutorials/pillar_channel.py) combines cylindrical walls, a depth-integrated flow calculation, attached founder lineages, and released daughters:
+
+```console
+uv run microsimulator view --model examples/tutorials/pillar_channel.py --seed 7 --dt 0.01 --backend metal --open
+```
+
+## Geometry and flow
+
+Mechanics uses continuous cylinders. The transport mask classifies voxel centers against those same cylinders, without shrinking their radii. Staircase walls approximate the physical geometry; curved-wall error and narrow gaps require refinement. Fluid-only, face-connected interpolation keeps solid concentrations out of cell exchange. A cell with no fluid sampling support is rejected rather than given artificial access through a wall.
+
+The channel's six-micrometer depth is represented by two three-micrometer fluid layers. The shallow solver uses one pressure per x/y column and depth-integrated conductance `H*m`, with base `m` proportional to `H²`. It routes flux around the pillars and lifts it conservatively to the signal grid:
+
+```python
+field, report = solve_flow_field(
+    grid, mean_inlet_speed=FLOW_SPEED,
+    mobility=gap_mobility(grid), simulation=simulation,
+)
+grid.velocity_field = field
+```
+
+The shallow approximation requires contiguous depth columns with a common floor. It does not resolve wall shear, flow within the gap, or the in-plane viscous stresses of a full Brinkman model. CPU, Metal, and CUDA implement the numerical solve independently, retaining accelerator vectors on device.
+
+## Attached cells and prescribed shedding
+
+Fixed founders remain attached. At division, the daughter nearer its lineage's adhesion site keeps the anchor and the farther daughter is released. This is a prescribed attachment rule, not a detachment prediction. Free daughters move at the sampled fluid velocity with a finite-aspect Jeffery orientation approximation, then undergo contact relaxation and model-defined outlet removal.
+
+Drift uses midpoint substeps bounded by grid displacement and angular change. The angular limit no longer clips the total rotation. The full growth/drift/contact split is still first order, so reduce the outer timestep and check mechanics convergence for quantitative motion.
+
+Only attached cells contribute stationary porous resistance:
+
+```python
+mobility = colony_mobility(
+    GRID, (cell for cell in step.cells if cell.fixed),
+    base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT,
+    averaging_radius=4.0,
+)
+field, report = solve_flow_field(
+    GRID, mean_inlet_speed=FLOW_SPEED,
+    mobility=mobility, simulation=step.simulation,
+)
+step.simulation.set_velocity_field(field)
+```
+
+The amount-conserving smoothing radius has physical units and stays fixed under grid refinement. The density is capped only inside the empirical resistance formula; biomass itself is conserved. Three small anchors should not be presented as a validated bulk biofilm blockage experiment. Freely advected daughters do not form a stationary matrix.
+
+Nutrient uptake in the tutorial equals the actual increment of biochemical biomass `B = pi*r²*(length + 2*r)` divided by the chosen yield. Backward Euler treats transport and affine loss implicitly; explicit cell uptake still requires an affordable step, and a rejected native biological step rolls back growth and chemistry.
+
+## Numerical evidence
+
+```console
+uv run python scripts/run_flow_benchmarks.py --backend cpu
+uv run python scripts/run_flow_benchmarks.py --backend metal
+uv run python scripts/run_flow_benchmarks.py --backend cpu --fine
+```
+
+The analytic suite checks plane and square ducts, a two-layer Brinkman channel, shallow routing, and agreement in a common thin-gap regime. Resolved flow uses flexible GMRES on a fixed velocity-pressure operator; inexact momentum solves are preconditioners. Its report includes freshly computed momentum and block residuals and physical divergence. Fine binary32 grids may need an explicitly looser tolerance, as the Brinkman benchmark documents.
+
+Use `solve_stokes_field` for resolved profiles when the mesh resolves the gap. `min_gap_voxels` identifies poorly resolved passages but does not certify accuracy. The [nutrient validation study](nutrient-validation.md) provides a separate, quantitative attached-population example with nutrient balance and grid, timestep, and flow-refresh sensitivity. It establishes numerical behavior under stated parameters, not biological calibration.
diff --git a/docs/tutorials/microfluidics.md b/docs/tutorials/microfluidics.md
new file mode 100644
index 0000000..a8d4d61
--- /dev/null
+++ b/docs/tutorials/microfluidics.md
@@ -0,0 +1,190 @@
+# Microfluidic devices: walls, flow, and washout
+
+This tutorial builds models that live inside devices: geometry that confines cells, blocks
+chemistry, and carries media. Four examples cover the range:
+
+| Model | Device | Demonstrates |
+| --- | --- | --- |
+| [`examples/culture_dish.py`](../../examples/culture_dish.py) | round dish | one inside-cylinder constraint as a dish |
+| [`examples/microfluidic_trap.py`](../../examples/microfluidic_trap.py) | trap + channel | flow, obstacles, drift, washout |
+| [`examples/tutorials/danino_clock.py`](../../examples/tutorials/danino_clock.py) | trap + channel | the full quorum clock in a device |
+| [`examples/tutorials/biopixel_trap.py`](../../examples/tutorials/biopixel_trap.py) | biopixel array trap | reported cavity, CAD layout, monolayer model |
+
+Run any of them live:
+
+```console
+uv run microsimulator view --model examples/microfluidic_trap.py --seed 42 --dt 0.02 --backend metal --open
+```
+
+## Walls that cells and chemistry both respect
+
+Mechanical walls are typed external constraints: infinite planes, spheres, axis-aligned
+boxes, and z-aligned cylinders, each with an inside or outside permitted region. A round
+culture dish is a single inside cylinder whose barrel is the wall and whose caps confine the
+monolayer:
+
+```python
+dish = CylinderConstraintInit()
+dish.radius = 30.0
+dish.half_height = 1.0
+dish.allowed_region = ConstraintRegion.INSIDE
+simulation.add_cylinder_constraint(dish)
+```
+
+Constraints alone are invisible to signals. The signal grid's obstacle mask closes every
+lattice face between fluid and solid voxels, so diffusion and advection stop at walls, and
+sampling near a wall renormalizes over fluid sites. Keeping the mask consistent with the
+constraints is an authoring concern, which the device helpers handle.
+
+## Devices from one description
+
+`microsimulator.microfluidics.TrapChannelDevice` describes an open-sided trap fed by a
+straight channel and projects that one description into every engine input:
+
+```python
+from microsimulator.microfluidics import TrapChannelDevice
+
+DEVICE = TrapChannelDevice(mean_flow_speed=20.0)
+DEVICE.add_constraints(simulation)                     # box walls for mechanics
+DEVICE.apply_to_grid(
+    grid,
+    inlet_values=[10.0],
+    outlet_values=[0.0],
+    simulation=simulation,
+)
+```
+
+`apply_to_grid` materializes the solid mask, fixed inlet and outlet boundaries on the y axis, and the numerically solved steady device flow on the grid's face-staggered velocity field (see the next section). Passing the model's `Simulation` makes the solve execute through the backend selected by the runner. Flow runs through the channel, circulates weakly at the open trap face, and the dead-end trap exchanges with the channel chiefly by diffusion in this model.
+
+## Flow on signals and on cells
+
+The velocity field advects every signal with conservative upwind face fluxes under both
+integrators. Cells feel the same field through explicit drift: with
+`MechanicsConfig(flow_drift=True)`, the controller advects each non-fixed cell by the fluid
+velocity sampled at its endpoints before contact relaxation, so escaped cells travel down the
+channel and rods rotate in shear. Contact relaxation then resolves any overlap the drift
+produced against walls or neighbors.
+
+## Washout
+
+Cells that reach the end of the channel leave the system through plan removals:
+
+```python
+def _regulate(step):
+    divisions = DIVISION.requests(step)
+    washed = tuple(cell.id for cell in step.cells if abs(cell.position.y) > WASHOUT_Y)
+    if washed:
+        DIVISION.forget(step, washed)
+        divisions = tuple(r for r in divisions if r.parent_id not in washed)
+    return StepPlan(updates=..., divisions=divisions, removals=washed)
+```
+
+Removal keeps stable identifiers and lineage history, so analysis can count washout events
+and trace removed cells' ancestry from checkpoints.
+
+## Numerical flow and stationary resistance
+
+`microsimulator.flow` solves `div_xy(H*m*grad_xy(p)) = 0` with one pressure per depth column. The default mobility is proportional to H², so integrated flux has the required H³ gap dependence. Harmonic face conductance and conservative lifting provide the face field used by transport. The shallow solver requires contiguous columns above a common floor; full three-dimensional obstructions require `microsimulator.stokes`. Both normalize velocity to a prescribed inlet speed and execute on the selected native backend.
+
+```python
+from microsimulator.flow import colony_mobility, solve_flow_field
+
+mobility = colony_mobility(
+    grid, (cell for cell in cells if cell.fixed),
+    drag_coefficient=100.0, averaging_radius=4.0,
+)
+field, report = solve_flow_field(
+    grid, mean_inlet_speed=20.0, mobility=mobility, simulation=simulation,
+)
+simulation.set_velocity_field(field)
+```
+
+Stationary resistance uses explicitly attached cells. A trapped but freely moving population is not automatically a stationary porous matrix. The four examples apply that distinction; the pillar model supplies attached founders. Biomass is deposited conservatively over a fixed physical radius, and only the resistance formula caps density. The drag coefficient and smoothing radius require calibration; grid size and refresh interval require numerical sensitivity checks.
+
+For resolved wall profiles use `solve_stokes_field`, with the same field interface and adequately resolved gaps. It solves a fixed Stokes-Brinkman block system by flexible GMRES and reports true momentum/block residuals and divergence. `min_gap_voxels` is diagnostic rather than a guarantee. The [flow tutorial](flow-solvers.md) gives the solver assumptions and analytic checks, and [nutrient validation](nutrient-validation.md) measures spatial growth, conservation, and refinement effects.
+
+## A source-backed Prindle biopixel example
+
+The [`prindle.dwg` and `prindle.dxf` files](devices) supplied with this tutorial are associated with the sensing-array project reported by Prindle et al. in [Nature 481, 39–44 (2012)](https://www.nature.com/articles/nature10722). Their provenance is recorded beside the files. The repository does not assert that this drawing is the exact fabrication revision used for the published experiments.
+
+The example deliberately separates three kinds of information:
+
+| Basis | Values used or observed | Role in the example |
+| --- | --- | --- |
+| Published methods | trapping region 100 x 85 x 1.65 micrometers; 25-micrometer trap spacing; nominal arrays of 500 and 12,000 biopixels | source of the modeled cavity dimensions and context for the array scale |
+| Supplied CAD | 496 matching model-space `Layer-2` outlines in a 16 x 31 layout; raw outline size 0.110 x 0.100 drawing units; raw row pitch 0.125 | validates the supplied layout and its source-specific scale, but does not define cavity walls or layer thicknesses |
+| Model choices | one 100 x 85 x 1.65 cavity beside a 100 x 10 x 300 micrometer channel; 10-micrometer numerical walls; mean inlet speed 20 micrometers per model time unit; chosen nutrient, drag, and re-solve parameters | defines a qualitative single-trap simulation, not a calibrated reconstruction of the experimental device |
+
+The trapping-region dimensions and spacing come from the [published supplementary methods](https://media.springernature.com/original/springer-static/esm/art%3A10.1038%2Fnature10722/MediaObjects/41586_2012_BFnature10722_MOESM313_ESM.pdf), not from subtracting a guessed wall inset from the CAD. `BiopixelTrapDevice` therefore defaults to a 100 x 85 x 1.65 micrometer cavity. Its channel dimensions, wall thickness, and flow speed remain ordinary constructor parameters:
+
+```python
+from microsimulator.microfluidics import BiopixelTrapDevice
+
+DEVICE = BiopixelTrapDevice(mean_flow_speed=20.0)
+```
+
+### Reading the supplied CAD layout
+
+`microsimulator.masks` is a bounded, data-only reader for model-space `LWPOLYLINE` geometry. It returns drawing coordinates unchanged unless the caller provides an explicit, source-specific `unit_scale`:
+
+```python
+from microsimulator.masks import extract_rectangles, load_mask_polylines, match_rectangles
+
+polylines = load_mask_polylines("docs/tutorials/devices/prindle.dxf")
+raw_rectangles = extract_rectangles(polylines, layer="Layer-2")
+raw_traps = match_rectangles(raw_rectangles, 0.110, 0.100, tolerance=0.001)
+
+rectangles_um = extract_rectangles(polylines, layer="Layer-2", unit_scale=1000.0)
+traps_um = match_rectangles(rectangles_um, 110.0, 100.0, tolerance=1.0)
+```
+
+For this file, treating one drawing unit as one millimeter is an inference corroborated by the publication: the raw 0.100 outline dimension maps to the reported 100-micrometer trap dimension, and the raw 0.125 row pitch maps to that dimension plus the reported 25-micrometer spacing. The DXF also stores `$INSUNITS=1`; [Autodesk documents `INSUNITS` as automatic insertion-scaling metadata and code 1 as inches](https://help.autodesk.com/cloudhelp/2026/ENU/AutoCAD-Core/files/GUID-A58A87BB-482B-4042-A00A-EEF55A2B4FD8.htm), which does not reconcile with these feature sizes. The reader therefore does not infer physical units from this header or impose the conversion on other drawings.
+
+Both raw and scaled queries yield 496 outlines in a 16 x 31 layout. Their centers span 2.4 x 3.75 millimeters after the inferred conversion; the rows have a 125-micrometer pitch, while column pitches are 135, 160, or 172.5 micrometers. The published device is described nominally as having 500 biopixels, so the documented result preserves the distinction between the paper's nominal count and this file's exact count.
+
+With `include_blocks=True`, the reader also exposes geometry in unplaced block definitions and records each block name. It does not apply `INSERT` transforms. The supplied file contains substantial `Layer-5` block geometry, but without a process map the tutorial does not assign that layer a physical role or infer cross-layer registration from it.
+
+The executable example loads and checks this layout, then simulates one cavity using the independently published dimensions. That single-trap reduction assumes one selected local inlet condition; it does not assert uniform flow across the array, reproduce the array manifold, or include inter-trap coupling. Run it live:
+
+```console
+uv run microsimulator view --model examples/tutorials/biopixel_trap.py --seed 5 --dt 0.02 --backend metal --open
+```
+
+## Units and timescales
+
+Model lengths are expressed in micrometers. Only the 100 x 85 x 1.65 trapping region is taken from the published methods; the table above identifies the remaining geometry and transport inputs as model choices.
+
+Time is a model growth scale. `growth_rate` is the exponential rate of cell length, so `BASE_GROWTH_RATE = 1.0` doubles cylindrical length in `ln 2 ≈ 0.69` model time units. Biochemical biomass includes an end contribution and therefore does not obey that exact exponential law. Mapping that doubling to a biological duration, such as 30 minutes, is illustrative and would make one model time unit about 43 minutes; it is not a calibration performed by this example. Nutrient and AHL levels are dimensionless concentration scales set by their inlet values and coupling parameters.
+
+For the biopixel example's configured channel values, `U = 20`, `L = 100`, and `D = 40` give a nominal channel-scale Péclet number `U L / D = 50`. That number characterizes this model only. Velocity is nonuniform, flow inside the dead-end cavity is much weaker, and no experimental flow or diffusivity measurements are fitted here, so the example makes no claim of experimental Péclet-number fidelity.
+
+The model also does not reproduce an experimentally established separation between transport and growth timescales. Its initial signal field is primed with inlet media, and its transport coefficients are chosen for a tractable tutorial run. Quantitative comparison with an experiment would require measured boundary conditions and material properties, grid and timestep convergence, and sensitivity analysis over the channel, transport, drag, and feedback parameters.
+
+## Numerical guidance
+
+- Choose `dt` so the largest per-step drift, `max_speed * dt`, stays below a cell radius;
+  `solve_flow_field` reports `max_speed`, and the trap examples use `dt = 0.02` with a mean
+  channel speed of 20.
+- Forward Euler enforces its stability bound from the per-site advective outflow; the trap
+  models select Crank-Nicolson.
+- The implicit solve's relative tolerance is the accuracy the step delivers: it asks for
+  that reduction of the residual the step starts with, so a model gets what it asked for
+  regardless of its concentration scale. These models keep the engine defaults.
+- Let the lattice of site centers cover every position a cell can reach, with about a voxel
+  of margin past each wall: contact relaxation lets a crowded cell press slightly into a
+  wall, and sampling outside the lattice is an error.
+- Keep the mechanics walls enclosing the solid mask. The device helpers voxelize
+  conservatively — a site is solid only when its whole voxel lies inside a wall — so the
+  voxel holding any reachable position stays fluid and a cell against a wall always has a
+  fluid site to sample. A hand-built mask needs the same rule; the
+  [pillar channel](flow-solvers.md) shows it for curved walls.
+- A sampling position whose whole stencil is solid raises an error rather than returning
+  zero.
+
+## Uptake and time integration
+
+The tutorials consume the actual biochemical biomass increment `Delta B / yield`, where `B = pi*r²*(length + 2*r)`. Division conserves this amount, which is distinct from geometric capsule volume. `growth_rate * B` is not the realized biomass rate under the length-growth law.
+
+Backward Euler is the baseline for their stiff transport and affine losses. Cellular sinks remain explicit; native biological failure restores growth, species, signals, time, and the prior solver report. Controller regulation, division callbacks, and mechanics are separate operations. No rejected step may be counted as successful growth. The Danino circuit is a qualitative example: its AHL secretion uses intracellular concentration times B, and its AiiA loss uses a conservatively smoothed enzyme amount. Its parameters and oscillations are not experimentally calibrated.
+
+The biopixel model explicitly uses signal absolute residual tolerance `1e-5` because binary32 noise at concentration 10 and its fine depth spacing prevents reliable convergence at `1e-6`. This is a model-scale numerical choice, not a biological accuracy claim.
diff --git a/docs/tutorials/nutrient-validation.md b/docs/tutorials/nutrient-validation.md
new file mode 100644
index 0000000..bc0b197
--- /dev/null
+++ b/docs/tutorials/nutrient-validation.md
@@ -0,0 +1,26 @@
+# Nutrient penetration and attached-population growth
+
+Run the controlled numerical study with:
+
+```console
+uv run python scripts/run_nutrient_benchmarks.py --backend cpu --output build/nutrient-cpu.json
+uv run python scripts/run_nutrient_benchmarks.py --backend metal --output build/nutrient-metal.json
+```
+
+The experiment places 45 fixed cells in a 20 x 40 x 4 micrometer channel. Initially nutrient is absent; inlet concentration is one and outlet concentration zero. Diffusion is 4 square micrometers per model time, mean inlet speed 0.3 micrometers per model time, maximum cylindrical-length growth rate 0.2 per model time, Monod constant 0.2, and biomass yield 0.5. Time is an illustrative model scale. There is no division, mechanics, or detachment in this experiment. These choices isolate nutrient delivery, growth, and stationary resistance.
+
+## Conservative coupling
+
+Each attached cell has conserved biochemical amount `B_i = pi*r_i²*(l_i + 2*r_i)`. A separable tent kernel K_i with physical radius 4 micrometers is integrated exactly over voxels and normalized so `sum_j K_ij V_j = 1`. Its support stays fixed as the grid is refined. Nutrient seen by a cell is `cbar_i = sum_j K_ij c_j V_j`.
+
+At step n, define `a_i = mu*pi*r_i²*l_i / (Y*(K_M + cbar_i))`. The grid loss is `lambda_j = sum_i a_i K_ij`. Native backward Euler advances conservative transport with this nonnegative loss, so the uptake assigned to cell i is `U_i = dt*a_i*sum_j K_ij*c_j^(n+1)*V_j`. The cell then gains `Delta B_i = Y*U_i`, with length updated accordingly. This is a first-order, semi-implicit Monod approximation with a lagged denominator and cylinder amount. It avoids the mesh-dependent singularity of a point sink and makes biomass gain equal the implicit nutrient loss to solver and rounding error. The four interactive tutorials separately exercise the native realized-growth rate instruction with explicit cell exchange.
+
+The flow resistance uses the same physical kernel and only attached biomass. Its coefficient is 40 in the normalized shallow closure. Flow is initially loaded and refreshed at the stated physical interval. All runs use the same physical domain, population, kernel, and parameters.
+
+## What is measured
+
+The script compares spacings 2, 1, and 0.5 micrometers, timesteps 0.04 and 0.02, and flow-refresh intervals 0.4 and 0.2 over eight model time units. It records total biomass gain, the first downstream crossing of half the inlet nutrient concentration, spatial growth, and the balance `nutrient remaining + biomass gained / yield = net boundary supply`. Boundary supply uses the engine's discrete ghost-center concentration convention and backward-Euler end-of-step fluxes. The chosen inlet boundary is therefore also refined with the grid.
+
+The JSON stores parameters, source commit, every profile, sensitivity comparisons, and pass/fail gates. CPU and Metal development runs gave about 31.00 biomass-volume units gained and a 5.71 micrometer half-concentration depth at spacing 1. Halving spacing changed biomass gain by 2.8% and penetration by 3.5%; halving the timestep changed gain by 0.44%; halving flow-refresh time changed gain by less than 0.001%. Maximum nutrient-balance error was below 0.008% of net supply. Regenerate these figures from the recorded commit before using them as release evidence.
+
+These results support a numerical demonstration of spatially limited growth. They do not calibrate nutrient yield, physical time, or the resistance law. The weak refresh sensitivity is specific to this slowly changing attached population, not a universal refresh recommendation. First-order upwind numerical diffusion remains: at spacing 1, its nominal scale `U*h/2 = 0.15` is about 3.8% of physical diffusion; it halves with spacing. More advective applications need their own refinement study or a higher-order transport method.
diff --git a/docs/tutorials/simbol.md b/docs/tutorials/simbol.md
index 17a4c70..9b3908d 100644
--- a/docs/tutorials/simbol.md
+++ b/docs/tutorials/simbol.md
@@ -116,7 +116,7 @@ These choices change trajectories relative to the generated callback scripts. A
 uv run microsimulator view \
   --model examples/tutorials/danino_clock.py \
   --seed 42 \
-  --dt 0.01 \
+  --dt 0.005 \
   --open
 ```
 
@@ -126,46 +126,68 @@ The example includes:
 - shared extracellular AHL and nutrient fields;
 - AHL-activated production with a third-order Hill response;
 - LuxI-dependent AHL production and AiiA-dependent AHL removal;
-- an AHL sink in the channel, nutrient replenishment in the trap, nutrient decay in the channel, and nutrient-limited growth;
+- a flow-fed channel that delivers nutrient, carries secreted AHL downstream, and washes out escaped cells;
+- nutrient-limited growth from the sampled local field;
 - stochastic daughter perturbations; and
-- the finite trap/channel obstacle geometry, expressed with typed plane and outside-sphere constraints.
+- the device geometry, flow field, obstacle mask, and inlet/outlet built from one `TrapChannelDevice` description in `microsimulator.microfluidics`.
 
 The biological motif is based on Danino et al., “A synchronized quorum of genetic clocks,” Nature 463, 326–330 (2010), as cited by the SimBOL model. The example equations and constants are a tutorial realization, not a reproduction of the paper’s experimental parameter inference.
 
-### Spatial field reactions
+### What the clock needs to run
 
-`CM_Danino.py` subclasses the legacy grid to apply an x-dependent AHL sink and an x-dependent nutrient source/decay field. It then reads nutrient to regulate growth. MicroSimulator represents those terms with the optional affine reaction field on `SignalGridSpec`:
+Three of the model's constants exist to make the clock a clock, and each is set against something measurable rather than by taste.
 
-```text
-dc[channel, x, y, z] / dt += source_rate[channel, x, y, z]
-                              - loss_rate[channel, x, y, z] * c[channel, x, y, z].
-```
+The Hill threshold must sit below the AHL the circuit can reach. LuxI and AiiA are driven by the same activation term, so their ratio, and with it the AHL where production balances enzymatic removal, is pinned by their decay constants at `2 * 0.3 / 1.2`. A threshold above that half is unreachable at any cell density and for any run length: the circuit sits at its basal state forever. `AHL_THRESHOLD` is set inside the window where the response is steep enough to oscillate.
+
+The rate scale sets the period. Growth defines the model's unit of time, so what matters is the clock's period relative to a doubling; `CLOCK_RATE` scales every rate constant together, which leaves the circuit's fixed points untouched and divides its period.
+
+AHL's diffusivity sets whether the trap oscillates as one. A patch of colony stays in phase with its neighbours only within about `sqrt(D * period)` of them, so with a period near one time unit and a trap 120 micrometers deep, `AHL_DIFFUSION` has to reach the order of ten thousand. Below that the trap breaks into independent patches.
 
-The coefficient arrays use the same signal-major, then x/y/z-major order as the concentration field. The model builds them once from the physical lattice coordinate `origin.x + x * spacing.x`; no Python callback runs during signal integration. With `outside = x < -60`, the declared coefficients are:
+AiiA's removal of AHL is a loss proportional to the AHL already present, which is a property of the field rather than of the cell, so the model rasterizes AiiA into an affine grid reaction each step and hands it to transport with `set_signal_reaction`. Transport takes a loss into its implicit diagonal and stays positive while the loss times the step is under two; the same removal scattered from the cells is explicit and needs half that step. What a synchronized pulse reaches in a packed trap is what sets `--dt` here.
 
-| Field and region             | source rate | loss rate |
-| ---------------------------- | ----------: | --------: |
-| AHL, inside trap             |           0 |         0 |
-| AHL, outside in channel      |           0 |         5 |
-| nutrient, inside trap        |          20 |         2 |
-| nutrient, outside in channel |           0 |       0.5 |
+A run at seed 42 measures the result: the trap is quiet through the colony's growth, first pulses once it holds about four thousand cells, and then pulses every 1.03 time units - 1.5 doubling times, the fast end of the 1.5 to 3 the paper reports - for as long as the run continues. The front and back halves of the trap rise and fall together, correlated 0.87 at zero lag. That is the quorum the circuit is named for: not a clock that each cell keeps, but one the population only starts once it is dense enough to talk to itself.
 
-Thus the inside nutrient equation is `dN/dt += 2(10 - N)`. The outside reaction is `dN/dt += -0.5 N`, and the outside AHL reaction is `dA/dt += -5 A`. Both fields begin at zero, so the nutrient reservoir develops dynamically rather than being installed as an initial condition.
+### Device flow and washout
 
-Before each biological step, the controller samples nutrient channel 1 at each cell center and applies the saturating growth law:
+`CM_Danino.py` subclasses the legacy grid to fake the channel with an x-dependent AHL sink and
+nutrient source field. The MicroSimulator model expresses the channel physically: a
+`TrapChannelDevice` projects one geometry description into box wall constraints, a signal-grid
+obstacle mask, a numerically solved steady flow field along the channel, and fixed inlet
+and outlet boundaries; as the colony packs the trap, the model re-solves the flow with the
+colony's Brinkman drag and swaps the field into the running simulation. Nutrient enters at
+the inlet at concentration 10 and is carried past the trap mouth; AHL secreted by the colony
+diffuses out of the trap and is advected downstream; walls block both diffusion and
+advection.
+
+Cells feel the same flow: the controller enables `flow_drift`, so a cell that escapes the
+trap is carried along the channel, and the regulation step removes any cell past the washout
+boundary with a `StepPlan` removal, forgetting its division target first.
+
+Before each biological step, the controller samples nutrient channel 1 at each cell center
+and applies the saturating growth law:
 
 ```text
 growth = nutrient / (5 + nutrient).
 ```
 
-The affine coefficients are immutable grid configuration and exact checkpoint state. CPU, Metal, and CUDA evaluate the same focused native operator; the model does not inject or compile an arbitrary voxel function at runtime.
+Growth and consumption are one loop: the coupled rate plan returns a nutrient sink of
+`growth_rate * cell_volume / NUTRIENT_YIELD` per cell, so a cell consumes in proportion to
+the growth the sampled field allows, and a packed trap draws down the field that feeds it.
+The yield sets the coupling strength; nutrient is an abstract limiting substrate in the
+inlet's concentration units.
+
+The device starts flooded with media, matching how a physical device is loaded before flow
+begins.
 
 ### Numerical interpretation
 
 Forward Euler evaluates transport, affine field reaction, and cell scatter from the old field and commits them together. Its preflight stability bound includes the largest local loss rate for each signal. This model selects Crank–Nicolson: spatial losses enter the implicit diagonal, fixed sources enter both trapezoidal halves, and cell-scattered AHL exchange remains an old-field explicit source. A converged negative result is still rejected because Crank–Nicolson is not positivity preserving for arbitrarily stiff steps.
 
+Crank–Nicolson accepts a step on its residual, and this model's cell exchange — AHL secreted into the field, nutrient drawn out of it — is small next to a nutrient background of ten. Convergence is therefore judged against the residual the step starts with rather than against the field, so the threshold does not grow with the background and a cell's contribution cannot fall under it. The model keeps the engine's default tolerances.
+
 ## Exercises
 
+- Vary `mean_flow_speed` and measure how the trap's AHL retention, and therefore the clock's synchronization, responds.
 - Run an inducer sweep with a data-only run manifest and compare reporter concentration at a fixed physical time and cell count.
 - Compare BBa_0004 with a self-repressed LacI equation derived directly from the summarized SBOL topology. Treat it as a different model, not a bug-free rerun of the generated script.
 - Restore the larger BBa_0003 grid and perform a grid-extent convergence check.
diff --git a/examples/microfluidic_trap.py b/examples/microfluidic_trap.py
new file mode 100644
index 0000000..2d812b7
--- /dev/null
+++ b/examples/microfluidic_trap.py
@@ -0,0 +1,178 @@
+"""A cell trap fed by a flowing channel.
+
+Fresh nutrient enters at the channel inlet, is carried past the trap mouth by
+the numerically solved steady device flow, and reaches the colony by diffusion
+through the open trap face. Only explicitly fixed cells contribute stationary
+flow resistance; this freely moving trap population is not a porous matrix. Cell
+growth follows Monod kinetics on the local nutrient level and consumes
+nutrient at a fixed yield, so the colony's growth pattern reflects the balance
+between flow supply, diffusion into the trap, and consumption by the cells
+already there.
+"""
+
+from __future__ import annotations
+
+from dataclasses import replace
+
+from microsimulator import (
+    CellInit,
+    CellUpdate,
+    ControllerStep,
+    CoupledRatePlan,
+    DivisionEvent,
+    GridShape,
+    MechanicsConfig,
+    ModelContext,
+    NativeController,
+    RatePlanBuilder,
+    SignalGridSpec,
+    SignalIntegrationKind,
+    Simulation,
+    StepPlan,
+    UniformLengthDivision,
+    Vec3,
+)
+from microsimulator.checkpoint import CheckpointBundle, JSONValue
+from microsimulator.flow import colony_mobility, gap_mobility, solve_flow_field
+from microsimulator.microfluidics import TrapChannelDevice
+
+MODEL_ID = "examples.microfluidic-trap"
+MODEL_VERSION = 4
+DIVISION = UniformLengthDivision(3.2, 3.8, jitter_z=False)
+
+FLOW_SPEED = 20.0
+DEVICE = TrapChannelDevice(mean_flow_speed=FLOW_SPEED)
+CELL_RADIUS = 0.5
+NUTRIENT_INLET = 10.0
+BASE_GROWTH_RATE = 1.0
+NUTRIENT_K = 5.0
+# Nutrient uses arbitrary concentration units. Each accepted step consumes
+# the actual increase of B = pi*r^2*(length + 2*r), divided by this yield.
+# The value is illustrative; penetration and growth require refinement checks.
+NUTRIENT_YIELD = 0.5
+WASHOUT_Y = DEVICE.channel_half_length - 10.0
+
+# Resistance feedback uses only fixed (attached) cells, with a physical
+# smoothing radius independent of the grid. Free cells do not form a matrix.
+RESOLVE_INTERVAL = 100
+DRAG_COEFFICIENT = 100.0
+
+
+def _grid(simulation: Simulation | None = None) -> SignalGridSpec:
+    shape = GridShape()
+    shape.x, shape.y, shape.z = 64, 72, 6
+    grid = SignalGridSpec()
+    grid.signal_count = 1
+    grid.shape = shape
+    # Two z layers span the trap's six-micrometer depth exactly, so its fluid
+    # volume is the device's rather than the half voxel of slack a coarser
+    # lattice would leave on each side, and the lattice still reaches a voxel
+    # past the walls: contact relaxation can press a crowded cell into a wall
+    # and briefly out through it, and sampling outside the lattice is an error.
+    grid.origin = Vec3(-140.0, -144.0, -7.5)
+    grid.spacing = Vec3(4.0, 4.0, 3.0)
+    grid.diffusion = [40.0]
+    grid.advection = [Vec3()]
+    grid.integration = SignalIntegrationKind.BACKWARD_EULER
+    device = DEVICE if simulation is not None else replace(DEVICE, mean_flow_speed=0.0)
+    device.apply_to_grid(
+        grid,
+        inlet_values=[NUTRIENT_INLET],
+        outlet_values=[0.0],
+        simulation=simulation,
+    )
+    return grid
+
+
+GRID = _grid()
+GAP_MOBILITY = gap_mobility(GRID)
+
+
+def _rate_plan() -> CoupledRatePlan:
+    rates = RatePlanBuilder()
+    uptake = -rates.cell_volume_change_rate() / NUTRIENT_YIELD
+    return rates.coupled_plan(0, 1, (), (uptake,))
+
+
+def _primed_levels(grid: SignalGridSpec) -> list[float]:
+    # The device is loaded flooded with fresh media before flow starts.
+    return [
+        NUTRIENT_INLET if solid == 0 else 0.0
+        for solid in grid.obstacles
+    ]
+
+
+def _nutrient_growth(simulation: Simulation, position: Vec3) -> float:
+    nutrient = max(0.0, simulation.sample_signals(position)[0])
+    return BASE_GROWTH_RATE * nutrient / (NUTRIENT_K + nutrient)
+
+
+def _regulate(step: ControllerStep) -> StepPlan:
+    if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0:
+        mobility = colony_mobility(
+            GRID, (cell for cell in step.cells if cell.fixed),
+            base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT
+        )
+        field, _ = solve_flow_field(
+            GRID,
+            mean_inlet_speed=FLOW_SPEED,
+            mobility=mobility,
+            simulation=step.simulation,
+        )
+        step.simulation.set_velocity_field(field)
+    divisions = DIVISION.requests(step)
+    washed = tuple(cell.id for cell in step.cells if abs(cell.position.y) > WASHOUT_Y)
+    if washed:
+        DIVISION.forget(step, washed)
+        divisions = tuple(request for request in divisions if request.parent_id not in washed)
+    return StepPlan(
+        updates=tuple(
+            CellUpdate(cell.id, growth_rate=_nutrient_growth(step.simulation, cell.position))
+            for cell in step.cells
+            if cell.id not in washed
+        ),
+        divisions=divisions,
+        removals=washed,
+    )
+
+
+def _divided(step: ControllerStep, event: DivisionEvent) -> None:
+    DIVISION.on_division(step, event)
+
+
+def build(context: ModelContext) -> NativeController:
+    simulation = context.simulation(reserved_capacity=10_000)
+    grid = _grid(simulation)
+    simulation.configure_signal_grid(grid, _primed_levels(grid))
+    simulation.set_coupled_rate_plan(_rate_plan())
+    DEVICE.add_constraints(simulation)
+
+    founder = CellInit()
+    founder.position = Vec3(DEVICE.trap_back_x - 5.0, 0.0, 0.0)
+    founder.length = 3.5
+    founder.radius = CELL_RADIUS
+    founder.growth_rate = 1.0
+    founder_id = simulation.add_cell(founder)
+    state: dict[str, JSONValue] = {"scope": "microfluidic-trap"}
+    DIVISION.initialize(state, context.rng, (founder_id,))
+    return NativeController(
+        simulation,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        rng=context.rng,
+        regulate=_regulate,
+        on_division=_divided,
+        mechanics=MechanicsConfig(flow_drift=True),
+        state=state,
+    )
+
+
+def resume(context: ModelContext, checkpoint: CheckpointBundle) -> NativeController:
+    del context
+    return NativeController.from_checkpoint(
+        checkpoint,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        regulate=_regulate,
+        on_division=_divided,
+    )
diff --git a/examples/tutorials/biopixel_trap.py b/examples/tutorials/biopixel_trap.py
new file mode 100644
index 0000000..900c089
--- /dev/null
+++ b/examples/tutorials/biopixel_trap.py
@@ -0,0 +1,208 @@
+"""One modeled biopixel from the Prindle sensing-array study.
+
+The supplied CAD contains 496 matching Layer-2 outlines in a 16 by 31 layout;
+the Nature article describes a nominal 500-biopixel device. The supplemental
+methods, rather than an inferred CAD wall inset, supply this model's 100 by 85
+by 1.65 micrometer trapping region. Loading the DXF validates its layout and a
+source-specific unit conversion but does not determine the cavity walls.
+
+The adjacent 100 by 10 by 300 micrometer channel, mean flow speed, numerical
+wall thickness, transport parameters, and Brinkman feedback are explicit model
+choices. The example simulates one trap under one chosen local boundary
+condition; it does not model hydraulic variation or coupling across the array.
+"""
+
+from __future__ import annotations
+
+from dataclasses import replace
+from pathlib import Path
+
+from microsimulator import (
+    CellInit,
+    CellUpdate,
+    ControllerStep,
+    CoupledRatePlan,
+    DivisionEvent,
+    GridShape,
+    MechanicsConfig,
+    ModelContext,
+    NativeController,
+    RatePlanBuilder,
+    SignalGridSpec,
+    SignalIntegrationKind,
+    Simulation,
+    StepPlan,
+    UniformLengthDivision,
+    Vec3,
+)
+from microsimulator.checkpoint import CheckpointBundle, JSONValue
+from microsimulator.flow import colony_mobility, gap_mobility, solve_flow_field
+from microsimulator.masks import (
+    MaskError,
+    MaskRectangle,
+    extract_rectangles,
+    load_mask_polylines,
+    match_rectangles,
+)
+from microsimulator.microfluidics import BiopixelTrapDevice
+
+MODEL_ID = "tutorials.biopixel-trap"
+MODEL_VERSION = 7
+DIVISION = UniformLengthDivision(3.2, 3.8, jitter_z=False)
+
+_MASK = Path(__file__).resolve().parents[2] / "docs" / "tutorials" / "devices" / "prindle.dxf"
+_MASK_UNIT_SCALE = 1000.0
+_MASK_OUTLINE = (110.0, 100.0)
+FLOW_SPEED = 20.0
+
+
+def _load_prindle_layout() -> tuple[MaskRectangle, ...]:
+    polylines = load_mask_polylines(_MASK)
+    rectangles = extract_rectangles(
+        polylines,
+        layer="Layer-2",
+        unit_scale=_MASK_UNIT_SCALE,
+    )
+    traps = match_rectangles(rectangles, *_MASK_OUTLINE, tolerance=1.0)
+    if len(traps) != 496:
+        raise MaskError(f"expected 496 Prindle layout outlines, found {len(traps)}")
+    return traps
+
+
+TRAP_OUTLINES = _load_prindle_layout()
+DEVICE = BiopixelTrapDevice(mean_flow_speed=FLOW_SPEED)
+CELL_RADIUS = 0.5
+WASHOUT_Y = DEVICE.channel_half_length - 10.0
+
+NUTRIENT_INLET = 10.0
+BASE_GROWTH_RATE = 1.0
+NUTRIENT_K = 5.0
+# Nutrient uses arbitrary concentration units. Each accepted step consumes
+# the actual increase of B = pi*r^2*(length + 2*r), divided by this yield.
+# The value is illustrative; penetration and growth require refinement checks.
+NUTRIENT_YIELD = 0.5
+
+# Resistance feedback uses only fixed (attached) cells, with a physical
+# smoothing radius independent of the grid. Free cells do not form a matrix.
+RESOLVE_INTERVAL = 100
+DRAG_COEFFICIENT = 100.0
+
+
+def _grid(simulation: Simulation | None = None) -> SignalGridSpec:
+    shape = GridShape()
+    shape.x, shape.y, shape.z = 42, 60, 14
+    grid = SignalGridSpec()
+    grid.signal_count = 1
+    grid.shape = shape
+    # The lattice of site centers covers every position a cell can reach, with
+    # a margin of one voxel past the floor and the far channel wall: contact
+    # relaxation lets a crowded cell press slightly into a wall, and sampling
+    # outside the lattice is an error. Two z layers span the cavity exactly;
+    # the resulting gap-height mobility ratio belongs to this model geometry.
+    grid.origin = Vec3(-100.0, -147.5, -0.4125)
+    grid.spacing = Vec3(5.0, 5.0, 0.825)
+    grid.diffusion = [40.0]
+    grid.advection = [Vec3()]
+    grid.integration = SignalIntegrationKind.BACKWARD_EULER
+    # Fine depth spacing at concentration 10 gives binary32 residual noise
+    # above 1e-6. State the absolute tolerance explicitly for this model scale.
+    grid.solver.absolute_tolerance = 1.0e-5
+    device = DEVICE if simulation is not None else replace(DEVICE, mean_flow_speed=0.0)
+    device.apply_to_grid(
+        grid,
+        inlet_values=[NUTRIENT_INLET],
+        outlet_values=[0.0],
+        simulation=simulation,
+    )
+    return grid
+
+
+GRID = _grid()
+GAP_MOBILITY = gap_mobility(GRID)
+
+
+def _rate_plan() -> CoupledRatePlan:
+    rates = RatePlanBuilder()
+    uptake = -rates.cell_volume_change_rate() / NUTRIENT_YIELD
+    return rates.coupled_plan(0, 1, (), (uptake,))
+
+
+def _primed_levels(grid: SignalGridSpec) -> list[float]:
+    # The device is loaded flooded with fresh media before flow starts.
+    return [NUTRIENT_INLET if solid == 0 else 0.0 for solid in grid.obstacles]
+
+
+def _nutrient_growth(simulation: Simulation, position: Vec3) -> float:
+    nutrient = max(0.0, simulation.sample_signals(position)[0])
+    return BASE_GROWTH_RATE * nutrient / (NUTRIENT_K + nutrient)
+
+
+def _regulate(step: ControllerStep) -> StepPlan:
+    if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0:
+        mobility = colony_mobility(
+            GRID, (cell for cell in step.cells if cell.fixed),
+            base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT
+        )
+        field, _ = solve_flow_field(
+            GRID,
+            mean_inlet_speed=FLOW_SPEED,
+            mobility=mobility,
+            simulation=step.simulation,
+        )
+        step.simulation.set_velocity_field(field)
+    divisions = DIVISION.requests(step)
+    washed = tuple(cell.id for cell in step.cells if abs(cell.position.y) > WASHOUT_Y)
+    if washed:
+        DIVISION.forget(step, washed)
+        divisions = tuple(request for request in divisions if request.parent_id not in washed)
+    return StepPlan(
+        updates=tuple(
+            CellUpdate(cell.id, growth_rate=_nutrient_growth(step.simulation, cell.position))
+            for cell in step.cells
+            if cell.id not in washed
+        ),
+        divisions=divisions,
+        removals=washed,
+    )
+
+
+def _divided(step: ControllerStep, event: DivisionEvent) -> None:
+    DIVISION.on_division(step, event)
+
+
+def build(context: ModelContext) -> NativeController:
+    simulation = context.simulation(reserved_capacity=20_000)
+    grid = _grid(simulation)
+    simulation.configure_signal_grid(grid, _primed_levels(grid))
+    simulation.set_coupled_rate_plan(_rate_plan())
+    DEVICE.add_constraints(simulation)
+
+    founder = CellInit()
+    founder.position = Vec3(DEVICE.trap_depth * 0.5, 0.0, DEVICE.trap_height * 0.5)
+    founder.length = 3.5
+    founder.radius = CELL_RADIUS
+    founder.growth_rate = 1.0
+    founder_id = simulation.add_cell(founder)
+    state: dict[str, JSONValue] = {"scope": "biopixel-trap"}
+    DIVISION.initialize(state, context.rng, (founder_id,))
+    return NativeController(
+        simulation,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        rng=context.rng,
+        regulate=_regulate,
+        on_division=_divided,
+        mechanics=MechanicsConfig(flow_drift=True, passes=2),
+        state=state,
+    )
+
+
+def resume(context: ModelContext, checkpoint: CheckpointBundle) -> NativeController:
+    del context
+    return NativeController.from_checkpoint(
+        checkpoint,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        regulate=_regulate,
+        on_division=_divided,
+    )
diff --git a/examples/tutorials/danino_clock.py b/examples/tutorials/danino_clock.py
index ac8b409..31c88bc 100644
--- a/examples/tutorials/danino_clock.py
+++ b/examples/tutorials/danino_clock.py
@@ -1,11 +1,21 @@
-"""SimBOL's Danino quorum-sensing clock, nutrient field, and trap geometry."""
+"""SimBOL's Danino quorum-sensing clock in a flow-fed microfluidic trap.
+
+Media flows along the channel with the numerically solved steady device flow:
+it delivers nutrient, carries secreted AHL downstream, and washes out cells
+that escape the trap. The colony feeds back on the flow: at a fixed cadence
+the model rasterizes the packed cells into a Brinkman drag field, re-solves
+the flow, and swaps the field into the running simulation.
+"""
 
 from __future__ import annotations
 
-import math
+from collections.abc import Sequence
+from dataclasses import replace
 
+import numpy as np
 from microsimulator import (
     CellInit,
+    CellSnapshot,
     CellUpdate,
     ControllerStep,
     CoupledRatePlan,
@@ -14,74 +24,120 @@
     MechanicsConfig,
     ModelContext,
     NativeController,
-    PlaneConstraintInit,
     RatePlanBuilder,
     SignalGridAffineReaction,
     SignalGridSpec,
     SignalIntegrationKind,
     Simulation,
-    SphereConstraintInit,
-    SphereRegion,
     StepPlan,
     UniformLengthDivision,
     Vec3,
 )
 from microsimulator.checkpoint import CheckpointBundle, JSONValue
+from microsimulator.flow import (
+    colony_mobility,
+    colony_species_density,
+    gap_mobility,
+    solve_flow_field,
+)
+from microsimulator.microfluidics import TrapChannelDevice
 
 MODEL_ID = "tutorials.danino-clock"
-MODEL_VERSION = 2
+MODEL_VERSION = 8
 DIVISION = UniformLengthDivision(3.2, 3.8, jitter_z=False)
 
-TRAP_OPEN_X = -60.0
-TRAP_BACK_X = 60.0
-TRAP_HALF_Y = 15.0
-TRAP_HALF_Z = 3.0
-CHANNEL_FAR_X = -100.0
-CHANNEL_HALF_LENGTH = 120.0
+FLOW_SPEED = 20.0
+DEVICE = TrapChannelDevice(mean_flow_speed=FLOW_SPEED)
 CELL_RADIUS = 0.5
+# Cells leave the analysis region well before the channel's end. Flow here is
+# compressed relative to growth, so a cell swept the full channel length would
+# divide several times on the way out and the downstream population would grow
+# without bound; the trap itself spans only the first fifteen micrometers.
+WASHOUT_Y = 40.0
 
-AHL_SINK_RATE = 5.0
-NUTRIENT_TARGET = 10.0
-NUTRIENT_SUPPLY_RATE = 2.0
-NUTRIENT_DECAY_RATE = 0.5
+# The clock's rate constants share one scale. Growth sets the model's unit of
+# time, so the scale is what places the clock's period relative to a doubling:
+# this value gives about two doubling times, the order Danino et al. report.
+CLOCK_RATE = 25.0
+# AiiA's removal of AHL, per unit of each. This is a loss proportional to the
+# AHL already there, so the model hands it to the grid as an affine reaction
+# rather than scattering it from the cells: transport takes a loss field into
+# its implicit diagonal, which stays positive while the loss times the step is
+# under two, where an explicit cell source of the same strength needs half
+# that step. The loss a synchronized pulse reaches in a packed trap is what
+# sets the time step this model runs at.
+AHL_REMOVAL = 1.0
+# The AHL concentration at which the Hill response is half activated. LuxI and
+# AiiA respond to the same activation, so their ratio - and with it the AHL
+# where production balances removal - is fixed by their decay constants at
+# 8 / AHL_REMOVAL times 0.3 / 1.2. A threshold above that is unreachable at any
+# cell density, for any run length, and the clock never starts; this one sits
+# where the response is steep enough to oscillate.
+AHL_THRESHOLD = 4.0
+# AHL crosses the trap in about the square of its width over this coefficient.
+# Below roughly ten thousand that exchange is slower than the clock's period
+# and the trap oscillates in independent patches instead of as one quorum.
+AHL_DIFFUSION = 10_000.0
+
+NUTRIENT_INLET = 10.0
 BASE_GROWTH_RATE = 1.0
 NUTRIENT_K = 5.0
+# Nutrient uses arbitrary concentration units. Each accepted step consumes
+# the actual increase of B = pi*r^2*(length + 2*r), divided by this yield.
+# The value is illustrative; penetration and growth require refinement checks.
+NUTRIENT_YIELD = 0.5
+
+# Resistance feedback uses only fixed (attached) cells, with a physical
+# smoothing radius independent of the grid. Free cells do not form a matrix.
+RESOLVE_INTERVAL = 400
+# How often AiiA is rasterized into the grid's AHL loss. The field follows the
+# clock, so refreshing it a few dozen times a period keeps it current while
+# leaving the per-step cost of building it in the noise.
+REMOVAL_INTERVAL = 10
+DRAG_COEFFICIENT = 100.0
 
 
-def _grid() -> SignalGridSpec:
+def _grid(simulation: Simulation | None = None) -> SignalGridSpec:
     shape = GridShape()
-    shape.x, shape.y, shape.z = 64, 72, 4
+    shape.x, shape.y, shape.z = 64, 72, 6
     grid = SignalGridSpec()
     grid.signal_count = 2
     grid.shape = shape
-    grid.origin = Vec3(-140.0, -144.0, -8.0)
-    grid.spacing = Vec3(4.0, 4.0, 4.0)
-    grid.diffusion = [40.0, 20.0]
+    # Two z layers span the trap's six-micrometer depth exactly, so its fluid
+    # volume is the device's rather than the half voxel of slack a coarser
+    # lattice would leave on each side, and the lattice still reaches a voxel
+    # past the walls: contact relaxation can press a crowded cell into a wall
+    # and briefly out through it, and sampling outside the lattice is an error.
+    grid.origin = Vec3(-140.0, -144.0, -7.5)
+    grid.spacing = Vec3(4.0, 4.0, 3.0)
+    grid.diffusion = [AHL_DIFFUSION, 20.0]
     grid.advection = [Vec3(), Vec3()]
-    grid.integration = SignalIntegrationKind.CRANK_NICOLSON
-    grid.solver.absolute_tolerance = 1.0e-12
-
-    site_count = shape.x * shape.y * shape.z
-    source_rates = [0.0] * (2 * site_count)
-    loss_rates = [0.0] * (2 * site_count)
-    for x in range(shape.x):
-        outside = grid.origin.x + x * grid.spacing.x < TRAP_OPEN_X
-        for y in range(shape.y):
-            for z in range(shape.z):
-                site = x * shape.y * shape.z + y * shape.z + z
-                if outside:
-                    loss_rates[site] = AHL_SINK_RATE
-                    loss_rates[site_count + site] = NUTRIENT_DECAY_RATE
-                else:
-                    source_rates[site_count + site] = NUTRIENT_SUPPLY_RATE * NUTRIENT_TARGET
-                    loss_rates[site_count + site] = NUTRIENT_SUPPLY_RATE
-    reaction = SignalGridAffineReaction()
-    reaction.source_rates = source_rates
-    reaction.loss_rates = loss_rates
-    grid.reaction = reaction
+    grid.integration = SignalIntegrationKind.BACKWARD_EULER
+    device = DEVICE if simulation is not None else replace(DEVICE, mean_flow_speed=0.0)
+    device.apply_to_grid(
+        grid,
+        inlet_values=[0.0, NUTRIENT_INLET],
+        outlet_values=[0.0, 0.0],
+        simulation=simulation,
+    )
     return grid
 
 
+GRID = _grid()
+GAP_MOBILITY = gap_mobility(GRID)
+
+
+def _primed_levels(grid: SignalGridSpec) -> list[float]:
+    # The device is loaded flooded with fresh media before flow starts; AHL
+    # starts at zero everywhere.
+    site_count = grid.shape.x * grid.shape.y * grid.shape.z
+    levels = [0.0] * (2 * site_count)
+    for site, solid in enumerate(grid.obstacles):
+        if solid == 0:
+            levels[site_count + site] = NUTRIENT_INLET
+    return levels
+
+
 def _rate_plan() -> CoupledRatePlan:
     rates = RatePlanBuilder()
     luxi = rates.maximum(rates.species(0), 0.0)
@@ -89,102 +145,45 @@ def _rate_plan() -> CoupledRatePlan:
     gfp = rates.maximum(rates.species(2), 0.0)
     ahl = rates.maximum(rates.signal(0), 0.0)
     ahl_cubed = ahl**3.0
-    hill = ahl_cubed / (8.0 + ahl_cubed)
-    activated = 0.02 + 8.0 * hill
+    hill = ahl_cubed / (AHL_THRESHOLD**3.0 + ahl_cubed)
+    activated = CLOCK_RATE * (0.02 + 8.0 * hill)
     return rates.coupled_plan(
         3,
         2,
         (
-            activated - 1.2 * luxi,
-            activated - 0.3 * aiia,
-            activated - 0.5 * gfp,
+            activated - CLOCK_RATE * 1.2 * luxi,
+            activated - CLOCK_RATE * 0.3 * aiia,
+            activated - CLOCK_RATE * 0.5 * gfp,
+        ),
+        (
+            CLOCK_RATE * 8.0 * luxi * rates.cell_volume(),
+            -rates.cell_volume_change_rate() / NUTRIENT_YIELD,
         ),
-        (8.0 * luxi - 4.0 * aiia * ahl, rates.constant(0.0)),
     )
 
 
-def _add_plane(
-    simulation: Simulation,
-    point: tuple[float, float, float],
-    normal: tuple[float, float, float],
-) -> None:
-    plane = PlaneConstraintInit()
-    plane.point = Vec3(*point)
-    plane.inward_normal = Vec3(*normal)
-    plane.coefficient = 1.0
-    simulation.add_plane_constraint(plane)
-
-
-def _add_sphere(simulation: Simulation, center: tuple[float, float, float]) -> None:
-    sphere = SphereConstraintInit()
-    sphere.center = Vec3(*center)
-    sphere.radius = CELL_RADIUS
-    sphere.coefficient = 1.0
-    sphere.allowed_region = SphereRegion.OUTSIDE
-    simulation.add_sphere_constraint(sphere)
-
-
-def _wall(
-    simulation: Simulation,
-    start: tuple[float, float, float],
-    end: tuple[float, float, float],
-) -> None:
-    delta = tuple(right - left for left, right in zip(start, end, strict=True))
-    length = math.sqrt(sum(value * value for value in delta))
-    count = max(2, math.ceil(length / CELL_RADIUS) + 1)
-    for index in range(count):
-        fraction = index / (count - 1)
-        center = (
-            start[0] + fraction * delta[0],
-            start[1] + fraction * delta[1],
-            start[2] + fraction * delta[2],
-        )
-        _add_sphere(simulation, center)
+_FLUID = np.asarray(GRID.obstacles, dtype=np.uint8) == 0
+_NO_SOURCES = [0.0] * (2 * GRID.site_count)
 
 
-def _add_trap(simulation: Simulation) -> None:
-    setback = 3.0
-    radius = CELL_RADIUS
-    _wall(
-        simulation,
-        (TRAP_OPEN_X + setback, -TRAP_HALF_Y - radius, 0.0),
-        (TRAP_BACK_X, -TRAP_HALF_Y - radius, 0.0),
-    )
-    _wall(
-        simulation,
-        (TRAP_OPEN_X + setback, TRAP_HALF_Y + radius, 0.0),
-        (TRAP_BACK_X, TRAP_HALF_Y + radius, 0.0),
-    )
-    _wall(
-        simulation,
-        (TRAP_BACK_X + radius, -TRAP_HALF_Y, 0.0),
-        (TRAP_BACK_X + radius, TRAP_HALF_Y, 0.0),
-    )
-    _add_plane(simulation, (CHANNEL_FAR_X, 0.0, 0.0), (1.0, 0.0, 0.0))
-    _wall(
-        simulation,
-        (TRAP_OPEN_X, -CHANNEL_HALF_LENGTH + 3.0, 0.0),
-        (TRAP_OPEN_X, -TRAP_HALF_Y, 0.0),
-    )
-    _wall(
-        simulation,
-        (TRAP_OPEN_X, TRAP_HALF_Y, 0.0),
-        (TRAP_OPEN_X, CHANNEL_HALF_LENGTH - 3.0, 0.0),
-    )
-    for y in (-TRAP_HALF_Y, TRAP_HALF_Y):
-        outer_y = y - radius if y < 0.0 else y + radius
-        _wall(
-            simulation,
-            (TRAP_OPEN_X + setback, outer_y, 0.0),
-            (TRAP_OPEN_X, outer_y, 0.0),
-        )
-        _wall(
-            simulation,
-            (TRAP_OPEN_X, outer_y, 0.0),
-            (TRAP_OPEN_X, y, 0.0),
-        )
-    _add_plane(simulation, (0.0, 0.0, TRAP_HALF_Z), (0.0, 0.0, -1.0))
-    _add_plane(simulation, (0.0, 0.0, -TRAP_HALF_Z), (0.0, 0.0, 1.0))
+def _ahl_removal_field(cells: Sequence[CellSnapshot]) -> SignalGridAffineReaction:
+    """Rasterize AiiA into the grid's first-order AHL loss.
+
+    Every cell removes AHL in proportion to its AiiA and to the AHL around it.
+    Concentration times conserved biochemical volume gives enzyme amount. A
+    fixed physical kernel distributes that amount conservatively, giving a loss
+    rate per unit time on the AHL field, which is what an affine reaction
+    carries. Nutrient takes no field reaction; its uptake follows growth and
+    stays a cell source.
+    """
+
+    aiia = np.asarray(colony_species_density(GRID, cells, species=1))
+    loss = np.zeros(2 * GRID.site_count, dtype=np.float64)
+    loss[: GRID.site_count] = np.where(_FLUID, CLOCK_RATE * AHL_REMOVAL * aiia, 0.0)
+    reaction = SignalGridAffineReaction()
+    reaction.source_rates = _NO_SOURCES
+    reaction.loss_rates = loss.tolist()
+    return reaction
 
 
 def _nutrient_growth(simulation: Simulation, position: Vec3) -> float:
@@ -193,12 +192,33 @@ def _nutrient_growth(simulation: Simulation, position: Vec3) -> float:
 
 
 def _regulate(step: ControllerStep) -> StepPlan:
+    if step.completed_steps % REMOVAL_INTERVAL == 0:
+        step.simulation.set_signal_reaction(_ahl_removal_field(step.cells))
+    if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0:
+        mobility = colony_mobility(
+            GRID, (cell for cell in step.cells if cell.fixed),
+            base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT
+        )
+        field, _ = solve_flow_field(
+            GRID,
+            mean_inlet_speed=FLOW_SPEED,
+            mobility=mobility,
+            simulation=step.simulation,
+        )
+        step.simulation.set_velocity_field(field)
+    divisions = DIVISION.requests(step)
+    washed = tuple(cell.id for cell in step.cells if abs(cell.position.y) > WASHOUT_Y)
+    if washed:
+        DIVISION.forget(step, washed)
+        divisions = tuple(request for request in divisions if request.parent_id not in washed)
     return StepPlan(
         updates=tuple(
             CellUpdate(cell.id, growth_rate=_nutrient_growth(step.simulation, cell.position))
             for cell in step.cells
+            if cell.id not in washed
         ),
-        divisions=DIVISION.requests(step),
+        divisions=divisions,
+        removals=washed,
     )
 
 
@@ -213,12 +233,13 @@ def _divided(step: ControllerStep, event: DivisionEvent) -> None:
 
 def build(context: ModelContext) -> NativeController:
     simulation = context.simulation(reserved_capacity=5_000, species_count=3)
-    simulation.configure_signal_grid(_grid())
+    grid = _grid(simulation)
+    simulation.configure_signal_grid(grid, _primed_levels(grid))
     simulation.set_coupled_rate_plan(_rate_plan())
-    _add_trap(simulation)
+    DEVICE.add_constraints(simulation)
 
     founder = CellInit()
-    founder.position = Vec3(TRAP_BACK_X - 5.0, 0.0, 0.0)
+    founder.position = Vec3(DEVICE.trap_back_x - 5.0, 0.0, 0.0)
     founder.length = 3.5
     founder.radius = CELL_RADIUS
     founder.growth_rate = 1.0
@@ -233,7 +254,7 @@ def build(context: ModelContext) -> NativeController:
         rng=context.rng,
         regulate=_regulate,
         on_division=_divided,
-        mechanics=MechanicsConfig(),
+        mechanics=MechanicsConfig(flow_drift=True),
         state=state,
     )
 
diff --git a/examples/tutorials/pillar_channel.py b/examples/tutorials/pillar_channel.py
new file mode 100644
index 0000000..b3e7340
--- /dev/null
+++ b/examples/tutorials/pillar_channel.py
@@ -0,0 +1,249 @@
+"""Colonies seeded on a pillar array in a flowing channel.
+
+The device is a monolayer channel crossed by a staggered array of cylindrical
+pillars - geometry with no analytic flow profile, so the field comes from the
+numerical solve: `solve_flow_field` routes the media around every pillar with
+per-voxel mass conservation, and the same solve re-runs at a fixed cadence
+with attached-cell resistance. The three anchors illustrate prescribed shedding;
+they do not establish a quantitative biofilm blockage or detachment model. Founder
+cells are adhered (fixed) in pillar wakes; each division keeps the mother
+attached and releases the daughter into the stream, which carries it between
+the pillars and washes it out at the end of the channel - a biofilm shedding
+cells into flow.
+"""
+
+from __future__ import annotations
+
+import math
+
+from microsimulator import (
+    BoxConstraintInit,
+    CellInit,
+    CellUpdate,
+    ConstraintRegion,
+    ControllerStep,
+    CoupledRatePlan,
+    CylinderConstraintInit,
+    DivisionEvent,
+    GridBoundaryKind,
+    GridShape,
+    MechanicsConfig,
+    ModelContext,
+    NativeController,
+    RatePlanBuilder,
+    SignalGridSpec,
+    SignalIntegrationKind,
+    Simulation,
+    StepPlan,
+    UniformLengthDivision,
+    Vec3,
+)
+from microsimulator.checkpoint import CheckpointBundle, JSONValue
+from microsimulator.flow import colony_mobility, gap_mobility, solve_flow_field
+
+MODEL_ID = "tutorials.pillar-channel"
+MODEL_VERSION = 3
+DIVISION = UniformLengthDivision(3.2, 3.8, jitter_z=False)
+
+CHANNEL_HALF_WIDTH = 40.0
+CHANNEL_HALF_LENGTH = 120.0
+CHANNEL_HALF_HEIGHT = 3.0
+PILLAR_RADIUS = 10.0
+PILLARS = ((-20.0, -60.0), (20.0, -60.0), (0.0, 0.0), (-20.0, 60.0), (20.0, 60.0))
+
+FLOW_SPEED = 20.0
+CELL_RADIUS = 0.5
+WASHOUT_Y = CHANNEL_HALF_LENGTH - 10.0
+# Adhesion sites in pillar wakes; the anchored cell of each lineage stays
+# within a cell length of its site.
+FOUNDER_SITES = ((-20.0, -46.0), (20.0, -46.0), (0.0, 14.0))
+
+NUTRIENT_INLET = 10.0
+BASE_GROWTH_RATE = 1.0
+NUTRIENT_K = 5.0
+# Nutrient uses arbitrary concentration units. Each accepted step consumes
+# the actual increase of B = pi*r^2*(length + 2*r), divided by this yield.
+# The value is illustrative; penetration and growth require refinement checks.
+NUTRIENT_YIELD = 0.5
+
+# Resistance feedback uses only fixed (attached) cells, with a physical
+# smoothing radius independent of the grid. Free cells do not form a matrix.
+RESOLVE_INTERVAL = 100
+DRAG_COEFFICIENT = 100.0
+
+
+def _in_pillar_core(px: float, py: float) -> bool:
+    # Classify centers against the same cylinder used by contact mechanics.
+    return any((px - x) ** 2 + (py - y) ** 2 < PILLAR_RADIUS**2 for x, y in PILLARS)
+
+
+def _grid(simulation: Simulation | None = None) -> SignalGridSpec:
+    shape = GridShape()
+    shape.x, shape.y, shape.z = 22, 60, 4
+    grid = SignalGridSpec()
+    grid.signal_count = 1
+    grid.shape = shape
+    grid.origin = Vec3(-42.0, -118.0, -4.5)
+    grid.spacing = Vec3(4.0, 4.0, 3.0)
+    grid.diffusion = [40.0]
+    grid.advection = [Vec3()]
+    grid.integration = SignalIntegrationKind.BACKWARD_EULER
+    obstacles = [0] * grid.site_count
+    for x in range(shape.x):
+        px = grid.origin.x + grid.spacing.x * x
+        for y in range(shape.y):
+            py = grid.origin.y + grid.spacing.y * y
+            for z in range(shape.z):
+                pz = grid.origin.z + grid.spacing.z * z
+                solid = (
+                    abs(px) >= CHANNEL_HALF_WIDTH
+                    or abs(pz) >= CHANNEL_HALF_HEIGHT
+                    or _in_pillar_core(px, py)
+                )
+                if solid:
+                    obstacles[(x * shape.y + y) * shape.z + z] = 1
+    grid.obstacles = obstacles
+    for name in ("y_lower", "y_upper"):
+        boundary = getattr(grid, name)
+        boundary.kind = GridBoundaryKind.FIXED
+        boundary.values = [NUTRIENT_INLET if name == "y_lower" else 0.0]
+        setattr(grid, name, boundary)
+    if simulation is not None:
+        field, _ = solve_flow_field(
+            grid,
+            mean_inlet_speed=FLOW_SPEED,
+            mobility=gap_mobility(grid),
+            simulation=simulation,
+        )
+        grid.velocity_field = field
+    return grid
+
+
+GRID = _grid()
+GAP_MOBILITY = gap_mobility(GRID)
+
+
+def _add_walls(simulation: Simulation) -> None:
+    chamber = BoxConstraintInit()
+    chamber.center = Vec3(0.0, 0.0, 0.0)
+    chamber.half_extents = Vec3(
+        CHANNEL_HALF_WIDTH, CHANNEL_HALF_LENGTH, CHANNEL_HALF_HEIGHT
+    )
+    chamber.coefficient = 1.0
+    chamber.allowed_region = ConstraintRegion.INSIDE
+    simulation.add_box_constraint(chamber)
+    for x, y in PILLARS:
+        pillar = CylinderConstraintInit()
+        pillar.center = Vec3(x, y, 0.0)
+        pillar.radius = PILLAR_RADIUS
+        pillar.half_height = CHANNEL_HALF_HEIGHT + 1.0
+        pillar.coefficient = 1.0
+        pillar.allowed_region = ConstraintRegion.OUTSIDE
+        simulation.add_cylinder_constraint(pillar)
+
+
+def _rate_plan() -> CoupledRatePlan:
+    rates = RatePlanBuilder()
+    uptake = -rates.cell_volume_change_rate() / NUTRIENT_YIELD
+    return rates.coupled_plan(0, 1, (), (uptake,))
+
+
+def _primed_levels(grid: SignalGridSpec) -> list[float]:
+    # The device is loaded flooded with fresh media before flow starts.
+    return [NUTRIENT_INLET if solid == 0 else 0.0 for solid in grid.obstacles]
+
+
+def _nutrient_growth(simulation: Simulation, position: Vec3) -> float:
+    nutrient = max(0.0, simulation.sample_signals(position)[0])
+    return BASE_GROWTH_RATE * nutrient / (NUTRIENT_K + nutrient)
+
+
+def _regulate(step: ControllerStep) -> StepPlan:
+    if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0:
+        mobility = colony_mobility(
+            GRID, (cell for cell in step.cells if cell.fixed),
+            base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT
+        )
+        field, _ = solve_flow_field(
+            GRID,
+            mean_inlet_speed=FLOW_SPEED,
+            mobility=mobility,
+            simulation=step.simulation,
+        )
+        step.simulation.set_velocity_field(field)
+    divisions = DIVISION.requests(step)
+    washed = tuple(cell.id for cell in step.cells if abs(cell.position.y) > WASHOUT_Y)
+    if washed:
+        DIVISION.forget(step, washed)
+        divisions = tuple(request for request in divisions if request.parent_id not in washed)
+    return StepPlan(
+        updates=tuple(
+            CellUpdate(cell.id, growth_rate=_nutrient_growth(step.simulation, cell.position))
+            for cell in step.cells
+            if cell.id not in washed
+        ),
+        divisions=divisions,
+        removals=washed,
+    )
+
+
+def _site_distance(position: Vec3) -> float:
+    return min(math.hypot(position.x - x, position.y - y) for x, y in FOUNDER_SITES)
+
+
+def _divided(step: ControllerStep, event: DivisionEvent) -> None:
+    DIVISION.on_division(step, event)
+    # Daughters inherit adhesion. The daughter nearer the adhesion site stays
+    # attached and the other is released into the stream; anchoring by site,
+    # not by daughter order, keeps the attached lineage at its wake instead of
+    # random-walking with every division (fixed cells are never moved by
+    # mechanics, so a walking anchor would end up inside a pillar).
+    if event.parent.fixed:
+        released = (
+            event.second
+            if _site_distance(event.first.position) <= _site_distance(event.second.position)
+            else event.first
+        )
+        step.simulation.set_cell_fixed(released.id, False)
+
+
+def build(context: ModelContext) -> NativeController:
+    simulation = context.simulation(reserved_capacity=10_000)
+    grid = _grid(simulation)
+    simulation.configure_signal_grid(grid, _primed_levels(grid))
+    simulation.set_coupled_rate_plan(_rate_plan())
+    _add_walls(simulation)
+
+    founder_ids: list[int] = []
+    for x, y in FOUNDER_SITES:
+        founder = CellInit()
+        founder.position = Vec3(x, y, 0.0)
+        founder.direction = Vec3(0.0, 1.0, 0.0)
+        founder.length = 3.5
+        founder.radius = CELL_RADIUS
+        founder.growth_rate = 1.0
+        founder.fixed = True
+        founder_ids.append(simulation.add_cell(founder))
+    state: dict[str, JSONValue] = {"scope": "pillar-channel"}
+    DIVISION.initialize(state, context.rng, tuple(founder_ids))
+    return NativeController(
+        simulation,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        rng=context.rng,
+        regulate=_regulate,
+        on_division=_divided,
+        mechanics=MechanicsConfig(flow_drift=True),
+        state=state,
+    )
+
+
+def resume(context: ModelContext, checkpoint: CheckpointBundle) -> NativeController:
+    del context
+    return NativeController.from_checkpoint(
+        checkpoint,
+        model_id=MODEL_ID,
+        model_version=MODEL_VERSION,
+        regulate=_regulate,
+        on_division=_divided,
+    )
diff --git a/python/tests/test_masks.py b/python/tests/test_masks.py
index 2f23021..5abe653 100644
--- a/python/tests/test_masks.py
+++ b/python/tests/test_masks.py
@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import math
+from itertools import pairwise
 from pathlib import Path
 
 import pytest
@@ -12,6 +13,8 @@
     match_rectangles,
 )
 
+_PRINDLE = Path(__file__).resolve().parents[2] / "docs" / "tutorials" / "devices" / "prindle.dxf"
+
 
 def test_rectangle_extraction_is_selective_and_explicitly_scaled() -> None:
     polylines = (
@@ -74,3 +77,65 @@ def test_block_definitions_are_opt_in_and_retain_their_name(tmp_path: Path) -> N
     block = next(polyline for polyline in with_blocks if polyline.block is not None)
     assert block.block == "FEATURE"
     assert block.vertices[0] == (10.0, 20.0)
+
+
+def test_prindle_mask_yields_the_documented_layout() -> None:
+    polylines = load_mask_polylines(_PRINDLE)
+    raw_rectangles = extract_rectangles(polylines, layer="Layer-2")
+    raw_traps = match_rectangles(raw_rectangles, 0.110, 0.100, tolerance=0.001)
+
+    # The supplemental methods report a 100-micrometer trap dimension with
+    # 25-micrometer spacing. In this particular drawing, its raw 0.100 outline
+    # dimension and 0.125 row pitch therefore corroborate a conversion from
+    # one drawing unit to one millimeter; this is evidence about this file,
+    # not a convention imposed on other DXF inputs.
+    assert len(raw_traps) == 496
+    assert all(
+        math.isclose(trap.width, 0.110, abs_tol=0.001)
+        and math.isclose(trap.height, 0.100, abs_tol=0.001)
+        for trap in raw_traps
+    )
+
+    rectangles = extract_rectangles(polylines, layer="Layer-2", unit_scale=1000.0)
+    traps = match_rectangles(rectangles, 110.0, 100.0, tolerance=1.0)
+    xs = sorted({round(trap.center[0], 1) for trap in traps})
+    ys = sorted({round(trap.center[1], 1) for trap in traps})
+
+    assert len(traps) == 496
+    assert len(xs) == 16
+    assert len(ys) == 31
+    assert math.isclose(ys[1] - ys[0], 125.0, abs_tol=0.1)
+    assert math.isclose(ys[-1] - ys[0], 3750.0, abs_tol=1.0)
+    column_pitches = sorted({round(right - left, 1) for left, right in pairwise(xs)})
+    assert column_pitches == [135.0, 160.0, 172.5]
+    assert math.isclose(xs[-1] - xs[0], 2400.0, abs_tol=1.0)
+
+
+def test_prindle_block_traversal_exposes_unplaced_layer_geometry() -> None:
+    polylines = load_mask_polylines(_PRINDLE, include_blocks=True)
+    blocks = {polyline.block for polyline in polylines if polyline.block is not None}
+    assert len(blocks) >= 2
+
+    layer5 = [
+        polyline
+        for polyline in polylines
+        if polyline.layer == "Layer-5" and polyline.block is not None
+    ]
+    assert len(layer5) > 3000
+    large_outlines = [
+        polyline
+        for polyline in layer5
+        if min(
+            max(vertex[0] for vertex in polyline.vertices)
+            - min(vertex[0] for vertex in polyline.vertices),
+            max(vertex[1] for vertex in polyline.vertices)
+            - min(vertex[1] for vertex in polyline.vertices),
+        )
+        >= 0.9
+    ]
+    assert len(large_outlines) >= 40
+
+    # These are unplaced block definitions. Without an accompanying process
+    # map, the test deliberately makes no claim about their fabrication role.
+    model_space = load_mask_polylines(_PRINDLE)
+    assert all(polyline.block is None for polyline in model_space)
diff --git a/python/tests/test_microfluidics.py b/python/tests/test_microfluidics.py
index 7e1a4d8..ede1437 100644
--- a/python/tests/test_microfluidics.py
+++ b/python/tests/test_microfluidics.py
@@ -3,9 +3,22 @@
 from __future__ import annotations
 
 import math
-
-from microsimulator import GridShape, SignalGridSpec, Vec3
+from pathlib import Path
+
+import pytest
+from microsimulator import (
+    BackendKind,
+    GridShape,
+    ModelContext,
+    SignalGridSpec,
+    SimulationController,
+    Vec3,
+    backend_available,
+)
 from microsimulator.microfluidics import BiopixelTrapDevice, TrapChannelDevice
+from microsimulator.runner import build_model
+
+_EXAMPLES = Path(__file__).resolve().parents[2] / "examples"
 
 
 def _grid() -> SignalGridSpec:
@@ -134,3 +147,80 @@ def test_planar_mask_wall_error_is_bounded_by_half_spacing() -> None:
             upper = spec.origin.x + (columns[-1] + 0.5) * h
             assert abs(lower - device.channel_far_x) <= h / 2 + 1e-5
             assert abs(upper - device.trap_back_x) <= h / 2 + 1e-5
+
+
+def test_trap_example_builds_steps_and_transports_nutrient() -> None:
+    model, _ = build_model(
+        _EXAMPLES / "microfluidic_trap.py",
+        ModelContext(BackendKind.CPU, 0, seed=11),
+    )
+    assert isinstance(model, SimulationController)
+    for _ in range(20):
+        model.step(0.02)
+
+    simulation = model.simulation
+    device = TrapChannelDevice()
+    channel_x = (device.channel_far_x + device.trap_open_x) * 0.5
+    upstream = simulation.sample_signals(Vec3(channel_x, -100.0, 0.0))[0]
+    trap_interior = simulation.sample_signals(Vec3(0.0, 0.0, 0.0))[0]
+    assert upstream > 5.0
+    assert trap_interior > 5.0
+    assert upstream >= trap_interior - 1.0e-3
+    with pytest.raises(ValueError, match="inside a grid obstacle"):
+        simulation.sample_signals(Vec3(0.0, 100.0, 0.0))
+    assert len(simulation.cells()) >= 1
+
+
+@pytest.mark.parametrize("backend", [BackendKind.METAL, BackendKind.CUDA])
+def test_trap_example_builds_its_initial_flow_on_the_selected_backend(
+    backend: BackendKind,
+) -> None:
+    if not backend_available(backend):
+        pytest.skip(f"{backend.name} backend is unavailable")
+    model, _ = build_model(
+        _EXAMPLES / "microfluidic_trap.py",
+        ModelContext(backend, 0, seed=13),
+    )
+    assert isinstance(model, SimulationController)
+    assert model.simulation.backend_info.kind == backend
+    checkpoint = model.simulation._checkpoint()
+    assert checkpoint.signal_grid is not None
+    assert checkpoint.signal_grid.spec.velocity_field is not None
+
+
+def test_biopixel_model_uses_reported_cavity_dimensions() -> None:
+    device = BiopixelTrapDevice(mean_flow_speed=20.0)
+
+    assert device.trap_width == 100.0
+    assert device.trap_depth == 85.0
+    assert device.trap_height == 1.65
+    assert device.channel_height == 10.0
+
+    # The CAD layout is tested independently in test_masks.py. These checks
+    # cover the published cavity size and the separately chosen model channel.
+    half = (2.5, 2.5, 0.825)
+    assert not device._solid(42.5, 0.0, 0.825, half)
+    assert device._solid(42.5, 0.0, 2.475, half)
+    assert not device._solid(-50.0, 0.0, 9.075, half)
+
+
+def test_biopixel_example_confines_a_monolayer_under_flow() -> None:
+    model, _ = build_model(
+        _EXAMPLES / "tutorials" / "biopixel_trap.py",
+        ModelContext(BackendKind.CPU, 0, seed=5),
+    )
+    assert isinstance(model, SimulationController)
+    # 110 steps crosses the model's Brinkman re-solve cadence at step 100, so
+    # the run exercises the colony-drag solve and the runtime field swap.
+    for _ in range(110):
+        model.step(0.02)
+
+    cells = model.simulation.cells()
+    assert len(cells) >= 2
+    for cell in cells:
+        assert 0.0 < cell.position.z < 1.65
+        assert -50.0 < cell.position.y < 50.0
+        assert cell.position.x < 95.0
+    checkpoint = model.simulation._checkpoint()
+    assert checkpoint.signal_grid is not None
+    assert checkpoint.signal_grid.spec.velocity_field is not None
diff --git a/python/tests/test_signal_model_runs.py b/python/tests/test_signal_model_runs.py
new file mode 100644
index 0000000..48ce826
--- /dev/null
+++ b/python/tests/test_signal_model_runs.py
@@ -0,0 +1,62 @@
+"""Multi-step runs of every model that integrates signals implicitly.
+
+A Crank-Nicolson step can fail long after a model starts: the solver's
+convergence threshold is compared against a residual whose floor rises with
+the magnitude of the field, so a model that converges from a near-empty grid
+can stop converging once its signals have grown. One step proves nothing about
+that. These runs advance each implicit model far enough for its field to
+develop, and require every step to converge and commit.
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import pytest
+from microsimulator import (
+    BackendKind,
+    ModelContext,
+    SimulationController,
+    build_model,
+)
+from microsimulator.checkpoint import JSONValue
+
+_ROOT = Path(__file__).resolve().parents[2]
+
+# One case per model that selects Crank-Nicolson, with the time step its
+# documentation recommends.
+_IMPLICIT_MODELS: tuple[tuple[str, dict[str, JSONValue], float], ...] = (
+    ("examples/tutorials/signaling.py", {"scenario": "communication"}, 0.02),
+    ("examples/tutorials/simbol_circuits.py", {"circuit": "bba_0003"}, 0.02),
+    ("examples/legacy/ex4_simpleCellCellSignaling.py", {}, 0.02),
+    ("examples/legacy/Tutorial_3/Tutorial_3.py", {}, 0.02),
+    ("examples/legacy/ACS2012/EdgeDetectorChamber.py", {}, 0.02),
+    ("examples/microfluidic_trap.py", {}, 0.02),
+    ("examples/tutorials/danino_clock.py", {}, 0.005),
+    ("examples/tutorials/pillar_channel.py", {}, 0.01),
+    ("examples/tutorials/biopixel_trap.py", {}, 0.02),
+)
+
+_STEPS = 200
+
+
+@pytest.mark.parametrize(("filename", "parameters", "dt"), _IMPLICIT_MODELS)
+def test_implicit_models_converge_over_a_long_run(
+    filename: str, parameters: dict[str, JSONValue], dt: float
+) -> None:
+    model, _ = build_model(
+        _ROOT / filename,
+        ModelContext(BackendKind.CPU, 0, seed=17, parameters=parameters),
+    )
+    assert isinstance(model, SimulationController)
+    assert model.simulation.has_signal_grid
+
+    for step in range(_STEPS):
+        model.step(dt)
+        report = model.simulation.last_signal_solve_report
+        assert report is not None, f"step {step} reported no signal solve"
+        assert report.converged, f"step {step} committed an unconverged field"
+
+    # The engine rejects a non-finite or negative field, so reaching here means
+    # every step committed a valid one.
+    assert model.simulation.cell_count > 0
diff --git a/python/tests/test_tutorials.py b/python/tests/test_tutorials.py
index 7d900b7..c23cab9 100644
--- a/python/tests/test_tutorials.py
+++ b/python/tests/test_tutorials.py
@@ -50,6 +50,8 @@
     ("plasmid_segregation.py", {"copies_per_cell": 10}, 0.001),
     ("conjugation.py", {"transfer_probability": 0.1}, 0.001),
     ("danino_clock.py", {}, 0.001),
+    ("biopixel_trap.py", {}, 0.001),
+    ("pillar_channel.py", {}, 0.001),
 )
 
 
@@ -116,7 +118,7 @@ def test_conjugation_tutorial_uses_current_contact_graph() -> None:
     assert model.simulation.cell(acceptor.id).cell_type == 2
 
 
-def test_danino_tutorial_declares_spatial_ahl_and_nutrient_reactions() -> None:
+def test_danino_tutorial_uses_device_flow_obstacles_and_washout() -> None:
     model, _ = build_model(
         _TUTORIALS / "danino_clock.py",
         ModelContext(BackendKind.CPU, 0, seed=31),
@@ -127,20 +129,48 @@ def test_danino_tutorial_declares_spatial_ahl_and_nutrient_reactions() -> None:
     checkpoint = model.simulation._checkpoint()
     assert checkpoint.signal_grid is not None
     spec = checkpoint.signal_grid.spec
-    assert spec.reaction is not None
-    sites = spec.site_count
-    outside = 0
-    inside = 20 * spec.shape.y * spec.shape.z
-    assert spec.origin.x + 19 * spec.spacing.x < -60.0
-    assert spec.origin.x + 20 * spec.spacing.x == -60.0
-    assert spec.reaction.source_rates[outside] == 0.0
-    assert spec.reaction.loss_rates[outside] == 5.0
-    assert spec.reaction.source_rates[sites + outside] == 0.0
-    assert spec.reaction.loss_rates[sites + outside] == 0.5
-    assert spec.reaction.source_rates[inside] == 0.0
-    assert spec.reaction.loss_rates[inside] == 0.0
-    assert spec.reaction.source_rates[sites + inside] == 20.0
-    assert spec.reaction.loss_rates[sites + inside] == 2.0
+    assert spec.reaction is None
+    assert spec.velocity_field is not None
+    assert any(value != 0.0 for value in spec.velocity_field.y_faces)
+    # The solved field is dominated by the axial channel flow; transverse
+    # components exist only as weak circulation at the trap mouth.
+    assert max(abs(value) for value in spec.velocity_field.x_faces) < max(
+        abs(value) for value in spec.velocity_field.y_faces
+    )
+    solid = sum(spec.obstacles)
+    assert 0 < solid < len(spec.obstacles)
+    assert spec.y_lower.values == [0.0, 10.0]
+    assert len(checkpoint.constraints.boxes) == 4
+
+
+def test_pillar_channel_anchors_sheds_and_washes_out() -> None:
+    model, _ = build_model(
+        _TUTORIALS / "pillar_channel.py",
+        ModelContext(BackendKind.CPU, 0, seed=7),
+    )
+    assert isinstance(model, SimulationController)
+    # 250 steps crosses the Brinkman re-solve cadence at step 100 and sheds
+    # daughters from every anchored lineage into the stream.
+    for _ in range(250):
+        model.step(0.02)
+
+    adhesion_sites = ((-20.0, -46.0), (20.0, -46.0), (0.0, 14.0))
+    cells = model.simulation.cells()
+    anchored = [cell for cell in cells if cell.fixed]
+    released = [cell for cell in cells if not cell.fixed]
+    assert len(anchored) == 3
+    assert len(released) > 3
+    for cell in anchored:
+        nearest = min(
+            math.hypot(cell.position.x - x, cell.position.y - y) for x, y in adhesion_sites
+        )
+        assert nearest < 4.0
+    # Released cells drift downstream of the anchors; the flow is doing work.
+    assert any(cell.position.y > 30.0 for cell in released)
+    for cell in cells:
+        assert cell.position.z == 0.0
+        assert abs(cell.position.x) < 40.0
+        assert abs(cell.position.y) < 120.0
 
 
 def test_plasmid_tutorial_resume_is_exact(tmp_path: Path) -> None:
diff --git a/scripts/run_nutrient_benchmarks.py b/scripts/run_nutrient_benchmarks.py
new file mode 100644
index 0000000..af0ddf7
--- /dev/null
+++ b/scripts/run_nutrient_benchmarks.py
@@ -0,0 +1,238 @@
+#!/usr/bin/env python3
+"""Quantify nutrient penetration, attached growth, and numerical sensitivity.
+
+A fixed population in a perfused 20 x 40 x 4 channel uses a physical averaging
+kernel for both nutrient uptake and resistance. Cells remain attached; this
+controlled experiment excludes division, mechanics, and detachment. All units
+are micrometers and model time, with nutrient in an arbitrary amount/volume
+scale. See docs/tutorials/nutrient-validation.md for equations and scope.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import math
+import subprocess
+from dataclasses import asdict, dataclass
+from pathlib import Path
+
+import numpy as np
+from microsimulator import (
+    BackendKind,
+    CellInit,
+    SignalGridAffineReaction,
+    SignalIntegrationKind,
+    Simulation,
+    Vec3,
+)
+from microsimulator.biomass import biomass_volume
+from microsimulator.flow import colony_mobility, colony_volume_fraction, solve_flow_field
+from microsimulator.flow_reference import duct_grid
+
+WIDTH, LENGTH, HEIGHT = 20.0, 40.0, 4.0
+DIFFUSION, SPEED = 4.0, 0.3
+MU, MONOD_K, YIELD = 0.2, 0.2, 0.5
+RADIUS, AVERAGING_RADIUS, DRAG = 0.5, 4.0, 40.0
+
+
+@dataclass
+class Result:
+    spacing: float
+    dt: float
+    refresh: float
+    duration: float
+    biomass_gain: float
+    penetration_half: float
+    upstream_growth: float
+    downstream_growth: float
+    nutrient_amount: float
+    boundary_supply: float
+    balance_relative_error: float
+    y: list[float]
+    nutrient_profile: list[float]
+    growth_y: list[float]
+    growth_rate: list[float]
+
+
+def run(h: float, dt: float, refresh: float, duration: float, backend: BackendKind) -> Result:
+    nx, ny = round(WIDTH / h), round(LENGTH / h)
+    spec = duct_grid(nx, ny, 1, (h, h, HEIGHT))
+    spec.origin = Vec3(h / 2, h / 2, HEIGHT / 2)
+    spec.diffusion, spec.integration = [DIFFUSION], SignalIntegrationKind.BACKWARD_EULER
+    spec.y_lower.values, spec.y_upper.values = [1], [0]
+    spec.solver.absolute_tolerance = 1e-7
+    spec.solver.relative_tolerance = 1e-6
+    sim = Simulation(backend)
+    field, _ = solve_flow_field(spec, mean_inlet_speed=SPEED, simulation=sim)
+    spec.velocity_field = field
+    sim.configure_signal_grid(spec, [0] * spec.site_count)
+    ids: list[int] = []
+    for x in (3, 6.5, 10, 13.5, 17):
+        for y in (3, 7, 11, 15, 19, 23, 27, 31, 35):
+            cell = CellInit()
+            cell.position = Vec3(x, y, HEIGHT / 2)
+            cell.length, cell.radius, cell.growth_rate, cell.fixed = 2, RADIUS, 0, True
+            ids.append(sim.add_cell(cell))
+    cells = sim.cells()
+    initial = np.array([biomass_volume(c.length, c.radius) for c in cells], dtype=np.float64)
+    # K_i integrates to one. The exact voxel-integrated physical kernel is
+    # independent of rod growth and spacing; only its attached amount changes.
+    kernels = np.stack(
+        [
+            colony_volume_fraction(spec, [c], averaging_radius=AVERAGING_RADIUS).ravel() / amount
+            for c, amount in zip(cells, initial, strict=True)
+        ]
+    )
+    volume = spec.voxel_volume
+    weights = kernels * volume
+    assert np.max(np.abs(weights.sum(axis=1) - 1)) < 1e-6
+    initial_y = [float(c.position.y) for c in cells]
+    boundary_supply = 0.0
+    next_refresh = 0.0
+    steps = round(duration / dt)
+    if not math.isclose(steps * dt, duration, abs_tol=1e-7):
+        raise ValueError("duration must be an integer number of steps")
+    for step in range(steps):
+        time = step * dt
+        cells = sim.cells()
+        if time + 1e-9 >= next_refresh:
+            mobility = colony_mobility(
+                spec, cells, drag_coefficient=DRAG, averaging_radius=AVERAGING_RADIUS
+            )
+            field, _ = solve_flow_field(
+                spec, mean_inlet_speed=SPEED, mobility=mobility, simulation=sim
+            )
+            sim.set_velocity_field(field)
+            next_refresh += refresh
+        old = np.asarray(sim.signal_levels, dtype=np.float64)
+        mean = weights @ old
+        cylinder = np.array([math.pi * c.radius**2 * c.length for c in cells])
+        coefficient = MU * cylinder / (YIELD * (MONOD_K + mean))
+        reaction = SignalGridAffineReaction()
+        reaction.source_rates = [0] * spec.site_count
+        reaction.loss_rates = (coefficient @ kernels).tolist()
+        sim.set_signal_reaction(reaction)
+        sim.step(dt)
+        new = np.asarray(sim.signal_levels, dtype=np.float64)
+        consumed = dt * coefficient * (weights @ new)
+        for cid, cell, amount in zip(ids, cells, consumed, strict=True):
+            length = cell.length + YIELD * float(amount) / (math.pi * cell.radius**2)
+            sim.set_cell_geometry(cid, cell.position, cell.direction, length)
+        # Exact discrete BE boundary flux convention used by engine transport:
+        # boundary values are ghost-center concentrations, one h from a site.
+        concentration = new.reshape(nx, ny)
+        y_faces = np.asarray(field.y_faces, dtype=np.float64).reshape(nx, ny + 1)
+        influx = DIFFUSION / h * (1 - concentration[:, 0]) + y_faces[:, 0]
+        outflux = DIFFUSION / h * concentration[:, -1] + y_faces[:, -1] * concentration[:, -1]
+        boundary_supply += dt * h * HEIGHT * float((influx - outflux).sum())
+    final_cells = sim.cells()
+    final = np.array([biomass_volume(c.length, c.radius) for c in final_cells])
+    nutrient = np.asarray(sim.signal_levels, dtype=np.float64).reshape(nx, ny)
+    profile = nutrient.mean(axis=0)
+    y = (np.arange(ny) + 0.5) * h
+    # Boundary sample follows the ghost-center discretization above.
+    # First downstream crossing, with no assumption of monotonicity farther on.
+    previous_y, previous_c = -h / 2, 1.0
+    penetration = LENGTH
+    for yy, cc in zip(y, profile, strict=True):
+        coordinate, concentration = float(yy), float(cc)
+        if concentration <= 0.5:
+            penetration = previous_y + (coordinate - previous_y) * (previous_c - 0.5) / (
+                previous_c - concentration
+            )
+            break
+        previous_y, previous_c = coordinate, concentration
+    growth = (final - initial) / (duration * initial)
+    gain = float((final - initial).sum())
+    nutrient_amount = float(nutrient.sum()) * volume
+    balance = abs(nutrient_amount + gain / YIELD - boundary_supply) / max(boundary_supply, 1e-12)
+    return Result(
+        h,
+        dt,
+        refresh,
+        duration,
+        gain,
+        penetration,
+        float(growth[np.array(initial_y) <= 11].mean()),
+        float(growth[np.array(initial_y) >= 27].mean()),
+        nutrient_amount,
+        boundary_supply,
+        balance,
+        y.tolist(),
+        profile.tolist(),
+        initial_y,
+        growth.tolist(),
+    )
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--output", type=Path, required=True)
+    parser.add_argument("--backend", choices=("cpu", "metal", "cuda"), default="cpu")
+    parser.add_argument("--duration", type=float, default=8)
+    args = parser.parse_args()
+    backend = {"cpu": BackendKind.CPU, "metal": BackendKind.METAL, "cuda": BackendKind.CUDA}[
+        args.backend
+    ]
+    cases = [
+        (2.0, 0.04, 0.4),
+        (1.0, 0.04, 0.4),
+        (0.5, 0.04, 0.4),
+        (1.0, 0.02, 0.4),
+        (1.0, 0.04, 0.2),
+    ]
+    results: list[Result] = []
+    for case in cases:
+        result = run(*case, args.duration, backend)
+        print(
+            f"h={result.spacing:g} dt={result.dt:g} refresh={result.refresh:g}: "
+            f"gain={result.biomass_gain:.6g}, penetration={result.penetration_half:.6g}, "
+            f"balance error={result.balance_relative_error:.3g}",
+            flush=True,
+        )
+        results.append(result)
+    base = results[1]
+    comparisons: dict[str, dict[str, float]] = {}
+    for label, other in [
+        ("grid_2_to_1", results[0]),
+        ("grid_1_to_half", results[2]),
+        ("dt_halved", results[3]),
+        ("refresh_halved", results[4]),
+    ]:
+        comparisons[label] = {
+            metric: abs(getattr(other, metric) - getattr(base, metric))
+            / max(abs(getattr(base, metric)), 1e-12)
+            for metric in ("biomass_gain", "penetration_half", "upstream_growth")
+        }
+    checks: dict[str, bool] = {
+        "mass_balance": max(r.balance_relative_error for r in results) < 2e-4,
+        "spatial_growth": base.upstream_growth > 3 * base.downstream_growth,
+        "grid_gain": comparisons["grid_1_to_half"]["biomass_gain"] < 0.1,
+        "time_gain": comparisons["dt_halved"]["biomass_gain"] < 0.02,
+        "refresh_gain": comparisons["refresh_halved"]["biomass_gain"] < 0.02,
+    }
+    payload = {
+        "source_commit": subprocess.check_output(["git", "rev-parse", "HEAD"], text=True).strip(),
+        "backend": args.backend,
+        "parameters": {
+            "diffusion": DIFFUSION,
+            "speed": SPEED,
+            "growth_rate": MU,
+            "monod_k": MONOD_K,
+            "yield": YIELD,
+            "averaging_radius": AVERAGING_RADIUS,
+            "drag": DRAG,
+        },
+        "runs": [asdict(r) for r in results],
+        "relative_changes": comparisons,
+        "checks": checks,
+    }
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    args.output.write_text(json.dumps(payload, indent=2) + "\n")
+    if not all(checks.values()):
+        raise SystemExit(f"Nutrient validation failed: {checks}")
+
+
+if __name__ == "__main__":
+    main()