Every row names its reproducer. Every caveat is part of the number.
The long form, including the post-mortems that produced these numbers,
is in archive/campaign1/BENCHMARKS.md.
Last full re-check: 2026-10-01 on an Apple-silicon laptop in the pinned
cellsim environment.
15 diverse cocrystals; structures and ligand SMILES fetched from RCSB at run time, so the benchmark is not circular on bundled inputs.
| Metric | Result | Gate |
|---|---|---|
| PoseBusters physical validity | 100 % (15/15) | ≥ 95 % |
| Pose recovery, best of top-3 ≤ 2 Å | 87 % (13/15) | ≥ 75 % |
| Pose recovery, top-1 ≤ 2 Å | 73 % (11/15) | ≥ 75 %, just under |
Reading: Vina's sampling finds the right pose; its scoring ranks it first less reliably. A UFF-strain re-rank does not recover the misses. Earlier numbers (the "2/3 mini-bench") were computed on molecules scrambled by a Meeko atom-order bug; biotin/1STP reported 2.02 Å then, and reads 0.56–0.69 Å now.
Reproduce: python scripts/run_blind_dock_bench.py benchmarks/pdbbind/blind_set.yaml
(results committed in benchmarks/pdbbind/blind_set_results.csv).
benchmarks/dock/streptavidin_calibration.yaml listed desthiobiotin at
K_d = 5×10⁻⁵ M (ΔG −5.9 kcal/mol). That was wrong by about six
orders of magnitude and inconsistent with the source cited on the entry
itself: Hirsch 2002 (Anal Biochem 308:343) is the standard reference for
~10⁻¹¹ M (ΔG −15.0), and Green 1990 (Methods Enzymol 184:51) puts
desthiobiotin 10²–10⁴ fold weaker than biotin, which brackets it. The
value is what makes desthiobiotin useful as a reversible affinity tag:
tight enough to capture, loose enough to elute with free biotin.
Corrected, and the consequences followed through rather than hidden.
| Before (wrong K_d) | After (corrected) | |
|---|---|---|
| Spearman ρ | +0.80 | +0.40 |
| MAE | 4.98 kcal/mol | 6.66 kcal/mol |
The claim this set used to support — that docking ranks biotin-site binders usably — does not hold. The class has been moved out of the "reliable for ranking" table in the tutorial into "pose filter only", alongside kinases.
What replaced it is a stronger test. tests/uq/test_calibration_smoke.py
asserted Spearman ≥ 0.8, which only passed because of the bad number.
It now pins the saturation itself: across four compounds spanning
12.0 kcal/mol of experimental affinity, the predictions span
0.33 kcal/mol, all between −7.1 and −7.4. That is Vina's tight-binder
saturation measured directly, it is reproducible, and it is the property
the reliability table's do_not_trust_absolute verdict rests on.
Downstream artifacts regenerated from the corrected data:
benchmarks/dock/uq_coverage.json and the ultra_tight_binder row of
benchmarks/dock/reliability_table.yaml.
Measured absolute error of docking ΔG against experiment, so a prediction can carry an accuracy estimate, not only seed scatter. Monte-Carlo over seeds understated biotin's real error about 40-fold.
| Target class | n | MAE (kcal/mol) | Verdict |
|---|---|---|---|
| Trypsin-like S1 pocket (3PTB) | 6 | 0.91 | absolute ΔG usable |
| Kinase ATP site (1M17, EGFR) | 6 | 2.16 | rank-order only, and even that is weak (Spearman −0.49) |
| Ultra-tight binders (1STP) | 4 | 4.98 | absolute ΔG meaningless; Vina saturates |
n < 19 per class, so these are typical errors seen so far, not guaranteed conformal bounds. Pooled split-conformal across classes gives ± 11.6 kcal/mol, which is honest and useless; calibrate per class.
Reproduce: python scripts/run_uq_coverage.py; table in
benchmarks/dock/reliability_table.yaml.
| Path | Literature set (aspirin, midazolam, diazepam) |
|---|---|
| xTB C–H BDE ranking | 1/3 |
| + heme-accessibility re-rank (default) | 2/3 |
Systematic blind spot: N-dealkylation sites (diazepam and many tertiary-amine substrates). The tool is advisory and is labelled so.
Reproduce: cellsim som-validate benchmarks/quantum/cyp3a4_som_validation.yaml.
FreeSolv-12 hydration, 11 windows × 50 ps per leg, CPU, tag
milestone-a-pilot-3 (report archived in
archive/campaign1/freesolv_pilot3/):
- 10 of 12 compounds converged (acetic acid and acetamide failed MBAR overlap on the hand-rolled sampler; replica exchange converges them).
- MAE 1.42 kcal/mol, Pearson 0.91, but only 1 of 10 residuals inside the MBAR error bar, and residuals grow with solute size: methane −0.2, ethane −1.6, propane −1.6, benzene −1.9, pyridine −3.1, toluene −3.7.
- Cause identified 2026-10-01: the solvated system is a packmol box run
in NVT; nothing in
cellsim/fepapplies a barostat or equilibrates the density. This is a systematic error, not sampling noise. - Binding ΔG (double decoupling) has never produced a sampled number on a real binder. PR #10 fixed 13 audited defects, including a ~10 kcal/mol sign error in the standard-state correction.
Status: kept as reference code with its smoke gates; not a product path.
Builds headless with cmake -S OLD -B build-old && cmake --build build-old
(no vcpkg needed for the validators). Re-run 2026-10-01:
| Validator | Seed 0x1 | Seed 0xDEADBEEF |
|---|---|---|
cellsim_blind_benchmarks (BM01–BM08) |
8/8 | 8/8 |
cellsim_p53_pulsing (gates a–h) |
PASS | — |
cellsim_cisplatin_arrest |
PASS | — |
cellsim_cisplatin_apoptosis |
PASS | — |
What those passes do and do not mean: the p53 axis responds, oscillates
(7 peaks / 24 h, period 3.03 h against 5.5 h published), arrests the
cycle (100 % G2 at 24 h, stronger than the literature's "dominant"), and
builds the PUMA/BAX/MOMP commitment signature. Cells do not die:
caspase-3 stays at 0.04 under saturated MOMP, so death counts are 0 under
cisplatin in every run. Units are "sim units" in several subsystems
(OLD/UNITS.md). Treat OLD/ as a reference implementation of the
pathway topology and the citations, not as a validated predictor.
Twelve gates in tests/cell/test_engine_smoke.py, all passing, each
asserting a phenotype rather than "it ran" (the two not tabulated below
cover cell-to-cell variability), plus five in
tests/cell/test_stream_smoke.py for the UI stream, including that a
6 h pulse washed out regrows (×4.0 at 72 h) where the same dose held
continuously kills the colony:
| Gate | Result |
|---|---|
| Untreated colony doubles at the line's quoted rate | A549 grows ×10.3 in 72 h vs ×9.7 expected for a 22 h doubling |
| Asynchronous phase fractions are physiological | G1 0.67 / S 0.18 / G2 0.13 / M 0.03 |
| p53 pulses under sustained damage, Purvis-like period | 4.9 h (published 5.5 h; the C++ prototype gave 3.0 h) |
| No pulsing without damage (specificity) | basal p53 stays below 0.12 |
| Damage executes death, not just commitment | caspase-3 reaches 1.0, cells are removed |
| p53-mutant line resists damage-induced death | HT-29 deaths ≪ A549 deaths at the same damage |
| Dose-response is monotone with a finite IC50 | holds |
| Fixed seed is reproducible | exact |
| State stays finite and bounded | holds |
| Cycle scale tracks doubling time | HCT116 (18 h) > MCF7 (29 h) |
The death-execution gate is the one that matters most: it is the defect that made the C++ prototype unusable (MOMP saturated while caspase-3 sat clamped at 0.04, so no cell ever died). The port fixes it by putting Smac in excess over XIAP, as measured (Rehm 2006).
Curve shape is not yet realistic. Each simulated cell carries log-normal variability (σ 0.3) in its Bcl-2 reserve and drug uptake, the documented source of fractional killing (Spencer 2009). That widens the dose window over which viability falls from 85 % to 15 % from 1.2× (identical cells) to about 2×, measured on a 10-point-per-decade grid. A Hill slope of 1 would give about 32×. So the simulated curves are far steeper than typical measured ones, and only IC50s are compared below.
Protocol: each drug's single fit_target constant is fitted by
bisection on the clean TP53 wild-type lines A549 and MCF7 (geometric
mean of the per-line fits) and applied unchanged to every line. HeLa and
both TP53-mutant lines are held out. Results in
benchmarks/cell/gdsc_validation_results.csv and …_summary.json.
Ten cell lines (five TP53 wild-type, five mutant), fitted on A549 and MCF7 only, so eight lines are held out. The table is the first ten-line run, before the apoptotic buffer and the efflux channel were built (see below); the current engine's figures follow it.
| Drug | Fitted constant | Engine in span | Null in span | log10 RMSE engine (held-out) | log10 RMSE null |
|---|---|---|---|---|---|
| cisplatin | k_damage_per_uM_h = 0.000954 |
9/10 | 7/10 | 0.54 (n=5) | 0.47 |
| doxorubicin | k_damage_per_uM_h = 0.00794 |
8/10 | 8/10 | 0.45 (n=8) | 0.46 |
| paclitaxel | partition = 0.153 |
10/10 | 8/10 | 0.31 (n=8) | 0.35 |
| held-out | 22/25 | 18/25 |
With the current engine (buffer and efflux on, library refitted
October 2026 to cisplatin 0.001181, doxorubicin 0.01749, paclitaxel
partition 0.212): held-out in span 20/25 against the null's 18/25;
held-out log10 RMSE 0.65 / 0.43 / 0.32 against the null's
0.47 / 0.46 / 0.35; paired test 10 : 11, p = 1.0
(benchmarks/cell/gdsc_validation_summary.json).
Re-measured at ten drugs in Phase 3 (below): 32 wins to 41 losses over 91 held-out predictions, p = 0.35. The conclusion holds, and the larger panel also found the one exception — a drug whose target is part of the modelled mechanism, where the engine does call lines correctly.
Paired per-line comparison against the null on held-out lines: the engine is closer on 10, the null is closer on 11, two-sided sign test p = 1.0. The engine does not yet tell cell lines apart better than predicting one number for all of them.
This corrects an earlier claim in this file. The gate first used to read "beat the null's held-out RMSE for 2 of 3 drugs", and on a five-line set it passed. That was wrong twice over: the set was too small, and the margin that passed it was 0.45 against 0.46, far inside the noise on eight points. The gate is now the paired sign test, which cannot pass on a difference that small, and the whole set is ten lines.
What is solid. The engine puts each drug's potency on the right scale from one fitted constant: 22 of 25 held-out predictions landed inside GDSC's replicate span against the null's 18 in this run (20 of 25 with the current engine), and A549 cisplatin predicts 9.0 µM against 9.8 measured. Paclitaxel is the best case, 10/10 in span and the lowest error of the three.
What is missing, now measured. Line-to-line variation. The only line-specific inputs are doubling time and TP53 status, and the reference data shows TP53 status only predicts sensitivity for one of the three drugs (table above). Doubling times themselves disagree between sources by up to 1.6× for these lines. So the engine has almost no information with which to separate lines, and it shows.
The p53-independent death route: helps, but not enough. ATM/c-Abl →
p73 → PUMA is the documented route by which TP53-mutant cells still die
of DNA damage (Gong 1999 Nature 399:806; Agami 1999 Nature 399:809).
Added as Params.p73_gain, default 0.1. Tested leave-one-out across all
five mutant lines (scripts/experiment_p73.py): the gain is chosen on
the other four and scored on the one left out, so the test line never
influences it. Held-out log10 error per mutant:
| Held out | p53-only | With route | Null |
|---|---|---|---|
| HT-29 | 0.388 | 0.332 | 0.410 |
| MDA-MB-231 | 0.524 | 0.325 | 0.405 |
| MDA-MB-468 | 0.980 | 0.337 | 0.321 |
| MIA-PaCa-2 | 1.100 | 0.590 | 0.422 |
| T47D | 0.215 | 0.716 | 0.739 |
The route is a large improvement over the p53-only model on four of five mutants, where that model is badly wrong (off by up to 12×). But it beats the null on only two, and it actively hurts T47D, the one mutant that is genuinely drug-resistant: the route kills it when the data says it survives.
That is the whole problem in one row. Treating p53 loss as "these cells still die" fixes the lines that are sensitive and breaks the line that is resistant, because what separates them is not p53 at all. The earlier two-mutant version of this table looked like a clean win and was simply too small.
The sign-test tie is not a calibration problem. The engine has no channel through which cell lines can differ enough to matter. Measured by sweeping each per-line field across its full observed range and recording how far the predicted IC50 moves:
| Per-line field | swept over | cisplatin | doxorubicin | paclitaxel |
|---|---|---|---|---|
| doubling time | 18–38 h | 1.0× | 1.1× | 1.0× |
| Bcl-2 reserve | 0.5–4× | 1.1× | 1.2× | 1.0× |
| DNA repair rate | t½ 32 h–2 h | 9.1× | 7.4× | 1.0× |
| observed between lines | 10.4× | 51.3× | 8.6× |
Three things follow.
- Doubling time, the engine's main line-specific input, has no leverage at all. Sweeping the entire range of the ten lines moves the IC50 by 1.0×. Since TP53 status is the only other per-line input and it is binary, the engine can produce at most two distinct answers per drug. Tying with a single constant is the expected result, not a surprise.
- DNA repair rate is the one channel that works, and it nearly spans what cisplatin needs (9.1× against 10.4×). This is where per-line information should enter for DNA-damage drugs — ERCC1 and the rest of the NER machinery, not proliferation rate.
- Paclitaxel has no line-specific channel whatsoever. Every field gives 1.0×, because its IC50 is set by the tubulin-occupancy arrest threshold, which is a property of the drug and not of the line. To make it line-specific the engine needs per-line efflux (ABCB1 acting on intracellular drug), which it does not currently model.
At ten lines nothing could be validated. Repeating the test on every
GDSC line DepMap has expression for — 156 to 683 lines rather than
6 to 10 — with genes pre-specified by mechanism before any
correlation was computed, and Bonferroni-corrected within each drug
(scripts/experiment_expression_scale.py):
| Drug | n lines | Survives correction, in the expected direction |
|---|---|---|
| doxorubicin | 683 | BCL2L1 ρ +0.19 (p 5e-06), ABCB1 ρ +0.17 (p 9e-05) |
| paclitaxel | 427 | BCL2L1 ρ +0.31 (p 9e-10), ABCB1 ρ +0.23 (p 1e-05), TUBB3 ρ +0.22 (p 3e-05) |
| cisplatin | 156 | nothing |
So real, mechanistically sensible signal does exist once the panel is large enough: the anti-apoptotic set-point (Bcl-xL), P-glycoprotein efflux, and class III β-tubulin for taxanes, which is the textbook taxane-resistance marker. Cisplatin gives nothing at all — and notably the repair genes do not predict cisplatin sensitivity (ERCC1 ρ −0.03, ERCC2 +0.03, XRCC1 −0.06), consistent with ERCC1's long record of failing to replicate as a platinum biomarker.
Put the two measurements side by side. Translating each significant correlation into the IC50 change implied across its 10th-to-90th percentile expression range:
| Channel | What the data implies | What the engine delivers |
|---|---|---|
| anti-apoptotic set-point (Bcl-xL) | 1.8–2.3× | 1.1× |
| efflux (ABCB1) | 1.5–2.1× | not modelled |
| class III β-tubulin | 1.9× | not modelled |
| DNA repair | no signal (cisplatin) | 7–9× |
| proliferation rate | not tested | 1.0× |
The engine has its leverage exactly where the data has no signal, and has little or no leverage exactly where the data does. That is the whole failure, stated in one table:
- Repair rate is the engine's only strong channel, and repair genes do not predict cisplatin sensitivity.
- The apoptotic set-point is the strongest real signal and the
engine already has the field for it (
CellLine.bcl2_level), but its MOMP switch is so sharp that changing that field 8-fold moves the IC50 by 1.1×. The channel exists and is throttled. - Efflux is the second strongest real signal and the engine does not model it at all.
A ceiling worth knowing before building. These effects are modest: roughly 1.5–2.3× each. Even used perfectly and independently they compose to maybe 3–5×, against the 8–51× by which these lines actually differ. So wiring them in should make the engine beat a constant, but it will not close the gap — most line-to-line variation in drug response is not captured by the canonical markers. That is a property of the biology, not of this implementation.
The first item was acted on. The Bcl-2 reserve was modelled as a rate penalty on Bax activation, which is why it barely mattered: a slower climb still crosses a fixed threshold inside 72 h, so the reserve changed when a cell died, not whether. It is now a stoichiometric buffer — only Bax in excess of the reserve can form pores — which is the mitochondrial-priming picture (Certo et al. 2006 Cancer Cell 9:351).
| Before | After | |
|---|---|---|
IC50 span over bcl2_level 0.5–4×, doxorubicin |
1.2× | 1.7× (data implies 1.8×) |
| doxorubicin held-out log10 RMSE vs null | 0.45 / 0.46 | 0.40 / 0.46 |
| paclitaxel held-out RMSE vs null | 0.31 / 0.35 | 0.30 / 0.35 |
| cisplatin held-out RMSE vs null | 0.54 / 0.47 | 0.65 / 0.47 (worse) |
| paired sign test (engine : null) | 10 : 11 | 11 : 10, p = 1.0 |
Predicted and confirmed. The ceiling above said wiring these in "should make the engine beat a constant, but will not close the gap". That is exactly what happened: the paired count flipped from losing to winning, and stayed far from significance. A 1.7× channel cannot resolve an 8–51× spread.
It also improved the model's honesty elsewhere. Treating the reserve as a buffer restored partial resistance in TP53-mutant lines — survival 0.08 against wild-type 0.00 under sustained damage, where the previous version killed mutants outright. That was flagged above as stronger than the literature supports, and is now pinned by a test.
Cisplatin got worse, which fits: it is the drug with no Bcl-xL signal, so the added mechanism only adds noise there.
ABCB1 was the remaining signal the engine could not express, and the only possible route to line-specific behaviour for paclitaxel, whose IC50 is set by a tubulin-occupancy threshold belonging to the drug. A saturable P-glycoprotein pump now acts on intracellular drug for substrate drugs only, scaled by the line's ABCB1.
CellLine.efflux_level is derived, not fitted: 2^(ABCB1 log2(TPM+1)
− panel median) from DepMap 24Q4, a fixed formula with no free
parameter, so held-out lines stay out of sample. At the pump strength
the measured ABCB1 range implies, leverage is 1.9× for doxorubicin and
2.0× for paclitaxel — matching the data's 1.5–2.1× — and exactly 1.0×
for cisplatin, which is not a substrate.
It did not improve prediction. The paired count moved from 11 : 10 to 10 : 11, a one-line change at p = 1.0 — noise, not degradation, but certainly not the improvement the population statistics suggested. The high-ABCB1 lines are the tell: HeLa carries ~10× the panel's median ABCB1 yet sits mid-pack for doxorubicin sensitivity, so the engine now predicts it too resistant.
That failure is not an implementation flaw, and the point generalises. Taking all eleven candidate markers together and fitting them to IC50 by least squares — far more freedom than any mechanistic model gets — across every GDSC line DepMap covers:
| Drug | n lines | best-case R² | spread it must explain |
|---|---|---|---|
| cisplatin | 156 | 0.09 | 4× |
| doxorubicin | 683 | 0.15 | 10× |
| paclitaxel | 427 | 0.20 | 7× |
80–92 % of why one cell line is more sensitive than another is not captured by the canonical markers at all. A correlation of ρ ≈ 0.2 is real at n = 683 and nearly useless for an individual line: it is a population trend, not a predictor. So a mechanistic model restricted to these inputs cannot pass the paired test, however faithfully it is built — and building it more faithfully, as the efflux channel was, does not help.
This is the honest answer to the question Phase 1 was posed to settle, and it redirects the project rather than ending it. Predicting which cell line is more sensitive is not a goal this engine can reach with available inputs. What it does well is the other thing: reasoning about one system over time — dose schedules and wash-out, resistance emerging under selection, the timing of combinations, the difference between arrest and kill. Those depend on mechanism the engine has and on no per-line marker at all.
Still to do, in order. Revised by the measurement above:
- Re-aim Phase 2 at within-line dynamics, which is where the engine's mechanism is load-bearing and where the validation data is not a ten-line IC50 table: schedule dependence, wash-out recovery, resistance under selection, combination timing.
- Keep the per-line channels — apoptotic set-point and efflux are correct biology and are needed for drug-level reasoning — but stop treating per-line IC50 as the headline metric.
- If per-line prediction is revisited, it needs inputs beyond canonical markers (pathway state, baseline apoptotic priming by BH3 profiling, or a much larger panel), not more mechanism.
- A shallower death response (curve shape above).
Re-run cellsim validate-gdsc after each; the gate is the paired sign
test, not the RMSE table.
The library's three fitted potency constants predated the apoptotic buffer and the efflux channel, although its docstring said re-running the fit reproduced them within 1 %. Refitting under the current defaults moved them 1.2–2.2× (cisplatin k_damage 0.000951 → 0.001181, doxorubicin 0.00794 → 0.01749, paclitaxel partition 0.153 → 0.212). That silently moved every Phase-2 dose relative to its drug's potency, and re-running the experiments showed which headlines were properties of the mechanism and which were properties of one dose or one random seed.
Two changes make that class of error structural rather than lucky.
Experiments now state doses as multiples of each drug's own engine IC50,
computed at run time (cellsim.cell.engine.ic50), so a refit cannot
change what a design means. Tests pin the potency their doses were
designed at; cisplatin and paclitaxel enter the engine only as potency ×
concentration, so for them this is exactly a rescaling of the dose.
| Claim as first published | Status | What re-measurement shows |
|---|---|---|
| Pulsing selects a pharmacokinetic escape (low uptake), not an apoptotic one | Retracted: backwards | It rested on ONE surviving lineage for cisplatin and two for doxorubicin. Measured exactly on a trait grid, continuous exposure selects on accumulation and pulsing shifts selection toward the apoptotic reserve, for all three drugs. |
| The sequence effect is a readout artefact (0.93 → 0.98 → 0.99) | Qualified | A single-seed reading; one 128-cell run carries ±0.02–0.03 on that ratio. Over five seeds there is no ordering effect with both drugs at their IC50, and ~12 % at 72 h with cisplatin at 1.25× its IC50, shrinking to ~3 % by 120–144 h. |
| Exposure-time ranking 242× / 17× / 7.6× (3 h ÷ 72 h survival) | Replaced | The ratio depends on the dose chosen, and once no drug acts within 3 h it reduces to 1/S(72 h). The iso-effect curve below replaces it; the ranking holds. |
| Antagonism 1.32–1.50× | Confirmed | 1.29–1.35 (±0.08 over seeds) with both drugs at their IC50, rising to 1.6–1.9 with more cisplatin. |
Everything here is one cell line compared against itself under different schedules, so the per-line ceiling above does not apply.
The measurement is the iso-effect curve: for each exposure time T,
the concentration C50(T) that halves the 72 h cell count when the drug
is washed out at T (A549; scripts/experiment_schedule.py). If killing
were governed by AUC, C50 × T would be constant and log C50 against
log T would have slope −1.
| Drug (72 h IC50) | 1 h | 3 h | 6 h | 12 h | 24 h | 72 h | slope, 1–6 h | shortest exposure reaching half kill |
|---|---|---|---|---|---|---|---|---|
| cisplatin (10.7 µM) | 326 | 120 | 59.1 | 30.6 | 16.3 | 10.4 | −0.95 | 1 h |
| doxorubicin (0.063 µM) | 2.06 | 0.68 | 0.34 | 0.17 | 0.089 | 0.060 | −1.01 | 1 h |
| paclitaxel (0.030 µM) | none | none | none | none | 0.037 | 0.030 | — | 24 h |
(C50 in µM; "none" = no concentration reaches half kill. These are single-seed figures and carry about ±20 % — repeating the run with a different seed moves each by 1.16–1.35×, for the reason set out under "Which published claims that ±20 % actually touches" below. Read them to one or two significant figures. The slope is not affected, because the seed shifts the whole curve together rather than tilting it, which is the point of that section.)
Cisplatin and doxorubicin are AUC-governed over short exposures: C50 × T stays at 326–367 µM·h for cisplatin from 1 to 12 h and at 2.0–2.1 µM·h for doxorubicin from 1 to 24 h, then rises as repair outpaces slow accumulation. Ozawa et al. measured exactly this C × T law for cisplatin (WiDr cells; Cancer Res 1989 49:3823) and time dependence instead for phase-specific agents.
Paclitaxel is time-governed: no exposure shorter than 24 h halves the colony at ANY concentration, because only cells that attempt mitosis while the drug is present are affected. At a fixed AUC both extremes therefore fail, for different reasons — too brief and few cells reach mitosis, too dilute and tubulin occupancy never crosses the arrest threshold — and the best split is interior (C50 × T is 0.88 µM·h at 24 h against 2.1 at 72 h). That is the documented clinical finding that taxane efficacy follows time above a threshold concentration rather than AUC (Gianni 1995 J Clin Oncol 13:180; Huizing 1993 J Clin Oncol 11:2127).
Against published exposure-duration studies, matched by assay:
| Study | Published | Engine | Verdict |
|---|---|---|---|
| Ozawa 1989, cisplatin, WiDr | C × T governs (slope −1) | slope −0.95 over 1–6 h | agrees |
| Liebmann 1993 (Br J Cancer 68:1104), 8 lines, clonogenic | 24 h kill saturates above ~50 nM | 0.46 / 0.43 / 0.40 / 0.49 surviving at 0.05 / 0.2 / 1 / 5 µM | agrees |
| Liebmann 1993 | 24 → 72 h exposure raises cytotoxicity 5–200× | C50 2.7 → 0.027 µM (100×) | agrees, after the October 2026 fix below |
| Georgiadis 1997 (Clin Cancer Res 3:449), 14 NSCLC lines, MTT, read at 120 h: 3 / 24 / 120 h exposure | IC50 > 32 / 9.4 / 0.027 µM | 147 / 2.7 / 0.027 µM (A549) | 3 h agrees; 120 h agrees but is inherited from the GDSC fit; 24 h now within 3.4×, having been 240× — see below |
The miss, and what fixed it (October 2026). Both studies say a 24 h paclitaxel exposure is far weaker than a 72–120 h one; the engine said it was nearly as good, by 240×. Two things were wrong, and finding the second depended on fixing the first.
Drug left the cell as fast as it entered. Uptake was symmetric —
(partition·C_out − C_in)/τ with one τ — so wash-out emptied a cell in
minutes. Paclitaxel does not behave that way: it is held on tubulin, and
after removal there is "little efflux", with an elimination half-life of
13–53 h. Drug.tau_efflux_h now gives exit its own time constant
(24 h for the taxanes, 0 = symmetric for everything else, so no other
drug moved). Without it the 24 h curve had an unreachable floor —
survival stalled at 0.42 from 0.1 µM to 30 µM, a 300-fold dose range in
which more drug did nothing, because a cell that mitosed after wash-out
met no drug at all.
The mitotic death rate was never determined by the fit. With the floor
gone the curve still crossed half-kill far too early. k_mitotic_death_per_h
was fitted to GDSC's 72 h screens — and a 72 h IC50 moves by under 10 %
while that constant varies a hundred-fold, because given three days
even a slow death rate kills every arrested cell. The fit could not see
it. Short exposures can:
| engine before | engine after | measured | |
|---|---|---|---|
| 3 h exposure | 73 µM | 147 µM | > 32 µM ✓ |
| 24 h exposure | 0.039 µM | 2.7 µM | 9.4 µM |
| continuous | 0.028 µM | 0.027 µM | 0.027 µM ✓ |
0.3 → 0.03. The 24 h error falls from 240× to 3.4×, the 3 h and continuous points are unchanged, and paclitaxel's median fold error against the GDSC panel stays at 2.4×.
Read honestly: the 24 h point is now calibrated on Georgiadis rather than predicted, so it is no longer available as a test. The 3 h point is an independent check and it passes. The remaining 3.4× is real and unexplained.
Left alone on purpose: docetaxel and vinorelbine still carry 0.3. They share paclitaxel's mechanism, so the same argument probably applies — but no short-exposure data for them has been found, and changing a constant with nothing to constrain it is how the first wrong value got there. The engine kills almost every cell that reaches mitosis inside the 24 h window. The likely causes are that its cells cycle almost in lockstep (real cycle times vary by 20–30 % and include slow-cycling cells that a 24 h window misses) and that cells slipping out of mitotic arrest re-enter the cycle instead of arresting as tetraploids. Drug retention after washout cannot be the cause, because it would make short exposures stronger, not weaker. Open; recorded so the schedule result is not over-read.
Gated by tests/cell/test_schedule_smoke.py (5 gates, 38 s): the C × T
slope for cisplatin, the 6 h paclitaxel pulse that cannot reach half
kill, the sub-threshold thin split, the interior optimum, and
continuous-beats-washout.
Each representative cell is a lineage carrying its own anti-apoptotic reserve (Bcl-2 multiplier) and drug accumulation (uptake multiplier), and a surviving lineage keeps both as it doubles, so selection emerges without being modelled.
Measured exactly, not sampled. cellsim.cell.selection runs a
13 × 13 grid of trait values (six cells per point, random cycle phases)
through each schedule, records survival as a function of trait, and
integrates that map against the population's log-normal trait
distribution. That gives the survivors' expected means however rare
survival is; a logistic fit gives the boundary's steepness along each
axis, and |b_uptake / b_reserve| says which trait decides survival.
Two weeks, A549, matched total exposure: continuous at the drug's own
72 h IC50 versus 3× IC50 for 24 h then 48 h drug-free
(scripts/experiment_resistance.py):
| Drug | Arm | Expected survival | Survivors' reserve | Survivors' uptake | |b_u / b_r| |
|---|---|---|---|---|---|
| cisplatin | continuous | 0.235 | 1.01 | 0.64 | 5.2 |
| pulsed 1 d / 2 d | 0.0016 | 1.88 | 0.50 | 1.35 | |
| doxorubicin | continuous | 0.070 | 1.21 | 0.56 | 1.8 |
| pulsed 1 d / 2 d | 0.0014 | 2.45 | 0.76 | 0.91 | |
| paclitaxel | continuous | 0.476 | 1.00 | 0.75 | uptake only |
| pulsed 1 d / 2 d | 0.45 | 1.14 | 0.98 | 0.39 |
(1.00 = no selection.)
Continuous exposure selects on drug accumulation. A sustained dose sets a damage steady state, and whether it stays under threshold depends on how much drug a cell takes up; for cisplatin the reserve barely moves (1.01). Pulsing shifts selection toward the apoptotic reserve. A brief high dose pushes nearly every cell over threshold for a while, and what decides survival is whether the reserve can absorb the transient — it can buffer a pulse but not a siege.
Pulsing kills far more of a concentration-driven drug: 143× fewer survivors for cisplatin and 49× for doxorubicin at matched exposure, but 1.06× for paclitaxel, which wastes a compressed exposure because it needs cells to transit mitosis while it is present — the same mechanism as the schedule result, measured independently.
A prediction an experimental lab can test. Pulsed and continuous selection are the two standard ways resistant lines are derived (McDermott et al. 2014 Front Oncol 4:40). The engine predicts that the continuously derived line resists mainly by accumulating less drug (e.g. raised efflux), and the pulse-derived one mainly by a larger apoptotic reserve (lower mitochondrial priming on BH3 profiling).
What this cannot show. Heterogeneity is drawn once, so this is selection from a pre-existing tail, not evolution of new resistance; lineages grow without a carrying capacity, so they never compete, and adaptive-therapy questions are out of reach until the spatial dish.
Gated by tests/cell/test_resistance_selection_smoke.py (5 gates,
18 s), on a 7 × 7 trait grid.
The engine carries one intracellular concentration per drug, so a
combination is a real two-drug run. Every arm gives each drug for 36 h
at the same concentration in the same line; only the timing differs.
Five seeds × 256 cells, mean ± sd over seeds
(scripts/experiment_combination.py):
| Both at 1× IC50 | Cisplatin 1.25× IC50 | |
|---|---|---|
| together ÷ independent expectation | 1.29 ± 0.08 | 1.64 ± 0.07 |
| cisplatin → paclitaxel ÷ independent | 1.34 ± 0.08 | 1.68 ± 0.08 |
| paclitaxel → cisplatin ÷ independent | 1.35 ± 0.06 | 1.92 ± 0.17 |
Antagonism is robust. Every arrangement kills less than independent action predicts, and more so with more cisplatin, because the platinum arrests the cycle through p53 → p21 and removes the cells paclitaxel needs in mitosis — the textbook interaction of a cytostatic with a phase-specific agent. Nothing was fitted to produce it.
The ordering effect is small, dose-dependent, and mostly timing. Cisplatin-first survival divided by paclitaxel-first survival:
| Readout | 72 h | 96 h | 120 h | 144 h |
|---|---|---|---|---|
| both at 1× IC50 | 0.99 ± 0.02 | 1.00 ± 0.02 | 1.00 ± 0.02 | 1.01 ± 0.02 |
| cisplatin 1.25× IC50 | 0.88 ± 0.06 | 0.93 ± 0.05 | 0.96 ± 0.03 | 0.97 ± 0.03 |
With both drugs at their IC50 there is no ordering effect at all. With more cisplatin the platinum-first order looks ~12 % better at 72 h, and most of that disappears once the slow drug given second has had time to act. So the engine does not reproduce a biological platinum/taxane sequence dependence; what it shows is how a fixed-endpoint assay can confound "the order matters" with "the second drug had less time".
Gated by tests/cell/test_combination_timing_smoke.py (4 gates, 41 s):
two-drug state, antagonism, no ordering effect at matched IC50 (3-seed
mean), and the fade of the above-IC50 ordering effect with readout.
The well-mixed engine says how many cells survive; the dish says which ones and where. Each lattice site holds at most one cell and each cell carries a full engine state row, so the validated biology runs per cell unchanged. Space adds contact inhibition (a dividing cell must push its neighbours to free space within a set reach, or its G1 growth signal drops to zero), diffusing oxygen and drug consumed by the cells they pass, hypoxic quiescence, necrosis, and clones that inherit their founder's traits. Two geometries: a monolayer under stirred medium and a 3-D spheroid whose fields are solved radially, which is exact for a symmetric aggregate and is how the spheroid literature models it.
Every environmental constant is an input with a source: oxygen diffusivity, medium pressure and consumption from DLD-1 spheroids (Grimes et al. 2014 J R Soc Interface 11:20131124: 2e-9 m²/s, 100 mmHg, 22.1 mmHg/s in packed tissue); oxygen-dependent proliferation from PhysiCell's documented defaults (Ghaffarizadeh et al. 2018); apoptotic clearance 8.6 h (Macklin et al. 2012); small-molecule interstitial diffusivity 1e-6 cm²/s (Nugent & Jain 1984). How much drug a cell sequesters, which sets how slowly drug penetrates tissue, is known only to an order of magnitude and is flagged as such per drug.
Two engine changes it needed, both inert by default. A per-cell G1
growth signal: lowering cyclin D synthesis alone does not arrest this
cycle, whose Rb/E2F/cyclin E feedback restarts itself, so the signal
scales progression through G1 only — cells past the restriction point
finish their cycle and their daughters wait (gs = 0 leaves exactly the
56 % of cells that were past G1 dividing once and nothing dividing
twice). And optional cycle-time variability (Params.cycle_cv). With
neither in use, every earlier result is bit-identical.
Numerics against closed forms (tests/cell/test_dish_smoke.py). The
oxygen solver reproduces the zero-order sphere profile to within 1 % of
the medium value and Grimes's 233 µm diffusion limit; with first-order
uptake the slab solver reproduces cellsim.cell.tissue's analytic
profile to 0.5 %; the transient drug solver relaxes to the steady one.
A sparse monolayer: the Cell Tracking Challenge HeLa movie. Both
46 h sequences re-run as closed fields of the same size, seeded with the
same counts (five seeds each; scripts/validate_dish.py, targets from
the curated tracks by scripts/ctc_reference.py):
| sequence 01 | sequence 02 | |
|---|---|---|
| cells at 46 h, measured | 43 → 137 | 125 → 363 |
| engine (cycle CV 0.25) | 140 (trajectory error 12 %) | 377 (15 %) |
| complete cycles, measured | 20.2 h, CV 0.22 | 18.5 h, CV 0.33 |
| engine | 26.4 h, CV 0.23 | 27.3 h, CV 0.22 |
| still undivided at 30 h, measured | 0.51 | 0.71 |
| engine | 0.41 | 0.45 |
The colony grows right; the single-cell picture does not. The movie's cells are a mix of fast cyclers (complete cycles ~19 h) and a large fraction that does not divide within the window — the proliferation / quiescence decision at mitotic exit (Spencer et al. 2013 Cell 155:369).
Now modelled (CellLine.quiescent_fraction, 0 for every line but
HeLa). The fix needed two changes together, which is what made the miss
informative:
Re-measured through scripts/validate_dish.py with the library's new
HeLa (five seeds per row), against the committed figures from before:
| undivided 20 h | 30 h | complete cycle | trajectory error | |
|---|---|---|---|---|
| measured, seq 01 | 0.77 | 0.51 | 20.2 h | — |
| 31 h cycle, no quiescence | 0.88 | 0.41 | 26.4 h | 0.12 |
| 19 h cycle + 50 % quiescent | 0.64 | 0.52 | 17.0 h | 0.06 |
| measured, seq 02 | 0.80 | 0.71 | 18.5 h | — |
| 31 h cycle, no quiescence | 0.94 | 0.45 | 27.3 h | 0.15 |
| 19 h cycle + 50 % quiescent | 0.69 | 0.57 | 17.8 h | 0.05 |
The trajectory error halves or better, the cycle times come down from 26–27 h to within a couple of hours of the measured 18–20 h, and sequence 01's plateau is reproduced almost exactly. It now slightly OVERSHOOTS at 20 h — the engine divides a little too fast early — where before it was far too slow.
This is a calibration, not a validation, and the distinction matters. The 19 h and the 0.50 were chosen by scanning against these very sequences, so "the engine matches the movie" is circular: it was tuned to. What survives the caveat is narrower:
- two parameters, eight targets — per sequence it reproduces the count trajectory, the complete-cycle mean, and the undivided fractions at 20 h and 30 h, with no per-sequence tuning;
- the STRUCTURAL claim, which is the real finding and does not depend on
the fitted values: no single cycle time can reproduce both a 46 h
count trajectory and a plateauing Kaplan-Meier curve, because they
demand different things of it. Needing two parameters is what exposed
that
doubling_time_hwas conflating them.
What would settle it is a different HeLa time-lapse, or another line with tracked single cells, held out entirely.
A held-out test of half of it (October 2026). The Cell Tracking Challenge also publishes DIC-C2DH-HeLa — different cultures, different microscopy, never used for any calibration here. At 84 frames of 10 min it spans only 14 h, too short to repeat the test above: a daughter born during the movie needs ~19 h to divide, so no complete cycle and no 30–40 h plateau can be measured. But it can test something the two parameterisations disagree on — what fraction of the cells present at the start divide within 14 h:
| founders dividing within 14 h | against 14 of 15 measured | |
|---|---|---|
| old HeLa: 31 h cycle | 0.56 ± 0.17 | rejected, p ≈ 0.002 |
| current HeLa: 19 h cycle + 50 % quiescent | 0.75 ± 0.12 | consistent, p ≈ 0.08 |
| measured, DIC-C2DH-HeLa (2 sequences) | 0.93 | — |
So independent data rejects the old 31 h cycle and supports the 19 h one. That is the first part of the HeLa change to be validated rather than calibrated.
It covers only half. Quiescence is decided at mitotic exit, so cells present at the start are never quiescent in the engine, and this observable tests the cycle time alone. The 50 % quiescent fraction acts on daughters, whose fate plays out beyond a 14 h window, and remains calibrated on one dataset and untested.
One thing worth recording without over-reading: the measurement sits on
the high side of even the new model. This field is sparse — about 10
founders, against 43–125 in the calibration movie — and sparse cultures
proliferate more freely. If that is the reason, the quiescent fraction
depends on culture density rather than being a fixed property of HeLa,
and CellLine.quiescent_fraction is the wrong home for it. Fifteen
cells cannot settle that; a long, sparse held-out movie could.
scripts/validate_hela_heldout.py downloads the movie, scores it and
re-runs both models, so the comparison is repeated whenever the engine
changes rather than resting on a number recorded once.
Quiescence alone was not enough: it fixed the plateau at 30 h and left the engine far too slow at 20 h. The cycling cells also had to cycle faster — and that exposed a conflation in the library worth stating plainly.
doubling_time_h is a single-cell cycle time, and the library fills
it with population doubling times. calibrate_cycle_scale uses the
field as the time one cell takes from early G1 to mitosis; ATCC and
Cellosaurus quote the time a culture takes to double. They agree only
when every cell cycles. HeLa's cells complete a cycle in ~19 h while its
counts double every 27–31 h, because about half of each generation stops
— so the 31 h in the library makes population growth right and every
single-cell cycle 1.6× too slow, which is precisely the miss above.
For the other lines the gap is small, because ATCC measures exponential-phase culture where quiescence is low. The entry is left unchanged for now rather than quietly corrected: splitting it means re-validating the GDSC predictions that rest on it, and the size of that re-validation is measurable rather than hypothetical.
What quiescence does to each drug class (A549, 72 h IC50, half of each generation quiescent):
| Drug | Mechanism | IC50 shift |
|---|---|---|
| gemcitabine | S-phase | 1.50× |
| sn-38 | TOP1, S-phase | 1.26× |
| cisplatin | DNA adduct | 1.06× |
| doxorubicin | TopII | 1.03× |
| paclitaxel | tubulin | 0.95× |
So changing HeLa's entry would move its antimetabolite predictions by about a third and leave the rest alone — a bounded, checkable impact rather than an unknown one.
The S-phase result is the clinical intuition that antimetabolites spare non-cycling cells, falling out of the mechanism rather than put in. The tubulin result was predicted wrongly before it was measured: a quiescent cell never reaches mitosis, so a taxane should have lost potency too. Looking at the arms shows why it does not. Under gemcitabine the treated colony shrinks less than the control (×0.57 against ×0.46) because quiescent cells escape a drug that kills and are enriched by it. Under paclitaxel both shrink together (×0.44 against ×0.46): a cell that was not going to divide loses nothing when division is blocked. Protection comes from escaping killing, not from escaping cytostasis. With every cell on one clock (cycle CV 0) no cell divides within 30 h of birth at all, so measured variability is clearly needed.
This also caught a provenance error. The library gave HeLa a 20 h doubling time, read off this movie's COMPLETE cycles, which a 46 h window biases short. Cellosaurus gives 1.3 days (PubMed 29156801; DSMZ ~48 h) and the movie's own counts double every 27–31 h, so it is now 31 h. HeLa's held-out GDSC predictions moved 3–7 % and still pass.
A spheroid: DLD-1, calibrated, not yet independently validated.
The dish has one free spatial constant, how far a dividing cell can
push its neighbours (DishParams.push_sites). Grimes's DLD-1 spheroids
grow ~15 µm/day in radius after their core turns anoxic (233 µm to
~400 µm between days 6 and 17), and from 3000 cells over 12 days:
| mechanical reach | necrosis starts | radial growth after that | Grimes |
|---|---|---|---|
| 1 site | day 6 (R 206 µm) | 15.3 µm/day | ~15 µm/day |
| 2 sites | day 4 (R 211 µm) | 28.2 µm/day | |
| 3 sites | day 4 (R 245 µm) | 36.3 µm/day | |
| unlimited | day 3 (R 277 µm) | 58.1 µm/day |
A reach of one site matches, and leaves ~80 % of viable cells quiescent with proliferation confined to the outer layer, consistent with Ki-67 rims; solid stress is known to limit spheroid proliferation this way (Helmlinger et al. 1997 Nat Biotechnol 15:778). Without a mechanical limit the spheroid grows four times too fast.
The necrotic core: a numerical bug, found and fixed (October 2026). Feeding the model's (radius, necrotic radius) pairs through Grimes's own anoxic-core relation should return their 233 ± 22 µm diffusion limit if the model's core tracks anoxia. It returned 173 µm — and, more tellingly, the value drifted down every day (196 → 189 → 181 → 173 → 169 → 163 → 157 µm), although a diffusion limit is a fixed physical property.
This page used to give two likely causes. Both were wrong, and so was a third tried along the way:
| Hypothesis | Test | Verdict |
|---|---|---|
| Formation overshoot: the core outgrows the anoxic region after it stops consuming | compared the dead radius with the currently anoxic radius | rejected — the dead core sits 2–13 µm inside the anoxic region |
| Surface roughening makes the equivalent radius understate the travel distance | compared the volume radius with the mean surface radius | rejected — the surface gets smoother as it grows, and the surface radius gives the same drift |
| Quiescent cells consume less oxygen (Freyer & Sutherland 1985) | swept their consumption | rejected as the cause — it lifts the level but leaves the drift at −6 µm/day |
The clue that cracked it: the engine's own constants give an analytic limit √(6Dc/q) of exactly 233 µm, with Km 1 % of the medium so the kinetics are effectively zero-order. The parameters were right; the solve was not delivering them. And oxygen at the spheroid's surface was falling from 88 to 32 mmHg as it grew — with the unstirred boundary layer set to zero.
The cause was where the medium began. The field applied the medium's oxygen at the distance of the single outermost occupied site. A thin tail of stray cells drifted out, so by day 12 the medium sat 159 µm beyond the spheroid, and all of its oxygen had to diffuse through that much nearly empty fluid first. Real medium is stirred and carries no such gradient. The artificial layer thickened as the spheroid grew — 32, then 93, then 159 µm — which is exactly what produced the drift.
The medium now begins at the equivalent-sphere radius, the same R that Grimes's relation is written in. That choice is made for consistency between the simulation and its own analysis, not tuned toward 233 — a percentile sweep showed the answer moves with the choice, which is why a principled definition was preferred over the best fit:
| Medium boundary at | Implied r_l | Drift | Surface O₂ over the run |
|---|---|---|---|
| outermost straggler (old) | 173 µm | −6.5 µm/day | 88 → 32 mmHg |
| equivalent radius (now) | 240 µm | +1.9 µm/day | 94 → 93 mmHg |
| analytic, from the engine's constants | 233 µm | 0 | — |
| measured, Grimes 2014 | 233 ± 22 µm | — | — |
The strongest evidence it is right is that the oxygen physics stopped
depending on the mechanics. Before the fix the implied limit varied
with push_sites (172, 188, 189, 219 µm for reach 1, 2, 3, ∞, a spread
of 47 µm); after it, 240, 242, 245 and 248 µm, a spread of 8 — physics
now depends only on physics.
What it cost: push_sites = 1 had been calibrated to Grimes's growth
rate under the broken boundary, and the two errors had partly
cancelled, which is how the bug survived. With the surface properly
oxygenated, growth after onset is now 17.1 µm/day against Grimes's
~15 (it was 15.3). Reach 1 is still the best of those tried (reach 2
gives 32), so the choice stands, about 14 % fast. Onset moved from day 6
at 206 µm to day 8 at 238 µm — which is now consistent with a 233 µm
diffusion limit, where the old 206 was not.
What the dish does not do yet. One cell per site, so no compression or multilayering beyond a fixed reach; necrotic debris never lyses; no glucose, lactate or pH; drug sequestration is order-of-magnitude; and spheroid drug response has not been checked against data (multicellular resistance to a 1 h paclitaxel exposure, Nicholson et al. 1997 Eur J Cancer 33:1291, cannot be used: the engine keeps no paclitaxel in a cell after wash-out, and kills nothing with exposures under 24 h).
Gated by tests/cell/test_dish_smoke.py (7 gates, 20 s).
The panel is ten drugs across eleven lines — 91 held-out predictions
against 25 in Phase 1 — chosen so that every mechanism the engine can
express is covered and each has GDSC screens on our lines. Each carries
exactly one fitted constant, fitted on the clean TP53 wild-type lines and
applied unchanged everywhere else (scripts/validate_gdsc.py).
Refitted in 1.5.1 under the p53 damage gate; before that, seven of these constants existed only in the fit and not in the library users ran (see What watching a single cell revealed).
| Drug | Mechanism | Fitted | In GDSC span | Sensitive/resistant call | Spearman |
|---|---|---|---|---|---|
| cisplatin | DNA adduct | k_damage 0.00127 |
8/11 | 5/11 | −0.26 |
| doxorubicin | TopII | k_damage 0.0186 |
9/11 | 11/11 | 0.58 |
| etoposide | TopII | k_damage 0.0102 |
6/10 | 10/10 | 0.54 |
| sn-38 | TOP1, S phase | k_damage 2.19 |
8/11 | 11/11 | −0.09 |
| gemcitabine | antimetabolite | k_damage 1.31 |
10/11 | 10/11 | 0.75 |
| 5-fluorouracil | antimetabolite | k_damage 0.00325 |
9/11 | 6/11 | 0.77 |
| paclitaxel | tubulin | partition 0.214 |
11/11 | 11/11 | 0.06 |
| docetaxel | tubulin | partition 0.274 |
7/11 | 11/11 | 0.00 |
| vinorelbine | tubulin | partition 0.301 |
7/11 | 10/11 | 0.50 |
| nutlin-3a | MDM2 | partition 0.016 |
11/11 | 11/11 (null 5/11) | −0.70 |
The Phase-1 conclusion survives the larger panel, and is now much better measured. On held-out lines the engine is closer than a constant-IC50 null on 31 and further on 42 (sign test p = 0.24), and lands inside GDSC's replicate span 69 times against the null's 68. Ten drugs and 91 predictions say the same thing three drugs and 25 said: the engine reproduces each drug's potency scale and does not rank cell lines better than a single number. That is now a result rather than a small-sample suspicion.
Those figures are medians over three seeds, and switching to them is a small result in itself. Each held-out prediction used to be a single run carrying the ±20 % the engine varies by; now it is a median. The fitted constants did not move at all (≤ 0.6 %), because the bisection's own tolerance was already coarser than the noise it was averaging. The paired comparison did move, and against the engine: 34 : 39 became 32 : 41, p 0.64 → 0.35. Removing a model's own noise is supposed to be able to do that — a fair comparison is the point, not a flattering one — and the negative result is now stated on cleaner numbers at three times the compute.
But the panel also found the exception, and it is the informative one. Nutlin-3a is the only drug whose line-to-line differences are mechanism rather than an unmeasured marker: it blocks MDM2, so p53 accumulates without any DNA damage, and whether that kills depends on whether functional p53 is there to accumulate. The engine gets all eleven lines right, nine of them held out, where the null gets five:
| Line | TP53 | Engine | GDSC | |
|---|---|---|---|---|
| A549, MCF7, HCT116, HT-1080, U-2-OS | wild-type | sensitive | sensitive | ✓ |
| HT-29, MDA-MB-231, MDA-MB-468, MIA-PaCa-2, T47D | mutant | resistant | resistant | ✓ |
| HeLa | wild-type | resistant | resistant | ✓ |
HeLa is the sharp one. Its TP53 is wild-type by sequence, so a model
reading the genotype would call it sensitive and be wrong. HPV18 E6 hands
p53 to the E6AP ligase instead of MDM2, and in cervical cancer cells that
switch is complete (Hengstermann et al. 2001 PNAS 98:1218), so blocking
MDM2 cannot rescue a protein that is no longer being degraded by MDM2.
That route is in the engine as CellLine.p53_mdm2_independent_deg_per_h,
and switching it off makes the engine predict HeLa sensitive — the
specificity check is in the test suite.
Nothing about this was fitted to the p53 split: nutlin's one fitted constant is its partition coefficient, fitted on two wild-type lines.
So the honest statement is narrower and more useful than either extreme. The engine cannot rank lines for cytotoxic drugs, because the differences there live in markers it cannot see. It can call which lines respond to a drug whose target is part of the modelled mechanism, including a line whose genotype points the wrong way. That is the case for building more mechanism rather than more markers.
Gated by tests/cell/test_mdm2_p53_smoke.py (4 gates, 8 s).
Palbociclib blocks cyclin D–CDK4/6, so Rb stays unphosphorylated and the cell stops in G1. Where there is no working Rb there is nothing to stop, and it "does not inhibit the growth of Rb-deficient cells" (Fry et al. 2004 Mol Cancer Ther 3:1427; breast lines lacking Rb were fully resistant in Finn et al. 2009 Breast Cancer Res 11:R77). It is the second drug, after nutlin, whose line-to-line differences are mechanism rather than a marker.
Two library lines lack Rb, for different reasons, and the engine now
knows (CellLine.rb_functional): MDA-MB-468 has lost RB1 (its
transcript is 4–8× below every other line in DepMap), and HeLa's
HPV18 E7 binds Rb and cuts its half-life from over 6 h to 2–3 h (Oh et
al. 2009 Virology 396:118) — the same virus whose E6 makes HeLa resist
nutlin. Neither line was used in the fit (A549 and MCF7, GDSC1).
| Line | Rb | Engine IC50 | GDSC | |
|---|---|---|---|---|
| HeLa | E7 | none | > 30 µM (both screens) | ✓ held out |
| MDA-MB-468 | lost | none | > 32 µM (both screens) | ✓ held out |
| A549 | yes | 1.5 µM | 0.68–7.5 µM | fit |
| MCF7 | yes | 2.6 µM | 1.5 µM | fit |
| HCT116, HT-29, MDA-MB-231, MIA-PaCa-2 | yes | 1.1–3.7 µM | 1.0–3.0 µM | ✓ held out |
| HT-1080 | yes | 4.4 µM | 1.2 µM | outside the span |
| U-2-OS | yes | 10.6 µM | 2.0 µM | ✗ called resistant |
| T47D | yes | > 100 µM | 1.2 µM | ✗ called resistant |
The mechanism calls are right; the values are fragile. Both Rb-less lines are predicted to ignore the drug at any dose and both did, with nothing fitted to them. The two misses have identifiable causes:
- T47D — at 72 h a cytostatic drug cannot push a count below the cells already committed to divide, and in this engine that floor sits at about half the untreated count (0.47–0.52 in every Rb-positive line at saturating dose: the cells past the restriction point divide once, and some later slip through as cyclin E–CDK2 builds without CDK4/6). An IC50 is the 50 % line, so whether one exists at all turns on a few per cent: MCF7 just crosses it, T47D just does not. GDSC shows the same fragility from the other side — GDSC1 found IC50s of 1–3 µM in most lines, while GDSC2, in the same 72-hour format, called most of them out of range. For a drug that stops cells rather than killing them, a 72-hour IC50 is close to the wrong readout, and the engine's floor may also be somewhat high; how deep and how lasting a CDK4/6 arrest is needs time-lapse data before it is tuned.
- U-2-OS — palbociclib is a P-gp substrate (Parrish et al. 2015 JPET 355:264; de Gooijer et al. 2015 Int J Cancer 137:2007), flagged before any fit was run, and U-2-OS has the panel's third-highest ABCB1 level (5.2×, after DLD-1 and HeLa). The same line was the one nutlin's P-gp flag broke. Either the line's efflux level is too high for drugs other than the ones it was set from, or the pump is less protective in these cells than the engine's saturable model makes it; one line cannot say which.
Panel totals with palbociclib: held-out span 75/100 (null 75), sensitive/resistant calls 86/100 (null 80), paired comparison 31 : 49 (p = 0.057). The ten drugs before it are unchanged in substance (span 69/91, calls 79/91 — one better than before, from the Rb correction to HeLa and MDA-MB-468 — paired 31 : 42); palbociclib's tight GDSC1 values (all 1–3 µM where measured) suit a constant null, and it adds seven of those losses.
The RB1 knockout in cellsim.cell.perturb now means what the line flag
means. Rb's ODE variable is the active, unphosphorylated fraction and
its "production" is dephosphorylation, so zeroing it only stopped Rb
being reactivated: under a CDK4/6 inhibitor nothing phosphorylated it,
the old Rb stayed active, and the "knockout" arrested like the parent
line. It now sets rb_functional=False, as a TP53 knockout sets
p53_functional=False.
Gated by tests/cell/test_cdk46_smoke.py (4 gates).
A prediction without an interval is not usable, and an interval that has
never been scored is not an interval
(scripts/calibrate_cell_uncertainty.py). The width comes from the
spread of held-out log10 errors; the test is leave-one-drug-out, so
each drug is scored with a width built only from the other nine. Scoring
a width on the errors that built it would read as well calibrated however
badly it generalised.
| Nominal | Engine interval | Engine coverage | Null interval | Null coverage |
|---|---|---|---|---|
| 68 % | ±2.6× | 64 % | ±2.7× | 70 % |
| 90 % | ±4.8× | 85 % | ±5.0× | 90 % |
| 95 % | ±6.5× | 95 % | ±6.8× | 93 % |
Every level lands within 15 points of nominal, so the stated interval means what it says. Two things to read with it. The intervals are wide — a factor of five at 90 % — and that is the honest width given that most line-to-line variation is not in the markers the engine can see; GDSC's own replicate screens of the same line and drug disagree by a median 5.3×, which is the floor nothing can honestly go below. And the null's intervals are just as narrow at the same coverage, which is the same finding as the paired test, now in the uncertainty picture: the engine's point predictions are not better than a constant per drug, so neither are its error bars.
Everything before this could select from the spread a population already
had: heterogeneity was drawn once, so every lineage that survived had
been resistant from the start. DishParams.mutation_rate lets a daughter
differ heritably from her mother, and the control that makes the claim
testable is a homogeneous start — identical cells, so selection has
nothing to act on (scripts/experiment_evolution.py, A549 + cisplatin,
seven days, fold against the untreated arm of the same condition):
| Start | Mutation | Schedule | Survivors | Accumulation | Reserve | IC50 |
|---|---|---|---|---|---|---|
| identical | off | continuous | 0 | — | — | — |
| identical | off | pulsed 3× | 0 | — | — | — |
| identical | on | continuous | 230 | 0.61 | 1.03 | 1.85× |
| identical | on | pulsed 3× | 0 (10 mutations first) | — | — | — |
| spread 0.3 | off | continuous | 2535 | 0.68 | 1.05 | 1.29× |
| spread 0.3 | on | pulsed 3× | 27 | 0.56 | 1.75 | 2.18× |
Four things follow, and the first is the control:
- Identical cells with no heritable change are eradicated. No resistance is possible when there is neither standing variation nor mutation, which is what makes the rest of the table mean something.
- Resistance arises where nothing was there to select. From identical cells, continuous exposure leaves survivors that accumulate 40 % less drug and need 1.85× the dose.
- A hard intermittent schedule can win the race. The same total exposure given as brief 3× pulses eradicates the identical population after only ten mutations have occurred — resistance never gets started.
- Where standing variation exists, pulsing leaves the most resistant remnant, and selects on the apoptotic reserve (1.75) rather than on accumulation. That is what the corrected Phase-2 selection measurement found by a completely different method.
Fold-resistance reaches 1.3–2.2×, against the two- to eight-fold that derived lines actually reach (McDermott et al. 2014 Front Oncol 4:40), so the assumed mutation rate is if anything conservative. It is a modelling choice, not a measurement: resistance here arises from many routes at once whose combined rate per division is not a published constant, and the script reports what a given rate produces rather than claiming one.
Gated by tests/cell/test_evolution_smoke.py (4 gates, 18 s), including
the control that identical cells without mutation cannot drift.
Everything above uses our cell lines and our fitted constants. A lab has its own line, its own reader and its own drug, and until their numbers can enter the simulator none of this is a tool for them. This is the first half of that door: get a plate in, normalised, fitted, with an error bar that has been checked. (Fitting the simulator's constants to those data is the second half, and it needs this uncertainty, so it comes after.)
cellsim plate --data reading.csv --map map.csv reads a reader's grid
export or a tidy table, subtracts blanks, sets the plate's own controls
to 1, fits a four-parameter logistic in log concentration, and reports
GR metrics when the plate carries a time-zero read.
Why GR as well as IC50. A 72 h viability assay confounds potency with division rate: a line that doubles twice in the assay has more to lose than one that doubles once, so the same drug looks more potent in the faster line. GR metrics divide by the control's own growth over the same window (Hafner et al. 2016 Nat Methods 13:521) and are comparable across lines and assay lengths. They need one extra read, at time zero.
The only way to know what a confidence interval is worth is to simulate
curves whose IC50 is known and count how often the interval contains
it (scripts/validate_plate_fit.py, 200 plates per row):
| Method | Coverage at a nominal 95 % |
|---|---|
| resample the raw residuals | 80 % |
| scale residuals by √(n/(n−p)) | 88–92 % |
| bias-corrected and accelerated (BCa) | 88–90 %, and no bias in the point estimate |
The first gap has a definite cause: residuals from a fitted curve are smaller than the true errors, because the fit has already absorbed p of the n degrees of freedom. The rest is what a small-sample bootstrap on a nonlinear parameter does, and ~90 % is where ten concentrations land. That number is printed next to every interval rather than left implied.
This is the part worth giving a biologist before they run the assay.
Concentrations per curve (noise sd 0.05):
| Points | Coverage | Interval width |
|---|---|---|
| 5 | 68 % | 1.85× |
| 7 | 86 % | 1.83× |
| 10 | 88 % | 1.75× |
| 12 | 84 % | 1.68× |
A five-point curve fits four parameters to five numbers. The interval it
returns looks no wider than a ten-point one and contains the truth two
times in three. fit_4pl therefore attaches an explicit note below seven
points, and summary() quotes 68 % rather than 90 % there.
Noise (ten points): coverage holds at 88–90 % up to a residual sd of 0.10, and falls to 84 % at 0.20 — where the interval has widened to 18× and is no longer telling anyone anything.
Range, which turned out to matter more than either (50 plates per row, since these fits are slow by construction):
| Tested range | Coverage | Interval width | Flagged extrapolated |
|---|---|---|---|
| 0.01–100 µM (brackets the IC50) | 86 % | 1.67× | 0/50 |
| 0.1–30 µM (brackets it) | 80 % | 1.76× | 0/50 |
| 0.001–0.5 µM (clips it from below) | 64 % | 121× | 20/50 |
| 10–1000 µM (clips it from above) | 44 % | 168× | 20/50 |
A range that sits entirely to one side gives a flat plate with no midpoint to find, and is refused rather than fitted — the earlier code fitted it anyway, taking minutes per call to extract a confident IC50 from noise.
The in-between case is the one to understand. A range that merely
clips the response still fits, and the interval does not hide the
problem: it blows up to a hundredfold and more, which is the method
saying it does not know. What it does not do is contain the truth —
even that enormous interval misses it between a third and a half of the
time, because the data constrain the curve's shape only on one side and
the extrapolation can go anywhere. So a wide interval here is not a
conservative answer, it is an absent one. Those fits are flagged
extrapolated, and the flag, not the width, is what to act on.
Fitting a constant so the simulated IC50 matches a measured one proves
nothing: the constant was chosen to make that number right. The claim
worth testing is the next one — that the calibrated model then says
something true about a condition it never saw — and
scripts/validate_calibration.py measures it with ground truth we
control. Perturb a drug's constant to a value the calibrator is never
told; simulate the noisy 72 h plate those cells would give; hand over
only the fitted IC50 and its interval; then ask about a 12 h pulse
washed out long before the readout, and check whether the predicted band
contains the truth.
It failed the first time, and the failure was the useful part. Calibrated bands contained the held-out truth in 3 of 7 trials. A band advertised as covering the answer was wrong more often than right, and nothing upstream was broken: the plate fit's coverage had been measured, the calibration's propagation had been measured, and each piece passed its own test. What had never been tested was the composition.
The cause is that the band counted the measurement's uncertainty and ignored the model's own.
Repeating the identical simulation with a different seed moves the engine's IC50 by 1.1–1.2×. That is not counting noise, and the measurement that settles it is switching the cell-to-cell heterogeneity off:
Params.het_sigma |
24 cells/dose | 96 cells/dose |
|---|---|---|
| 0.30 (default) | 1.187× | 1.108× |
| 0 (identical cells) | 1.009× | 1.008× |
With identical cells the engine is deterministic to under a per cent. With the default heterogeneity the spread is ten to twenty per cent and quadrupling the cells barely moves it, because a population's IC50 is set by which resistant lineages happened to be drawn, not by how many cells were counted — and a tail converges far more slowly than a mean.
That is a property of the model rather than a defect, and it is also the right biology: a real colony's survival is dominated by its resistant minority (Spencer et al. 2009 Nature 459:428, the same observation the heterogeneity was built from). But it has three consequences that matter to anyone using these numbers:
- A single-seed IC50 is good to about ±20 %, so differences smaller than that between two runs are not differences.
- Any calibration from one measurement inherits that, which is why the recovered constant scatters between 0.7× and 1.0× of the truth even with no measurement noise at all.
- More cells is the wrong remedy; more seeds is the right one.
predict_band therefore spans the calibrated interval and seeds,
and the hold-out is scored against the truth's own median over the same
seeds, since no band can honestly be narrower than the quantity it is
predicting.
A spread that large in every IC50 is alarming until it is checked against the claims that rest on it, so here is that check rather than a reassurance. Five seeds, A549:
| Claim | Across seeds | Verdict |
|---|---|---|
| Nutlin calls 11/11 lines correctly (Phase 3 headline) | 0 of 11 lines change call; sensitive lines give viability 0.00, resistant 0.85–1.11 | safe |
| Paclitaxel cannot halve the colony below 24 h at any dose | half kill reached in 0/5 seeds at 3, 6 and 12 h; 5/5 at 24 h | safe |
| Cisplatin follows C × T over 1–6 h (slope −1) | slope median −0.974, range −1.00 to −0.88; the −0.95 quoted above sits mid-range | safe |
| The C50 values themselves (e.g. 318 µM at 1 h) | spread 1.16–1.35× between seeds | quoted too precisely |
The pattern is consistent and it is worth stating as a working rule for this engine. A comparison made within one seed is far more robust than either number in it. The seed fixes which lineages were drawn, and that draw shifts a whole dose-response curve up or down together rather than tilting it — so a slope, a ratio between two arms, or a threshold with a wide margin survives, while the absolute concentration does not.
That is why the slope holds to ±0.06 while the C50s it is computed from move by a third, and it retroactively justifies how the Phase-2 and Phase-3 experiments were built: every one of them compares arms inside a run rather than quoting a number from one.
What it does NOT justify is the single C50 figures quoted in the schedule table above. They carried an unstated ±20 % until this audit; the table now says so, and they should be read to one or two significant figures.
Twenty trials, with the band spanning the calibrated interval and seeds:
| before the fix | after | |
|---|---|---|
| band contains the held-out truth | 3/7 (43 %) | 16/20 (80 %) |
| median band width | — | 0.435 surviving fraction |
| the truth's OWN spread across seeds | — | 0.173 |
| constant recovered | 0.79–0.83× | 0.84× median (0.69–0.98) |
Read the second and third rows together. The band is 0.44 wide and the
quantity it predicts moves by 0.17 on its own between runs, so roughly
40 % of the width is irreducible: no honest band here can be much
narrower. What that buys a user is an answer like "between a quarter
and two thirds of the cells survive this pulse" — enough to rule out
nearly all die and nearly all live, not enough to separate 35 % from
45 %. Two schedules whose bands overlap have not been shown to differ,
which is the whole reason predict_band returns a range.
One thing tried and rejected. If the engine's viability carries ±20 % run-to-run spread, averaging the calibration criterion over seeds ought to help. It does not: three seeds returned byte-identical constants in every case tried, at three times the cost, and moved the offset below by under half a per cent. The reason is that the offset is definitional, not noise, so averaging noise cannot reach it. The option is kept in the signature and the default is one seed. What does move with the seed is the answer itself — the same target calibrated under different seeds spans about 1.14× — and that is what the prediction band spans.
A residual bias, measured and left in place. The fitted constant
lands about 10 % below the truth, one-sidedly, and with the measurement
noise removed entirely it is still 0.88–0.92×. Two definitions of "the
IC50" are in play: calibration targets the concentration where median
viability is one half, while engine.ic50 interpolates a crossing on a
log grid and on a curve this steep lands ~1.1× low. Targeting the
estimator instead would cost a full concentration series per iteration —
the thing that makes calibration fast enough to use — to chase a bias
that sits well inside the ±20 % the engine varies by anyway. It is
recorded rather than removed.
Writing real plate files for the tests, rather than passing arrays directly, found both:
- A grid export pads its unused wells, and
emptywas being classified as a blank. Forty-two zero-reading filler wells were averaged into the background, which silently rescaled every viability on the plate.emptyis now its own role and is ignored. - Readers write
B07where plate maps writeB7. Unnormalised, every well fails to match and the error surfaces far from its cause; well names are now normalised, and an unmapped well raises instead of being dropped.
Gated by tests/cell/test_plate_smoke.py (11 gates, 9 s), which builds
plate files on disk in both layouts.
Added October 2026. All four default to OFF, and the first gate for each is that an untouched run is unchanged — four merged phases of results depend on that.
The engine already carries the pathway as state variables, so a genetic perturbation is a constraint on a mechanism that is already there rather than a new one: a knockout stops a node's production and lets what is there decay, a knockdown scales it, an overexpression raises it. Thirteen genes map to engine nodes; one the engine does not represent is REFUSED, because a knockout that silently does nothing looks like a result.
What makes these checkable is that each has a right answer to land on:
| Check | Result |
|---|---|
| TP53 knockout vs the library's p53-mutant lines | lands in the same place, via the same p53_functional route |
| PUMA vs BAX vs CASP3 knockout | identical to within 0.03 — they are one linear pathway, so losing any one must block it the same way |
| ABCB1 overexpression vs doxorubicin | protects (it is a P-gp substrate) |
| ABCB1 overexpression vs cisplatin | no change — ABCB1 does not transport it, and GDSC's own correlations agree |
The last pair is the real test. A perturbation framework that moves everything is not modelling anything.
One prediction worth recording, because it is the kind of thing the engine exists to produce: blocking apoptosis converts killing into arrest. With BAX gone, viability at the top cisplatin dose floors at 0.14 instead of reaching 0 — the cells stop dying and do not resume growing.
A second diffusing nutrient solved by the same quasi-steady reaction-diffusion code as oxygen, with its own diffusivity (Casciari 1988) and uptake (Casciari 1992). Below ~0.3 mM interior cells stop cycling; below ~0.06 mM they die. Whichever of oxygen and glucose is scarcer governs, because Freyer & Sutherland's central finding (1986 Cancer Res 46:3504, 46:3513) is that a spheroid's necrotic core can be set by glucose rather than oxygen — so neither can be the only route.
DLD-1 spheroid, 6 days from 3000 cells:
| Medium glucose | Live | Necrotic | Radius |
|---|---|---|---|
| 5.5 mM (physiological) | 10 840 | 105 | 207 µm |
| 0.4 mM | 10 830 | 105 | 207 µm |
| 0.2 mM | 6 903 | 0 | 177 µm |
| 0.1 mM | 3 488 | 615 | 149 µm |
Growth stalls before killing starts, which is the ordering starvation should have.
scripts/validate_glucose.py tests the directional claims of
Mueller-Klieser, Freyer & Sutherland (1986): lowering medium glucose
from 16.5 to 0.8 mM substantially thins the viable rim.
With glucose acting only by starvation, it failed. Across the whole published range the model gave byte-identical answers — 175 µm rim, 162 µm core at every concentration — because glucose penetrates ~260 µm and at a 311 µm radius the centre only fell to 0.47 mM, against a 0.30 mM threshold. Not one cell of 29 685 was starving.
The same papers report a second effect alongside the thinner rim: respiration rises as glucose falls. Casciari, Sotirchos & Sutherland 1992 (J Cell Physiol 151:386) measured EMT6/Ro oxygen consumption nearly doubling between 5.5 and 0.4 mM. Starved of glucose, cells burn more oxygen to make the same ATP — so the spheroid goes hypoxic sooner and the rim thins, without any cell running out of glucose at all. The effect runs through oxygen.
With that coupling (DishParams.respiration_boost):
| Medium glucose | Rim at 20 % O₂ | Rim at 5 % O₂ |
|---|---|---|
| 16.5 mM | 175 µm | 58 µm |
| 5.5 mM | 174 µm | 58 µm |
| 1.65 mM | 135 µm | 46 µm |
| 0.8 mM | 111 µm | 39 µm |
All five directional claims pass, within the published range, monotonically, at both oxygen levels.
Why this is a test and not a fit. The coupling's two anchor points (×1 at 5.5 mM, ×2 at 0.4 mM) come from Casciari 1992; the claims it is checked against come from Mueller-Klieser 1986 — different papers, different measurements. The only free choice is the curve's shape between the anchors, which is labelled as such in the code.
What is still not checked: the magnitudes. Mueller-Klieser's rim thicknesses are image-only tables in the scanned originals and could not be extracted, so the directions are validated and the numbers are not. Digitising those figures would settle it.
A random walk on the lattice, optionally biased up the oxygen gradient. This one has an exact answer to check against — a walk of hop rate k on spacing h has D = k·h²/4 in two dimensions:
| Hop rate | Simulated D | Theory |
|---|---|---|
| 0.5 /h | 43 µm²/h | 50 |
| 2.0 /h | 198 µm²/h | 200 |
| 5.0 /h | 414 µm²/h | 500 |
The first attempt gave 126 against a theoretical 200, because one hop per field interval capped the effective rate at 1/interval. Migration now runs in sub-rounds short enough that wanting two hops is rare.
Effector cells on their own layer, as T or NK cells crawl over a monolayer. Each random-walks biased toward tumour cells and delivers hits on contact.
The first version used the wrong mechanism. Each contact killed
outright — one lethal hit — which is the older picture. Live imaging by
Weigelin et al. 2021 (Nat Commun 12:5217) overturned it: a single
contact rarely kills. Each delivers a sublethal hit — perforin pores,
nuclear envelope rupture, DNA damage — that the target repairs, and
death needs about three hits within ~50 minutes of one another,
usually from several CTLs passing through. The engine now implements
that, with effector_hits_to_kill = 1 keeping the old picture available
for comparison.
What came out was not designed in:
| E:T | Lysis | Hits delivered | Hits per kill |
|---|---|---|---|
| 0.1 | 0 % | 418 | 418 |
| 0.3 | 2 % | 1310 | 164 |
| 1.0 | 40 % | 3886 | 17 |
| 3.0 | 100 % | 1890 | 4.1 |
At low density hundreds of hits are wasted on targets that repair between them. At high density the cost falls to about four hits per kill, close to the theoretical three. That steep nonlinearity is why CTLs swarm, and it emerges from the timing rule rather than from any fitted rate.
scripts/validate_killing.py checks four predictions that follow from
the mechanism and so cannot be tuned into passing — all four pass:
-
Cooperativity: hits per kill falls steeply with density.
-
The repair window matters, in the right direction:
Repair window Lysis at E:T 1 15 min 6 % 30 min 18 % 50 min (measured) 40 % 120 min 96 % -
One-hit killing is the old model: it gives 39 % lysis at E:T 0.1 where three-hit gives 0 %, and 87 % at E:T 0.3 where three-hit gives 2 %.
-
No free kills: kills never exceed hits ÷ 3.
What is still not validated is the absolute lysis at a given E:T against a measured chromium-release or impedance curve. None has been digitised here, so the hit rate and walk speed remain literature values and the magnitudes are unchecked. The mechanism is validated; the numbers are not.
Every check before this looked at one moment: viability, or an IC50, at
72 hours. The "Inside a cell" view (cellsim/cell/trace.py, web/lab/)
shows what each simulated cell does minute by minute, and the first
thing it showed was two time courses that are wrong although the 72-hour
endpoints they lead to are calibrated.
Median time to each step, A549, six cells:
| p53 doubles | PUMA rises | death | |
|---|---|---|---|
| cisplatin 20 µM | never | 7 h | 30.6 h |
| nutlin-3a 10 µM | 0.2 h | 0.5 h | 5.4 h |
| nutlin-3a 2 µM | 0.2 h | 0.5 h | 5.5 h |
Nutlin is roughly five to ten times too fast. p53 doubles in twelve minutes and cells die in five hours. In real cells p53 accumulates over several hours after nutlin, and the apoptosis it drives takes one to two days — in A549 it is mostly arrest. The identical timing at 2 and 10 µM is not itself a fault: nutlin binds MDM2 with a constant near 90 nM, so both doses saturate it fully. The speed is the fault, not the dose response.
Cisplatin barely moves p53. It never reaches twice its resting level, where DNA damage normally raises p53 several-fold. Cells still die on a plausible schedule, but much of the killing runs through the p53-independent route the engine models.
Neither is patched here. Both need measured time courses — p53 accumulation after nutlin and after cisplatin, and the timing of apoptosis — and a constant changed without data to constrain it is how the paclitaxel death rate was wrong for a month. Until then the Lab shows a Known limitation note on both drugs, beside the explanation the viewer is already reading, rather than letting a wrong timeline pass as a true one.
Looking for the measured time courses turned up three things, and the first was not a modelling question at all.
1. Seven of the ten drugs were running placeholder constants. The
GDSC fit (scripts/validate_gdsc.py) writes each drug's fitted constant
to benchmarks/cell/gdsc_validation_summary.json, and every score on
this page is computed with those. Everything a user runs — the Lab, the
CLI, cellsim.api — reads cellsim/cell/library.py, and only cisplatin,
doxorubicin and paclitaxel had ever been copied across. The other seven
carried the provisional values they were added with:
| Drug | Field | Library | Fitted | Library ÷ fitted |
|---|---|---|---|---|
| nutlin-3a | partition |
1.0 | 0.013 | 77× too potent |
| sn-38 | k_damage |
0.01 | 2.07 | 1/207 |
| gemcitabine | k_damage |
0.01 | 1.23 | 1/123 |
| etoposide | k_damage |
0.001 | 0.0096 | 1/10 |
| 5-fluorouracil | k_damage |
0.0005 | 0.0031 | 1/6 |
| vinorelbine | partition |
0.2 | 0.30 | 1/1.5 |
| docetaxel | partition |
0.2 | 0.27 | 1/1.4 |
With the library's nutlin, A549's IC50 is 0.063 µM; GDSC measures 4.2
and 6.9 µM. That alone is most of "nutlin is too fast": at the fitted
constant p53 rises 2.9-fold rather than 22-fold and the median death
moves from 6 h to 40 h. It is also what made 2 and 10 µM look identical —
the explanation given above, that both doses saturate MDM2, was wrong;
at the library constant every dose did. The validation page graded the
right numbers while the product ran the wrong ones, and nothing tied the
two together. tests/cell/test_library_constants_smoke.py now does: it
fails when a library constant differs from its fit by more than 2 %, and
when engine parameters have changed since the fit was made.
2. p53 killed without DNA damage, and it should not. At the right constant nutlin still killed 75 % of A549 and HCT116 cells by 72 h. Measured, it kills few of them:
| Measured | Source | |
|---|---|---|
| A549, 72 h nutlin-3 | 7.6 ± 1.1 % annexin-positive | Tovar et al. 2006 PNAS 103:1888 |
| HCT116, 72 h nutlin-3 | 8.8 ± 2.1 % | Tovar 2006 |
| SJSA-1 (MDM2 amplified), 72 h | 88 % | Tovar 2006 |
| HCT116 and U2OS, 24 h nutlin 5 µM | arrested in G1 and G2, alive | Shen et al. 2008 Cancer Res 68:8260 |
| A549, 5-25 µM nutlin-3 | "limited" apoptosis, attributed to the missing damage signal | Deben et al. 2015 Oncotarget 6:22666 |
| HCT116, cisplatin at its IC50 | about half the cells die, spread over 72 h | Paek et al. 2016 Cell 165:631 |
The last row is the one that decides the mechanism. The same HCT116
cells die in half when p53 rises after cisplatin and rarely when it rises
after nutlin, so the difference is in how p53 was raised, not in the
line. p53 raised by damage carries kinase marks — Ser46 among them — that
steer it to its apoptotic genes (Oda et al. 2000 Cell 102:849); p53
raised by blocking MDM2 lacks them. The engine now gates p53's
pro-apoptotic activity on the same ATM signal that engages p73
(Params.p53_apoptosis_damage_free): with no damage signal, p53 keeps
that fraction of its push towards PUMA; with full ATM activity, all of
it. Arrest through p21 is not gated.
The fraction is the one new fitted number, fitted to Tovar's A549 and
HCT116 values with nutlin's own constant refitted to GDSC at every step
(scripts/fit_p53_damage_gate.py):
| gate | A549, 10 µM | HCT116, 10 µM | A549, 5 µM | HCT116, 5 µM |
|---|---|---|---|---|
| 1.0 (ungated) | 69 % | 57 % | 13 % | 14 % |
| 0.6 | 40 % | 33 % | 8.8 % | 10 % |
| 0.4 | 23 % | 16 % | 3.9 % | 4.0 % |
| 0.3 | 12 % | 8.6 % | 0.8 % | 1.2 % |
| 0.2 | 3.1 % | 2.1 % | 0.1 % | 0.1 % |
With the gate at 0.3 and every drug refitted under it, the single-cell view reads (A549 and HCT116; traced cells and whole populations):
| Before | After | Measured | |
|---|---|---|---|
| nutlin 10 µM, dead by 24 h | 100 % | 0 % | HCT116 alive and arrested (Shen 2008) |
| nutlin 10 µM, dead at 72 h, A549 / HCT116 | 100 % / 100 % | 12 % / 8.6 % | 7.6 % / 8.8 % (Tovar 2006, fitted) |
| HCT116 phases at 24 h in nutlin | all dead | G1 7, S 4, G2 13 of 24 | G1 and G2 arrest, S emptied (Shen 2008) |
| nutlin p53 | ×22, doubling in 12 min | ×3.2, doubling in 30 min | — |
| cisplatin 12.5 µM, HCT116, dead at 72 h | — | 42 %, deaths spread over 24-72 h | about half, spread over 72 h (Paek 2016, held out) |
The last row is the check nothing was fitted to: the gate was fitted on nutlin, cisplatin's constant on A549's GDSC IC50, and HCT116's response to cisplatin in single cells lands near Paek's.
The GDSC panel did not get worse for it: held-out predictions inside the GDSC span 69 of 91 (67 before), sensitive/resistant calls 78 of 91 (unchanged), nutlin still 11 of 11, paired comparison against the null 31 : 42 (32 : 41 before, p = 0.24). The Phase-1 conclusion stands as it was.
0.3 assumes Tovar's dose was 10 µM of the active enantiomer. The full text could not be retrieved to confirm it; if it was racemic nutlin-3 (half active), the same fit gives 0.6. 0.3 is the choice that kills less, agrees with "limited apoptosis" in all three nutlin papers, and puts 10 µM — GDSC's top dose — where the measurements are. It is fitted, not validated: Tovar is the data it was fitted to.
3. Cisplatin's p53 was not the fault it was reported as. "DNA damage normally raises p53 several-fold" was written from memory. Measured across twelve p53 wild-type lines, DNA breaks raise p53 between 1.25- and 5-fold (Stewart-Ornstein & Lahav 2017 Sci Signal 10:eaah6671), with HCT116 accumulating steadily and A549 in a broad wave. The engine's 1.5- 1.7-fold is inside that range, at its low end. The Lab's note now says that instead of calling it a limitation.
A fourth, tried and reverted. Building drug look-up for the Lab
(cellsim/cell/compound.py), the P-glycoprotein "rule of fours" called
nutlin-3a a substrate where the library said it was not, and P-gp does
transport it (Michaelis et al. 2009 Cancer Res 69:416). Flagging it and
refitting moved U-2-OS — ABCB1 at 5.2 times the panel's typical level —
from a predicted 5.8 µM to 13.4 µM against GDSC1's 2.5 µM, and nutlin's
held-out error doubled. The same paper explains why: nutlin-3 also
blocks P-gp, at concentrations that leave viability untouched
(vincristine IC50s in P-gp-overexpressing lines fell 92- to 3,434-fold),
so at its active dose it disables its own efflux. The engine's pump has
no self-inhibition, so "substrate" would have meant "protected", which
the data says it is not. The flag stays off, and the rule's docstring
now says what it predicts: transport, not net efflux.
Extracted by scripts/gdsc_reference.py into
benchmarks/cell/gdsc_reference.csv from the public GDSC1 and GDSC2
fitted dose-response tables, release 8.4: 60 screens across the three
drugs and ten cell lines (HCT116 and SW480 are in the library but have
no screens for these drugs). The raw files are not vendored; the
script's docstring has the two-line download.
| Drug | Source | Screens | In tested range | Quantitative target? |
|---|---|---|---|---|
| paclitaxel | GDSC2 | 10 | 10 | Yes |
| doxorubicin | GDSC1 | 20 | 19 | Yes, with the replicate spread below |
| cisplatin | GDSC1 + GDSC2 | 30 | 10 | Partly. Two-thirds of its screens never reached 50 % kill, so those lines give lower bounds only. |
Camptothecin is not a usable stand-in for doxorubicin: in GDSC2 four of five lines have IC50 above its 0.1 µM top dose.
Does TP53 status predict sensitivity? Only for one drug. Geometric mean IC50 over the lines with an in-range reference (October 2026, ten lines: five TP53 wild-type, five mutant):
| Drug | p53 wild-type | p53 mutant | mutant / wild-type |
|---|---|---|---|
| doxorubicin | 0.059 µM (n=5) | 0.125 µM (n=5) | 2.1× more resistant |
| paclitaxel | 0.039 µM (n=5) | 0.037 µM (n=5) | 0.94× — no effect |
| cisplatin | 5.74 µM (n=4) | 2.46 µM (n=2) | 0.43× — mutants more sensitive |
Only doxorubicin behaves as the textbook p53 story predicts. Paclitaxel showing nothing is mechanistically right: it kills through mitotic arrest, not DNA damage, so p53 should not matter. Cisplatin runs the other way, on just two mutant lines, both of which happen to be sensitive. So a model is right to make p53 matter for DNA-damage drugs and not for tubulin drugs, and is wrong to treat p53 loss as blanket resistance.
Replicate spread sets the gate. GDSC1 screened doxorubicin twice per line. The two in-range IC50s of the same line differ by:
| Line | Spread |
|---|---|
| HeLa | 1.5× |
| A549 | 4.9× |
| HT-29 | 6.3× |
| MCF7 | 9.3× |
Median 5.6×. A model cannot be held to "within 3×" when the reference data disagrees with itself by more than that, so the Phase-1 gate is the replicate span (at least 3× either side of the geometric mean where there is one screen) plus the correct most/least-sensitive ordering.
The p53 split is not clean in this data and must not be assumed: cisplatin is less potent in p53-mutant HT-29 (122 µM) than in wild-type A549 (20 µM) as expected, but p53-mutant MDA-MB-231 (43 µM) is more sensitive than wild-type MCF7 (100 µM). The p53 effect is a measured outcome to report, never a gate.
22 standalone tests, each asserting an analytic limit (steady state
reduces to each component, transient converges to steady state, mass
conserved, Monte-Carlo agrees with interval arithmetic). No calibration
against cell-line data yet; that is Phase 1 of PLAN.md.
smoke.yml: 54 always-on steps. The 51 that existed before the October
wiring were replayed locally on the merged tree: 51 pass, 0 fail, 14 min.
The five newly wired steps (29 test files) pass in seconds each.