From b202446452dd0cdd3ad864d14a33a50234df5ad9 Mon Sep 17 00:00:00 2001 From: Mike Arpaia Date: Thu, 20 Aug 2026 11:01:49 -0600 Subject: [PATCH 1/2] Add microfluidic device and biology tutorials --- docs/tutorials/README.md | 2 + docs/tutorials/flow-solvers.md | 166 +++++++++++++++ docs/tutorials/microfluidics.md | 237 +++++++++++++++++++++ docs/tutorials/simbol.md | 60 ++++-- examples/microfluidic_trap.py | 168 +++++++++++++++ examples/tutorials/biopixel_trap.py | 195 +++++++++++++++++ examples/tutorials/danino_clock.py | 276 +++++++++++++------------ examples/tutorials/pillar_channel.py | 252 ++++++++++++++++++++++ python/tests/test_masks.py | 65 ++++++ python/tests/test_microfluidics.py | 76 ++++++- python/tests/test_signal_model_runs.py | 62 ++++++ python/tests/test_tutorials.py | 60 ++++-- 12 files changed, 1452 insertions(+), 167 deletions(-) create mode 100644 docs/tutorials/flow-solvers.md create mode 100644 docs/tutorials/microfluidics.md create mode 100644 examples/microfluidic_trap.py create mode 100644 examples/tutorials/biopixel_trap.py create mode 100644 examples/tutorials/pillar_channel.py create mode 100644 python/tests/test_signal_model_runs.py diff --git a/docs/tutorials/README.md b/docs/tutorials/README.md index 6a21867..27952eb 100644 --- a/docs/tutorials/README.md +++ b/docs/tutorials/README.md @@ -11,6 +11,8 @@ These tutorials introduce the CellModeller2 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 `cm` command, data-only checkpoints, and the standalone viewer. Each model selects its backend explicitly and can be run headlessly for batch experiments. diff --git a/docs/tutorials/flow-solvers.md b/docs/tutorials/flow-solvers.md new file mode 100644 index 0000000..3a7b5bd --- /dev/null +++ b/docs/tutorials/flow-solvers.md @@ -0,0 +1,166 @@ +# Solved flow: a pillar channel, Brinkman feedback, and the benchmarks + +The previous tutorial builds devices from the packaged helpers, whose grids and flow come +preassembled. This one goes a level down and demonstrates the flow machinery itself on +geometry no helper covers and no formula describes: a monolayer channel crossed by a +staggered array of cylindrical pillars, with colonies adhered in the pillar wakes shedding +daughters into the stream. Everything lives in one model: + +```console +uv run cm view --model examples/tutorials/pillar_channel.py --seed 7 --dt 0.01 --backend metal --open +``` + +[`examples/tutorials/pillar_channel.py`](../../examples/tutorials/pillar_channel.py) +demonstrates, in order: authoring walls and a solid mask directly, solving the flow through +an arbitrary geometry, adhesion with released daughters, Brinkman colony feedback on the +flow, and washout. The end state is a steady flow of cells: anchored lineages grow in the +wakes while the stream continuously carries their offspring between the pillars and out of +the channel. + +## Geometry the solver has to earn + +The channel is one inside-region box; each pillar is one outside-region z-cylinder: + +```python +pillar = CylinderConstraintInit() +pillar.center = Vec3(x, y, 0.0) +pillar.radius = PILLAR_RADIUS +pillar.allowed_region = ConstraintRegion.OUTSIDE +simulation.add_cylinder_constraint(pillar) +``` + +The signal grid needs the same geometry as a voxel mask, and curved walls raise a rule that +axis-aligned devices never surface: **the mechanics geometry must enclose the solid mask**. +Mechanics keeps cell centers outside the smooth cylinder; the mask is stair-stepped. If a +voxel whose center is barely inside the circle is marked solid, its corners poke out past +the cylinder wall, a cell hugging the wall can stand inside a solid voxel, and sampling +signals at its center is an error. The model therefore voxelizes conservatively — a voxel +is solid only when it lies *entirely* inside the pillar: + +```python +core = PILLAR_RADIUS - 0.5 * math.hypot(spacing.x, spacing.y) +solid = (px - x) ** 2 + (py - y) ** 2 < core * core +``` + +With that rule every reachable cell position is in fluid, and the stair-stepped flow +blockage errs on the small side by the same half-diagonal margin. + +## Solving flow where no profile exists + +An analytic profile for a pillar array does not exist; the field comes from the numerical +solve, exactly as in the device helpers: + +```python +field, report = solve_flow_field(grid, mean_inlet_speed=FLOW_SPEED, mobility=gap_mobility(grid)) +grid.velocity_field = field +``` + +The solved field is conservative per voxel and routes around every pillar. At a mean inlet +speed of 20 the plug away from the array runs at 20 as requested — the solve normalizes over +the open inlet faces, so blocked columns cannot inflate it — and the gaps beside the center +pillar carry ≈31, because the pillars take cross-section and the same flux has to fit +through what is left. Flow speeds up exactly where the physical device would. + +`report.max_speed` gives the number the `dt` bound needs. Drift is an explicit step, so a +cell must not cross more than about its own radius per step: keep `max_speed * dt` below +`CELL_RADIUS`. Here `max_speed` is 37.7, so `--dt 0.01` leaves a comfortable margin and +`--dt 0.02` would exceed it. + +## Adhesion: anchored mothers, shed daughters + +Founders are placed in the pillar wakes with `fixed = True` — mechanics and drift never +move them. Daughters inherit adhesion, so each division decides who stays: + +```python +def _divided(step, event): + DIVISION.on_division(step, event) + if event.parent.fixed: + released = ... # the daughter farther from the adhesion site + step.simulation.set_cell_fixed(released.id, False) +``` + +The daughter nearer the adhesion site keeps the anchor; the other is released and the +stream takes it. Anchoring by *site* rather than by daughter order matters: division +displaces both daughters by half a cell length, and a lineage that anchors whichever +daughter comes first random-walks away from its wake — and since fixed cells are never +pushed by mechanics, a walking anchor eventually stands inside a pillar. Site anchoring +keeps each attached cell within about a cell length of where its founder adhered, +indefinitely. + +Released cells drift with the local fluid velocity (`MechanicsConfig(flow_drift=True)`), +slowly in the wake, then fast in the gaps, and leave through plan removals at the channel +end — the same washout pattern as the trap models. Run at seed 7 for 1600 steps at +`dt = 0.01` and the population reaches a steady state: three anchored cells, on the order +of 140 in transit, and about 1200 washed out, with lineage recording every one. + +## The colony pushes back on the flow + +Like the trap models, the regulation step re-solves the flow at a fixed cadence with the +colony rasterized into Brinkman drag and swaps the field into the running simulation: + +```python +if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0: + mobility = colony_mobility(GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + step.simulation.set_velocity_field(field) +``` + +Here the feedback has a visible consequence: as a wake colony thickens it plugs its own +gap, the solve routes more of the flux through the neighboring gaps, and shed daughters +increasingly take the fast lanes around the crowd. + +## Checking the closure with the resolved solver + +The Hele-Shaw field the model runs on is a depth-averaged closure. The MAC Stokes–Brinkman +solver resolves the same problem with viscous boundary layers on every wall, and takes the +identical grid: + +```python +from cellmodeller2.stokes import solve_stokes_field + +resolved, report = solve_stokes_field(GRID, mean_inlet_speed=FLOW_SPEED) +``` + +Depth-averaging the resolved field reproduces the closure's flux split around obstacles to +under a percent in the thin-gap regime — that agreement is enforced continuously by the +benchmark suite, which validates both solvers against literature and exact references +(plane Poiseuille and the two-layer Brinkman channel against their exact solutions with +measured second-order convergence, the Shah–London square-duct peak-to-mean ratio 2.0962, +and the cross-solver thin-gap check): + +```console +uv run python scripts/run_flow_benchmarks.py # CI-gating benchmark table +uv run python scripts/run_flow_benchmarks.py --fine # doubled resolutions +``` + +### Which solver a grid deserves + +The two solvers are accurate in opposite regimes, and the grid decides which. + +The MAC solve resolves a no-slip profile only where it has voxels to resolve it in. Every +solve reports `min_gap_voxels`, the fluid voxels across its narrowest channel. One voxel +carries roughly two and a half times the flux the true parabolic profile would; four +voxels bring that within about ten percent and eight within a few percent. The pillar grid +here is one voxel deep in z, so its `min_gap_voxels` is 1 and the *depth-averaged* MAC +answer is the meaningful one — a comparison of in-plane flux splits, not of absolute +speeds. + +The Hele-Shaw closure has the mirror-image property. It carries the gap-height physics +analytically in its mobility, so it is accurate for shallow channels at any z resolution, +but it solves for the depth-averaged velocity: every z layer of a column gets the column's +mean. In a channel resolved across its depth, the cells near the floor drift at the mean +rather than at the slower speed the true profile gives them, and a rod sees no shear +across the gap. + +So: use the closure for device authoring and the in-model re-solve cadence, and reach for +`solve_stokes_field` when a study needs resolved wall shear or true cross-channel profiles +*and* the grid resolves the gap with at least four voxels. Refining z to reach that costs +grid sites in every other subsystem too. + +## What to watch in the viewer + +- Streamwise plug flow entering the array, cells accelerating visibly through the gaps. +- Wake colonies growing in place while a comet tail of released cells stretches + downstream from each pillar. +- Cells vanishing at the channel end as washout removes them; lineage keeps their + ancestry for analysis. diff --git a/docs/tutorials/microfluidics.md b/docs/tutorials/microfluidics.md new file mode 100644 index 0000000..565fe13 --- /dev/null +++ b/docs/tutorials/microfluidics.md @@ -0,0 +1,237 @@ +# 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 cm 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 + +`cellmodeller2.microfluidics.TrapChannelDevice` describes an open-sided trap fed by a +straight channel and projects that one description into every engine input: + +```python +from cellmodeller2.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]) +``` + +`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). 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: arbitrary geometry and colony feedback + +Device flow fields are solved, not authored: `cellmodeller2.flow` computes the steady +Hele-Shaw–Brinkman problem over the grid's fluid voxels and returns the same face-staggered +field the engine consumes. `apply_to_grid` runs this solve for every device, and it works +for any mask geometry — junctions, bends, pillars, a CAD-derived layout — not just straight +channels. The solver is also available directly for grids built without a device helper: + +```python +from cellmodeller2.flow import colony_mobility, solve_flow_field + +field, report = solve_flow_field(grid, mean_inlet_speed=20.0) # Stokes limit +grid.velocity_field = field +``` + +The solve is a variable-coefficient pressure problem (`div(m grad p) = 0`), so the returned +fluxes conserve mass per voxel and vanish on wall faces by construction; the flow-axis +boundaries must be `FIXED` to act as inlet and outlet, and the linear solution is rescaled to +the requested mean inlet speed. With uniform mobility this is the Stokes limit of the +depth-averaged closure — correct routing through any mask, plug profile across the channel +width (side-wall boundary layers, of order the gap height, are outside the closure). + +The mobility field is where Brinkman feedback enters: `colony_mobility` rasterizes the +colony's volume fraction and adds Kozeny–Carman style drag, so media diverts around a packed +trap and seeps through its edges. Because the field is data, regulation code re-solves as +the colony grows and swaps it into the running simulation — the trap models do this every +`RESOLVE_INTERVAL` steps: + +```python +def _regulate(step: ControllerStep) -> StepPlan: + if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0: + mobility = colony_mobility(GRID, step.cells, drag_coefficient=DRAG_COEFFICIENT) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + step.simulation.set_velocity_field(field) + ... +``` + +`Simulation.set_velocity_field` validates the replacement against the full grid +specification before swapping it; transport, drift, and checkpoints all use whichever field +is current. The re-solve cadence is a model choice, and the trade is cost against staleness: +the trap models re-solve every hundred steps, which at `dt = 0.02` is a couple of doublings +of colony growth, so the field the drift and transport see lags the colony by that much. +Shorten the interval where the blockage matters quantitatively. The drag coefficient is a +modeling parameter (how strongly a packed colony resists through-flow relative to the open +channel), not a measured constant. + +### Resolved flow: the MAC Stokes–Brinkman solver + +When a study needs the flow the closure cannot express — viscous boundary layers on side +walls, the true cross-channel profile, resolved wall shear — `cellmodeller2.stokes` solves +the full staggered-grid Stokes–Brinkman problem with the same call shape and returns the +same engine-ready field: + +```python +from cellmodeller2.stokes import colony_drag, solve_stokes_field + +field, report = solve_stokes_field(grid, mean_inlet_speed=20.0) +field, report = solve_stokes_field( + grid, mean_inlet_speed=20.0, drag=colony_drag(grid, cells, drag_coefficient=0.4) +) +``` + +It costs far more than the Hele-Shaw solve, so devices keep the closure for authoring and +in-loop feedback; the MAC solver anchors it where the grid resolves the gap. Each solve +reports `min_gap_voxels`, the fluid voxels across its narrowest channel: below about four +the MAC solve over-predicts that channel's flux and the depth-averaged closure is the more +accurate model, which is why the shallow device grids here stay on the closure. Both +solvers run against literature and exact references in `scripts/run_flow_benchmarks.py` — +plane Poiseuille and the two-layer Brinkman channel against their exact solutions with +measured second-order convergence, the Shah–London square-duct peak-to-mean ratio, and a +thin-gap cross-check in which the depth-averaged MAC solution reproduces the Hele-Shaw flux +split around a pillar: + +```console +uv run python scripts/run_flow_benchmarks.py # CI-gating benchmark table +uv run python scripts/run_flow_benchmarks.py --fine # doubled resolutions +``` + +The next tutorial, [Solved flow](flow-solvers.md), exercises all of this machinery on a +pillar-array channel built without any device helper. + +## 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 cellmodeller2.microfluidics import BiopixelTrapDevice + +DEVICE = BiopixelTrapDevice(mean_flow_speed=20.0) +``` + +### Reading the supplied CAD layout + +`cellmodeller2.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 cellmodeller2.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 cm 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` gives a doubling time of `ln 2 ≈ 0.69` model time units. 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. diff --git a/docs/tutorials/simbol.md b/docs/tutorials/simbol.md index f5fad4c..9e32a2c 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 cm 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 `cellmodeller2.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. CellModeller2 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 CellModeller2 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..dda4916 --- /dev/null +++ b/examples/microfluidic_trap.py @@ -0,0 +1,168 @@ +"""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. 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. 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 cellmodeller2 import ( + CellInit, + CellUpdate, + ControllerStep, + CoupledRatePlan, + DivisionEvent, + GridShape, + MechanicsConfig, + ModelContext, + NativeController, + RatePlanBuilder, + SignalGridSpec, + SignalIntegrationKind, + Simulation, + StepPlan, + UniformLengthDivision, + Vec3, +) +from cellmodeller2.checkpoint import CheckpointBundle, JSONValue +from cellmodeller2.flow import colony_mobility, gap_mobility, solve_flow_field +from cellmodeller2.microfluidics import TrapChannelDevice + +MODEL_ID = "examples.microfluidic-trap" +MODEL_VERSION = 3 +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 is one limiting substrate in arbitrary concentration units, fed at +# NUTRIENT_INLET. Uptake is tied to realized growth: a cell consumes +# growth_rate * volume / NUTRIENT_YIELD per unit time, so Monod-limited growth +# and consumption stay consistent. The yield sets the coupling strength, and +# this value makes a packed trap's uptake comparable to the diffusive supply +# through its mouth, so nutrient penetrates a few tens of micrometers and the +# colony behind that front grows more slowly. +NUTRIENT_YIELD = 0.5 +WASHOUT_Y = DEVICE.channel_half_length - 10.0 + +# Brinkman feedback: how often the colony's drag re-solves the device flow, +# and how strongly a packed voxel resists through-flow. +RESOLVE_INTERVAL = 100 +DRAG_COEFFICIENT = 100.0 + + +def _grid() -> 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.CRANK_NICOLSON + DEVICE.apply_to_grid(grid, inlet_values=[NUTRIENT_INLET], outlet_values=[0.0]) + return grid + + +GRID = _grid() +GAP_MOBILITY = gap_mobility(GRID) + + +def _rate_plan() -> CoupledRatePlan: + rates = RatePlanBuilder() + uptake = -(rates.growth_rate() * rates.cell_volume()) / 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, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT + ) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) + 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..785a7a0 --- /dev/null +++ b/examples/tutorials/biopixel_trap.py @@ -0,0 +1,195 @@ +"""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 pathlib import Path + +from cellmodeller2 import ( + CellInit, + CellUpdate, + ControllerStep, + CoupledRatePlan, + DivisionEvent, + GridShape, + MechanicsConfig, + ModelContext, + NativeController, + RatePlanBuilder, + SignalGridSpec, + SignalIntegrationKind, + Simulation, + StepPlan, + UniformLengthDivision, + Vec3, +) +from cellmodeller2.checkpoint import CheckpointBundle, JSONValue +from cellmodeller2.flow import colony_mobility, gap_mobility, solve_flow_field +from cellmodeller2.masks import ( + MaskError, + MaskRectangle, + extract_rectangles, + load_mask_polylines, + match_rectangles, +) +from cellmodeller2.microfluidics import BiopixelTrapDevice + +MODEL_ID = "tutorials.biopixel-trap" +MODEL_VERSION = 6 +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 is one limiting substrate in arbitrary concentration units, fed at +# NUTRIENT_INLET. Uptake is tied to realized growth: a cell consumes +# growth_rate * volume / NUTRIENT_YIELD per unit time, so Monod-limited growth +# and consumption stay consistent. The yield sets the coupling strength, and +# this value makes a packed trap's uptake comparable to the diffusive supply +# through its mouth, so nutrient penetrates a few tens of micrometers and the +# colony behind that front grows more slowly. +NUTRIENT_YIELD = 0.5 + +# Brinkman feedback: how often the colony's drag re-solves the device flow, +# and how strongly a packed voxel resists through-flow. +RESOLVE_INTERVAL = 100 +DRAG_COEFFICIENT = 100.0 + + +def _grid() -> 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.CRANK_NICOLSON + DEVICE.apply_to_grid(grid, inlet_values=[NUTRIENT_INLET], outlet_values=[0.0]) + return grid + + +GRID = _grid() +GAP_MOBILITY = gap_mobility(GRID) + + +def _rate_plan() -> CoupledRatePlan: + rates = RatePlanBuilder() + uptake = -(rates.growth_rate() * rates.cell_volume()) / 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, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT + ) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) + 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 7394c99..fb02fe6 100644 --- a/examples/tutorials/danino_clock.py +++ b/examples/tutorials/danino_clock.py @@ -1,11 +1,20 @@ -"""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 +import numpy as np from cellmodeller2 import ( CellInit, + CellSnapshot, CellUpdate, ControllerStep, CoupledRatePlan, @@ -14,74 +23,122 @@ MechanicsConfig, ModelContext, NativeController, - PlaneConstraintInit, RatePlanBuilder, SignalGridAffineReaction, SignalGridSpec, SignalIntegrationKind, Simulation, - SphereConstraintInit, - SphereRegion, StepPlan, UniformLengthDivision, Vec3, ) from cellmodeller2.checkpoint import CheckpointBundle, JSONValue +from cellmodeller2.flow import ( + colony_mobility, + colony_species_density, + gap_mobility, + solve_flow_field, +) +from cellmodeller2.microfluidics import TrapChannelDevice MODEL_ID = "tutorials.danino-clock" -MODEL_VERSION = 2 +MODEL_VERSION = 7 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 is one limiting substrate in arbitrary concentration units, fed at +# NUTRIENT_INLET. Uptake is tied to realized growth: a cell consumes +# growth_rate * volume / NUTRIENT_YIELD per unit time, so Monod-limited growth +# and consumption stay consistent. The yield sets the coupling strength, and +# this value makes a packed trap's uptake comparable to the diffusive supply +# through its mouth, so nutrient penetrates a few tens of micrometers and the +# colony behind that front grows more slowly. +NUTRIENT_YIELD = 0.5 + +# Brinkman feedback: how often the colony's drag re-solves the device flow, +# and how strongly a packed voxel resists through-flow. +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: 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 + DEVICE.apply_to_grid( + grid, + inlet_values=[0.0, NUTRIENT_INLET], + outlet_values=[0.0, 0.0], + ) 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 +146,44 @@ 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.growth_rate() * rates.cell_volume()) / 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. + Summed over the cells of a voxel and divided by its volume, that is 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,27 @@ 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, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT + ) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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 +227,12 @@ 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()) + 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 +247,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..a9ee084 --- /dev/null +++ b/examples/tutorials/pillar_channel.py @@ -0,0 +1,252 @@ +"""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 the colony's Brinkman drag so growing colonies divert the flow. 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 cellmodeller2 import ( + BoxConstraintInit, + CellInit, + CellUpdate, + ConstraintRegion, + ControllerStep, + CoupledRatePlan, + CylinderConstraintInit, + DivisionEvent, + GridBoundaryKind, + GridShape, + MechanicsConfig, + ModelContext, + NativeController, + RatePlanBuilder, + SignalGridSpec, + SignalIntegrationKind, + Simulation, + StepPlan, + UniformLengthDivision, + Vec3, +) +from cellmodeller2.checkpoint import CheckpointBundle, JSONValue +from cellmodeller2.flow import colony_mobility, gap_mobility, solve_flow_field + +MODEL_ID = "tutorials.pillar-channel" +MODEL_VERSION = 2 +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 is one limiting substrate in arbitrary concentration units, fed at +# NUTRIENT_INLET. Uptake is tied to realized growth: a cell consumes +# growth_rate * volume / NUTRIENT_YIELD per unit time, so Monod-limited growth +# and consumption stay consistent. The yield sets the coupling strength, and +# this value makes a packed trap's uptake comparable to the diffusive supply +# through its mouth, so nutrient penetrates a few tens of micrometers and the +# colony behind that front grows more slowly. +NUTRIENT_YIELD = 0.5 + +# Brinkman feedback: how often the colony's drag re-solves the device flow, +# and how strongly a packed voxel resists through-flow. +RESOLVE_INTERVAL = 100 +DRAG_COEFFICIENT = 100.0 + + +def _in_pillar_core(px: float, py: float, margin: float) -> bool: + # A voxel is solid only when it lies entirely inside the pillar (its + # center plus half the voxel diagonal stays within the radius). The + # mechanics cylinders therefore enclose every solid voxel, so a cell + # center can never sit inside the mask and signal sampling is always in + # fluid; the stair-stepped flow blockage is conservative by the same + # margin. + core = PILLAR_RADIUS - margin + if core <= 0.0: + return False + return any( + (px - x) * (px - x) + (py - y) * (py - y) < core * core for x, y in PILLARS + ) + + +def _grid() -> 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, -8.0) + grid.spacing = Vec3(4.0, 4.0, 4.0) + grid.diffusion = [40.0] + grid.advection = [Vec3()] + grid.integration = SignalIntegrationKind.CRANK_NICOLSON + margin = 0.5 * math.hypot(grid.spacing.x, grid.spacing.y) + 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, margin) + ) + 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) + field, _ = solve_flow_field( + grid, mean_inlet_speed=FLOW_SPEED, mobility=gap_mobility(grid) + ) + 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.growth_rate() * rates.cell_volume()) / 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, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT + ) + field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) + simulation.configure_signal_grid(GRID, _primed_levels(GRID)) + simulation.set_coupled_rate_plan(_rate_plan()) + _add_walls(simulation) + + founder_ids = [] + 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 a0c85e7..cdf600e 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 3cb6a23..fe56b94 100644 --- a/python/tests/test_microfluidics.py +++ b/python/tests/test_microfluidics.py @@ -3,9 +3,21 @@ from __future__ import annotations import math - -from cellmodeller2 import GridShape, SignalGridSpec, Vec3 +from pathlib import Path + +import pytest +from cellmodeller2 import ( + BackendKind, + GridShape, + ModelContext, + SignalGridSpec, + SimulationController, + Vec3, +) from cellmodeller2.microfluidics import BiopixelTrapDevice, TrapChannelDevice +from cellmodeller2.runner import build_model + +_EXAMPLES = Path(__file__).resolve().parents[2] / "examples" def _grid() -> SignalGridSpec: @@ -116,3 +128,63 @@ def test_biopixel_cavity_is_shallow_beside_the_model_channel() -> None: assert device._solid(42.5, 0.0, 2.475, half) assert not device._solid(-50.0, 0.0, 9.075, half) assert math.isclose(device.trap_height / device.channel_height, 0.165) + + +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 + + +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..d674cf4 --- /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 cellmodeller2 import ( + BackendKind, + ModelContext, + SimulationController, + build_model, +) +from cellmodeller2.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 68181fa..2b450e2 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: From 1ffe2d3a4d7a2ce4d33973feb63ba73144936638 Mon Sep 17 00:00:00 2001 From: Mike Arpaia Date: Sat, 29 Aug 2026 09:27:45 -0600 Subject: [PATCH 2/2] Route tutorial flow through selected backends --- docs/tutorials/flow-solvers.md | 31 ++++++++++++++------- docs/tutorials/microfluidics.md | 41 ++++++++++++++++++---------- examples/microfluidic_trap.py | 22 ++++++++++++--- examples/tutorials/biopixel_trap.py | 21 +++++++++++--- examples/tutorials/danino_clock.py | 17 +++++++++--- examples/tutorials/pillar_channel.py | 24 +++++++++++----- python/tests/test_microfluidics.py | 18 ++++++++++++ 7 files changed, 131 insertions(+), 43 deletions(-) diff --git a/docs/tutorials/flow-solvers.md b/docs/tutorials/flow-solvers.md index 3a7b5bd..dd53c17 100644 --- a/docs/tutorials/flow-solvers.md +++ b/docs/tutorials/flow-solvers.md @@ -51,15 +51,16 @@ An analytic profile for a pillar array does not exist; the field comes from the solve, exactly as in the device helpers: ```python -field, report = solve_flow_field(grid, mean_inlet_speed=FLOW_SPEED, mobility=gap_mobility(grid)) +field, report = solve_flow_field( + grid, + mean_inlet_speed=FLOW_SPEED, + mobility=gap_mobility(grid), + simulation=simulation, +) grid.velocity_field = field ``` -The solved field is conservative per voxel and routes around every pillar. At a mean inlet -speed of 20 the plug away from the array runs at 20 as requested — the solve normalizes over -the open inlet faces, so blocked columns cannot inflate it — and the gaps beside the center -pillar carry ≈31, because the pillars take cross-section and the same flux has to fit -through what is left. Flow speeds up exactly where the physical device would. +The solved field is conservative per voxel and routes around every pillar. Passing `simulation` selects its native CPU, Metal, or CUDA implementation; accelerator solves retain their pressure and Krylov vectors on the selected device. At a mean inlet speed of 20 the plug away from the array runs at 20 as requested; the solve normalizes over the open inlet faces, so blocked columns cannot inflate it. The gaps beside the center pillar carry approximately 31 because the pillars reduce the open cross-section and the same flux has to pass through what remains. `report.max_speed` gives the number the `dt` bound needs. Drift is an explicit step, so a cell must not cross more than about its own radius per step: keep `max_speed * dt` below @@ -101,7 +102,12 @@ colony rasterized into Brinkman drag and swaps the field into the running simula ```python if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0: mobility = colony_mobility(GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + field, _ = solve_flow_field( + GRID, + mean_inlet_speed=FLOW_SPEED, + mobility=mobility, + simulation=step.simulation, + ) step.simulation.set_velocity_field(field) ``` @@ -118,7 +124,11 @@ identical grid: ```python from cellmodeller2.stokes import solve_stokes_field -resolved, report = solve_stokes_field(GRID, mean_inlet_speed=FLOW_SPEED) +resolved, report = solve_stokes_field( + GRID, + mean_inlet_speed=FLOW_SPEED, + simulation=simulation, +) ``` Depth-averaging the resolved field reproduces the closure's flux split around obstacles to @@ -129,8 +139,9 @@ measured second-order convergence, the Shah–London square-duct peak-to-mean ra and the cross-solver thin-gap check): ```console -uv run python scripts/run_flow_benchmarks.py # CI-gating benchmark table -uv run python scripts/run_flow_benchmarks.py --fine # doubled resolutions +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 ``` ### Which solver a grid deserves diff --git a/docs/tutorials/microfluidics.md b/docs/tutorials/microfluidics.md index 565fe13..6b6c6e5 100644 --- a/docs/tutorials/microfluidics.md +++ b/docs/tutorials/microfluidics.md @@ -46,10 +46,15 @@ from cellmodeller2.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]) +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). 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. +`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 @@ -88,16 +93,15 @@ channels. The solver is also available directly for grids built without a device ```python from cellmodeller2.flow import colony_mobility, solve_flow_field -field, report = solve_flow_field(grid, mean_inlet_speed=20.0) # Stokes limit +field, report = solve_flow_field( + grid, + mean_inlet_speed=20.0, + simulation=simulation, +) # Stokes limit grid.velocity_field = field ``` -The solve is a variable-coefficient pressure problem (`div(m grad p) = 0`), so the returned -fluxes conserve mass per voxel and vanish on wall faces by construction; the flow-axis -boundaries must be `FIXED` to act as inlet and outlet, and the linear solution is rescaled to -the requested mean inlet speed. With uniform mobility this is the Stokes limit of the -depth-averaged closure — correct routing through any mask, plug profile across the channel -width (side-wall boundary layers, of order the gap height, are outside the closure). +The solve is a variable-coefficient pressure problem (`div(m grad p) = 0`), so the returned fluxes conserve mass per voxel and vanish on wall faces by construction; the flow-axis boundaries must be `FIXED` to act as inlet and outlet, and the linear solution is rescaled to the requested mean inlet speed. With uniform mobility this is the Stokes limit of the depth-averaged closure: correct routing through any mask and a plug profile across the channel width, with side-wall boundary layers outside the closure. The CPU implementation is C++, while Metal and CUDA execute independent MSL and CUDA kernels for the matrix-free operator and Krylov iterations; neither accelerator calls the CPU solver. The mobility field is where Brinkman feedback enters: `colony_mobility` rasterizes the colony's volume fraction and adds Kozeny–Carman style drag, so media diverts around a packed @@ -109,7 +113,12 @@ the colony grows and swaps it into the running simulation — the trap models do def _regulate(step: ControllerStep) -> StepPlan: if step.completed_steps and step.completed_steps % RESOLVE_INTERVAL == 0: mobility = colony_mobility(GRID, step.cells, drag_coefficient=DRAG_COEFFICIENT) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + field, _ = solve_flow_field( + GRID, + mean_inlet_speed=FLOW_SPEED, + mobility=mobility, + simulation=step.simulation, + ) step.simulation.set_velocity_field(field) ... ``` @@ -133,9 +142,12 @@ same engine-ready field: ```python from cellmodeller2.stokes import colony_drag, solve_stokes_field -field, report = solve_stokes_field(grid, mean_inlet_speed=20.0) +field, report = solve_stokes_field(grid, mean_inlet_speed=20.0, simulation=simulation) field, report = solve_stokes_field( - grid, mean_inlet_speed=20.0, drag=colony_drag(grid, cells, drag_coefficient=0.4) + grid, + mean_inlet_speed=20.0, + drag=colony_drag(grid, cells, drag_coefficient=0.4), + simulation=simulation, ) ``` @@ -151,8 +163,9 @@ thin-gap cross-check in which the depth-averaged MAC solution reproduces the Hel split around a pillar: ```console -uv run python scripts/run_flow_benchmarks.py # CI-gating benchmark table -uv run python scripts/run_flow_benchmarks.py --fine # doubled resolutions +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 next tutorial, [Solved flow](flow-solvers.md), exercises all of this machinery on a diff --git a/examples/microfluidic_trap.py b/examples/microfluidic_trap.py index dda4916..068acd8 100644 --- a/examples/microfluidic_trap.py +++ b/examples/microfluidic_trap.py @@ -13,6 +13,8 @@ from __future__ import annotations +from dataclasses import replace + from cellmodeller2 import ( CellInit, CellUpdate, @@ -61,7 +63,7 @@ DRAG_COEFFICIENT = 100.0 -def _grid() -> SignalGridSpec: +def _grid(simulation: Simulation | None = None) -> SignalGridSpec: shape = GridShape() shape.x, shape.y, shape.z = 64, 72, 6 grid = SignalGridSpec() @@ -77,7 +79,13 @@ def _grid() -> SignalGridSpec: grid.diffusion = [40.0] grid.advection = [Vec3()] grid.integration = SignalIntegrationKind.CRANK_NICOLSON - DEVICE.apply_to_grid(grid, inlet_values=[NUTRIENT_INLET], outlet_values=[0.0]) + 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 @@ -109,7 +117,12 @@ def _regulate(step: ControllerStep) -> StepPlan: mobility = colony_mobility( GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT ) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) @@ -133,7 +146,8 @@ def _divided(step: ControllerStep, event: DivisionEvent) -> None: def build(context: ModelContext) -> NativeController: simulation = context.simulation(reserved_capacity=10_000) - simulation.configure_signal_grid(GRID, _primed_levels(GRID)) + grid = _grid(simulation) + simulation.configure_signal_grid(grid, _primed_levels(grid)) simulation.set_coupled_rate_plan(_rate_plan()) DEVICE.add_constraints(simulation) diff --git a/examples/tutorials/biopixel_trap.py b/examples/tutorials/biopixel_trap.py index 785a7a0..1df83fc 100644 --- a/examples/tutorials/biopixel_trap.py +++ b/examples/tutorials/biopixel_trap.py @@ -14,6 +14,7 @@ from __future__ import annotations +from dataclasses import replace from pathlib import Path from cellmodeller2 import ( @@ -91,7 +92,7 @@ def _load_prindle_layout() -> tuple[MaskRectangle, ...]: DRAG_COEFFICIENT = 100.0 -def _grid() -> SignalGridSpec: +def _grid(simulation: Simulation | None = None) -> SignalGridSpec: shape = GridShape() shape.x, shape.y, shape.z = 42, 60, 14 grid = SignalGridSpec() @@ -107,7 +108,13 @@ def _grid() -> SignalGridSpec: grid.diffusion = [40.0] grid.advection = [Vec3()] grid.integration = SignalIntegrationKind.CRANK_NICOLSON - DEVICE.apply_to_grid(grid, inlet_values=[NUTRIENT_INLET], outlet_values=[0.0]) + 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 @@ -136,7 +143,12 @@ def _regulate(step: ControllerStep) -> StepPlan: mobility = colony_mobility( GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT ) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) @@ -160,7 +172,8 @@ def _divided(step: ControllerStep, event: DivisionEvent) -> None: def build(context: ModelContext) -> NativeController: simulation = context.simulation(reserved_capacity=20_000) - simulation.configure_signal_grid(GRID, _primed_levels(GRID)) + grid = _grid(simulation) + simulation.configure_signal_grid(grid, _primed_levels(grid)) simulation.set_coupled_rate_plan(_rate_plan()) DEVICE.add_constraints(simulation) diff --git a/examples/tutorials/danino_clock.py b/examples/tutorials/danino_clock.py index fb02fe6..fe4a876 100644 --- a/examples/tutorials/danino_clock.py +++ b/examples/tutorials/danino_clock.py @@ -10,6 +10,7 @@ from __future__ import annotations from collections.abc import Sequence +from dataclasses import replace import numpy as np from cellmodeller2 import ( @@ -100,7 +101,7 @@ DRAG_COEFFICIENT = 100.0 -def _grid() -> SignalGridSpec: +def _grid(simulation: Simulation | None = None) -> SignalGridSpec: shape = GridShape() shape.x, shape.y, shape.z = 64, 72, 6 grid = SignalGridSpec() @@ -116,10 +117,12 @@ def _grid() -> SignalGridSpec: grid.diffusion = [AHL_DIFFUSION, 20.0] grid.advection = [Vec3(), Vec3()] grid.integration = SignalIntegrationKind.CRANK_NICOLSON - DEVICE.apply_to_grid( + 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 @@ -198,7 +201,12 @@ def _regulate(step: ControllerStep) -> StepPlan: mobility = colony_mobility( GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT ) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) @@ -227,7 +235,8 @@ 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, _primed_levels(GRID)) + grid = _grid(simulation) + simulation.configure_signal_grid(grid, _primed_levels(grid)) simulation.set_coupled_rate_plan(_rate_plan()) DEVICE.add_constraints(simulation) diff --git a/examples/tutorials/pillar_channel.py b/examples/tutorials/pillar_channel.py index a9ee084..fc9fe9e 100644 --- a/examples/tutorials/pillar_channel.py +++ b/examples/tutorials/pillar_channel.py @@ -90,7 +90,7 @@ def _in_pillar_core(px: float, py: float, margin: float) -> bool: ) -def _grid() -> SignalGridSpec: +def _grid(simulation: Simulation | None = None) -> SignalGridSpec: shape = GridShape() shape.x, shape.y, shape.z = 22, 60, 4 grid = SignalGridSpec() @@ -122,10 +122,14 @@ def _grid() -> SignalGridSpec: boundary.kind = GridBoundaryKind.FIXED boundary.values = [NUTRIENT_INLET if name == "y_lower" else 0.0] setattr(grid, name, boundary) - field, _ = solve_flow_field( - grid, mean_inlet_speed=FLOW_SPEED, mobility=gap_mobility(grid) - ) - grid.velocity_field = field + 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 @@ -173,7 +177,12 @@ def _regulate(step: ControllerStep) -> StepPlan: mobility = colony_mobility( GRID, step.cells, base=GAP_MOBILITY, drag_coefficient=DRAG_COEFFICIENT ) - field, _ = solve_flow_field(GRID, mean_inlet_speed=FLOW_SPEED, mobility=mobility) + 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) @@ -213,7 +222,8 @@ def _divided(step: ControllerStep, event: DivisionEvent) -> None: def build(context: ModelContext) -> NativeController: simulation = context.simulation(reserved_capacity=10_000) - simulation.configure_signal_grid(GRID, _primed_levels(GRID)) + grid = _grid(simulation) + simulation.configure_signal_grid(grid, _primed_levels(grid)) simulation.set_coupled_rate_plan(_rate_plan()) _add_walls(simulation) diff --git a/python/tests/test_microfluidics.py b/python/tests/test_microfluidics.py index fe56b94..bb2c1b3 100644 --- a/python/tests/test_microfluidics.py +++ b/python/tests/test_microfluidics.py @@ -13,6 +13,7 @@ SignalGridSpec, SimulationController, Vec3, + backend_available, ) from cellmodeller2.microfluidics import BiopixelTrapDevice, TrapChannelDevice from cellmodeller2.runner import build_model @@ -152,6 +153,23 @@ def test_trap_example_builds_steps_and_transports_nutrient() -> None: 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)