Censored Weibull–Gaussian biosensor emission and palbociclib heritability analysis (#1016)

[aarmey]

Implements #1016: a right-censored Weibull + Gaussian biosensor emission, the Weibull M-step solver, and the heritability / dose-sweep / early-biomarker analysis. The analysis runs on two public Spencer-lab lineage datasets, and the findings below include where they fall short of what the issue hoped to answer.

⚠️ The data sources named in the issue do not exist

• zenodo.10498701 is an astrophysics thesis ("Identifying Nature's Largest Telescopes with the AbacusSummit N-body Simulations").
• zenodo.12479168 redirects to a record deleted as spam.
• github.com/sc-reporters/multiplex-imaging is a 404, and so is the sc-reporters user/org.
• Web search now returns this issue's own text as if it confirmed these accessions, so search hits for them are circular.

I found no public dataset that combines MCF7/T47D, a palbociclib dose series, lineage tracking, and a per-cell CDK2/DHB trace. The issue's wording most resembles Zikry et al. 2024 (Purvis lab; T47D at 0/10/100 nM palbociclib), but that is fixed-cell 4i imaging, with no time axis and no lineages.

Data used instead

A literature search turned up five public Spencer-lab deposits, all EBI BioStudies (S-BSST21, 167, 231, 314, 439). All are MCF10A with 12-min frames and the same tracedata/genealogy format. Their multi-day drug-escape movies (Yang 2021 Nat Commun, Hoffman 2023 Sci Signal, 90–96 h on drug) are "available upon request" only.

1. Palbociclib: Miller et al. 2018, Cell Reports, Fig. 3A (S-BSST167). DHB-mVenus, 248 frames, drug at frame 135; vehicle vs 1 µM palbociclib.

• One dose, not a sweep; MCF10A, not an ER+ line.
• ~22.6 h after drug, about one generation.
• The frame interval (12 min) is assumed; the deposit doesn't state it.
• The issue's 48 h forecast lead time is impossible in this window; I used 12 h.

2. MEK/ERK inhibitor exposure series: Min et al. 2020, Science, Fig. 1 (S-BSST314). DHB, 257 frames at a stated 12 min.

• The plate map gives MEKi or ERKi on for 1, 3, 6, 9 h, or left on, plus no-drug wells. That makes a six-level exposure series per drug.
• ~32.6 h after drug. The first frame after addition (0-based 94) is taken from a handling dip in the sensor that shows in every well.
• Deeper lineages: in the no-drug wells, 43% of cells born after the drug are granddaughters of the root or deeper, and 9% great-granddaughters or deeper. Palbociclib had essentially none.
• The movies hold 2k–18k lineages per condition (~530k cells). EM time is dominated by per-lineage overhead, so each condition is subsampled to 2,500 lineages (~10k cells).

What's added

Emission and fitting

• lineage/states/CensoredWeibullGaussian.py: N(x | μ, σ) · f_W(t)^δ · S_W(t)^(1−δ).
  • Observations are [x, t, δ]. NaN drops a factor; a negative t is the usual cross-validation mask. It also handles end-of-experiment censoring for simulation.
  • Optional truncation columns [x, t, δ, t_lo, t_hi]: the Weibull factor is divided by S(t_lo) − S(t_hi), for cells in the data only because their lifetime ended inside a window. NaN means no bound; further columns are ignored.
  • When several conditions are fit at once, emissions are shared, so a state means the same phenotype at every dose.
• lineage/states/weibullFit.py: the weighted censored Weibull MLE.
  • Without truncation, the scale is profiled out in closed form, and the shape solves a strictly decreasing score equation by Newton-Raphson on log κ, safeguarded with bisection.
  • With any truncated weight, L-BFGS-B over (log κ, log λ) with an analytic gradient. It starts from the better of the closed-form fit and the previous iterate, so EM stays monotone.

Baum-Welch options (BaumWelch.py, Analyze.py; defaults unchanged)

• shared_T=False: a separate transition matrix per condition.
• independent_T=True: all rows of T equal, i.e. the no-inheritance null.
• estimate_pi=True: the exact EM update for the root-state distribution.
• random_init=False: warm starts; restarts=: number of random restarts.
• An emission can set split_phases = False so that a six-column observation isn't read as G1/G2 phases.

The default (π = stationary distribution of T) is not an EM step and can lower the likelihood. That broke nesting in the likelihood-ratio tests, which is why the analysis uses estimate_pi=True.

Root selection. Every lineage root is in the data because it divided between drug addition and the end of the movie. build_lineages(..., roots=) handles this in one of three ways:

• "truncate" (default): truncates the root's lifetime to that window. With free π this is the exact likelihood given selection.
• "keep": takes root lifetimes at face value (the previous behavior).
• "drop": removes root lifetimes.

Analysis modules

• lineage/heritability.py:
  • persistence half-life ln ½ / ln T_kk; memory eigenvalue λ₂ and its half-life;
  • dose_sweep: per-dose, shared, and identical-row T fits with parallel, cross-seeded restarts, then LRTs for dose dependence and heritability;
  • a lineage bootstrap and exact per-edge mother–daughter posteriors.
• lineage/early_biomarker.py: out-of-fold forecasts of daughter escape from the mother's (x_p, t_p).
  • Folds hold out whole lineages. Forecasters: the mother's sensor, her lifetime, a per-dose logistic regression, and the tHMM.
  • AUCs come with lineage-bootstrap CIs, reported per dose.
• Loaders:
  • lineage/palbociclib_loader.py and lineage/mitogen_loader.py parse the raw .mat movies.
  • They ship per-cell tables (265 KB and 4.6 MB), so the 1.3 GB and 5.9 GB archives are only needed to rebuild them.
  • x = maximum G1 CDK2 activation rate; lost tracks are right-censored.
• lineage/palbociclib_analysis.py → figure 22 (plus sensitivity, which refits under all three root treatments).
• lineage/mitogen_analysis.py → figure 23. It fits each drug's exposure series, with a bootstrap (96 resamples), refits under the other two root treatments, and CDK2-rate correlations between relatives as a model-free check.

Tests: 128 total, all passing. They cover:

• the likelihood against scipy, including truncated densities that integrate to 1 and a gradient checked by finite differences;
• M-step optimality with and without truncation, and removal of selection bias;
• EM monotonicity with free π, nesting of the LRTs, and per-condition and identical-row T;
• a simulated drug time course in which truncated roots keep a slow state slow, while taking root lifetimes at face value collapses it onto the window;
• both loaders against synthetic movies.

ruff and ty are clean.

Findings

Palbociclib (S-BSST167)

Data: 6,697 control and 1,913 palbociclib lineages (20,709 / 5,761 cells). Two states:

• sensitive/arrested: CDK2 rate 0.067 ± 0.052 per hour, essentially never divides in the window;
• tolerant/cycling: 0.31 ± 0.20 per hour; κ = 3.7, λ = 25 h.

BIC keeps falling up to 4 states, which with ~26k cells is largely extra Gaussians fitting the skewed rate distribution, so I didn't interpret K > 2.

Main fit, roots truncated; 95% CIs from 200 lineage resamples:

	control	1 µM palbociclib
T₂₂ (tolerant → tolerant)	0.61 [0.57, 0.64]	0.12 [0.07, 0.24]
τ½ tolerant (generations)	1.41 [1.24, 1.57]	0.33 [0.26, 0.48]
T₁₁ (sensitive → sensitive)	0.99 [0.94, 1.00]	0.39 [0.19, 0.98]
λ₂ (memory eigenvalue)	0.60 [0.52, 0.64]	−0.49 [−0.70, 0.21]

Dose-dependence LRT: p = 5×10⁻⁴¹.

How roots are treated decides whether heritability shows up:

roots	heritability LRT p	λ₂ control	λ₂ palbociclib
truncate	1×10⁻²¹	0.60	−0.49
drop	3×10⁻⁸	0.62	0.22
keep (previous main fit)	0.10	0.17	0.09

• Control memory is robust. The two treatments that account for selection agree (λ₂ ≈ 0.6). The previous "not heritable" result came from taking selected root lifetimes at face value.
• Palbociclib memory is not identified. Its sensitive row rests on the few arrested cells that divide under drug. T₁₁'s CI spans 0.19–0.98, and λ₂ changes sign between root treatments.
• Truncation does not fix calibration. For cells born into drug, the fraction divided by 20 h is 15% observed vs 21% predicted in control, and 2.4% vs 9.5% under palbociclib; that is no better than "keep". So the remaining misfit is in the two-state emission, not root selection.

Early biomarker (palbociclib, n = 187 mother–daughter pairs). Pooled AUCs are confounded, because escape is 75% in control vs 26% under drug, so these are per dose:

• logistic on the mother's CDK2 rate and lifetime: AUC 0.66 [0.55, 0.76] (control 0.57);
• the mother's CDK2 rate alone is inversely associated, AUC 0.31 [0.21, 0.42]; absent in control (0.51);
• tHMM: 0.59 [0.48, 0.69], up from 0.47 with untruncated roots, but not distinguishable from chance.

MEK / ERK inhibitor exposure series (S-BSST314)

Emissions separate the same two states as in the palbociclib data: CDK2 rate 0.05 vs 0.29 per hour, and mean lifetime effectively infinite vs 18 h.

exposure	none	1 h	3 h	6 h	9 h	left on
MEKi T₂₂	0.99	0.90	0.86	0.78	0.59	0.18
MEKi λ₂	0.98	0.85	0.85	0.75	0.58	0.01 [−0.15, 0.15]
ERKi T₂₂	0.98	0.94	0.91	0.91	0.80	0.29
ERKi λ₂	0.97	0.83	0.78	0.79	0.75	0.05 [−0.79, 0.20]

T₂₂ bootstrap CIs are ±0.01–0.05 (figure 23c–d).

• Tolerant-state persistence falls steadily with exposure for both drugs; MEKi's half-life drops from 49 to 0.4 generations. This is the dose-dependence the issue asked for, with exposure duration as the dose.
• The state is strongly heritable (MEKi p = 5×10⁻⁹⁵, ERKi p = 7×10⁻¹⁰⁶), and this holds under all three root treatments (p ≤ 10⁻³⁸).
• Memory is lost only when the drug is left on: λ₂ ≈ 0 in the continuous-exposure wells, while any pulse of up to 9 h keeps λ₂ ≥ 0.58.
• Model-free check (no-drug wells, CDK2 activation rate):
  • sisters r = 0.52 and first cousins r = 0.26, i.e. correlation that decays over generations;
  • mother vs daughter r ≈ 0, which fits a state set at or around division (Min et al.'s mother-history model) better than the mother's own G1 rate carrying over directly;
  • sisters and cousins also share field and birth time, so this is not free of confounding.
• Calibration (divided by 20 h among cells born into the condition):
  • good at long exposures: 14% observed vs 14% predicted for MEKi left on;
  • underpredicts cycling in no-drug wells (85% vs 61%): one shared cycling state (mean 18 h) is a compromise across exposures, while no-drug cells cycle in ~14 h.
  • The refits with roots dropped give single-condition outliers (MEKi 3 h λ₂ = 0.15, left on −0.86). Those refits use fewer restarts (4 vs 8), so I read these as poor optima.

Not done

• Palbociclib dose series: no public data exists (see above). The sweep code takes any number of doses.
• Multi-generation drug escape (≥3 generations on drug): needs the on-request movies (Yang 2021, A375 dabrafenib; Arora 2017, MCF10A aphidicolin).
• K > 2 for S-BSST314, which would likely fix the no-drug calibration.
• Faster E step: vectorizing across lineages would allow fitting all ~530k S-BSST314 cells instead of a subsample.
• Figure checks: figures 22–23 need the committed JSON results, so there is no figure test.

🤖 Generated with Claude Code

https://claude.ai/code/session_01MTbZnG5sNwnZBR9uQerm2D