Model cell lifetime and fate as competing risks

[aarmey]

Addresses #648. Nothing is switched over — this adds the alternative emission and the means to compare it.

What prompted this

#648 asks how death times are distributed and notes they look exponential. Farnaz's follow-up identifies the blocker: death times are not independent of phase durations because of censoring, so fitting them directly is not straightforward.

That blocker is real, and it points at the fix rather than away from it. The dependence exists because death and division are two clocks racing each other and we only ever see the winner. Writing the likelihood that way makes the dependence part of the model instead of a problem to work around.

The gap in the current emission

StateDistributionGamma.logpdf has:

# Remove dead cells
ll[x[:, 0] == 0] = 0.0

A cell that dies contributes only its Bernoulli fate; its duration is zeroed, and the matching mask in estimator drops it from the Gamma fit. A cell dying at 2 h and one dying at 90 h are indistinguishable to the model.

There is a second, less visible problem. A time-censored live cell is scored with sp.gamma.logsf of the division clock, i.e. P(division > t). The correct statement is P(division > t and death > t) — that cell also did not die. Because a death is scored as a bare atom with no time attached while the censored branch ignores death entirely, the two overlap and the emission is not normalized:

horizon	Gamma/Bernoulli mass	competing-risks mass
4 h	1.098	1.000
12 h	1.074	1.000
24 h	1.028	1.000
96 h	1.000	1.000

(compare_emissions.outcome_mass, covered by a test.) The excess is largest for short-horizon cells — which is exactly the ~37% of cells censored in each phase.

What this adds

lineage/states/StateDistributionCR.py gives each phase two latent clocks — division T_D and death T_X — and observes min(T_D, T_X) plus which one fired:

transition seen at t    f_D(t) S_X(t)
death seen at t         f_X(t) S_D(t)
censored at t           S_D(t) S_X(t)

Consequences:

• Death times carry information instead of being discarded.
• The censored term is correct and the emission integrates to 1.
• P(divide) is derived, ∫ f_D(t) S_X(t) dt, rather than fit as a free Bernoulli parameter. The death fraction and the death timing are forced to agree.

The M step is no harder. In the complete-data likelihood the two clocks separate, so each is a weighted right-censored Gamma fit over the same cells — once with "transition observed" as the event indicator, once with "death observed". Both reuse the existing gamma_estimator.

params keeps the Gamma/GaPhs layout in its leading six entries, so figure code that indexes it positionally (common.py:956-960 is the only such place) works unchanged.

Choosing the death clocks from the data

Pooled over the lapatinib and gemcitabine conditions (7,716 cells, 1,724 lineages):

	n	mean	Weibull shape
G1 death	287	18.1 h	0.89
G1 transition	3,818	19.6 h	1.72
G2 death	252	41.6 h	2.72
G2 division	2,341	26.6 h	2.67

So the read in #648 holds for G1: shape ≈ 1, a memoryless constant hazard. Across the full 114-file set it is 1.03, and it falls toward and below 1 as dose rises (0.98 at lapatinib 250 µM, 0.69 at gemcitabine 30 µM).

G2 is a different process. Shape 2.6–4.0 in every well, with deaths arriving later than divisions (41.6 h vs 26.6 h) — a wear-out hazard, consistent with arrest-then-catastrophe rather than a constant risk of dying. Fitting one exponential to "death time" across both phases would wash this out.

G1 therefore gets a one-parameter exponential death clock and G2 a two-parameter one. G1 costs no extra degrees of freedom versus the Bernoulli it replaces; only G2 adds one, taking dof() per state from 6 to 7. New figureS18 documents this.

Evaluation

Setting the comparison up honestly

The two emissions are densities over different observation spaces — the competing-risks form puts a density on when a cell died, which the Gamma/Bernoulli form cannot express. Comparing raw likelihoods would charge the new model for predicting strictly more. (Raw training LL is indeed lower for CR, −17,226 vs −16,488 on lapatinib; that number is meaningless as a comparison.)

So the headline metric coarsens death timing away, leaving a space both models describe:

outcome	competing risks	Gamma/Bernoulli
transition at t	f_D(t) S_X(t)	p f_D(t)
death	1 − P(divide)	1 − p
censored at c	S_D(c) S_X(c)	S_D(c)

The transition branches carry identical total mass, so that part is like for like. The censored branch is where the Gamma/Bernoulli form claims probability it is not entitled to — so this comparison is conservative, biased in the old model's favour.

Held-out likelihood

25% of cells masked per replicate; both models fit on identical masked data and scored on identical held-out sets (lineage/compare_emissions.py). Higher is better; these are sums over ~1,300 held-out cells per replicate.

condition	states	Gamma/Bernoulli	competing risks	difference	replicates won
gemcitabine	2	−3,916.5	−3,847.2	+69.2	3/3
gemcitabine	3	−3,946.6	−3,885.7	+60.9	3/3
gemcitabine	4	−3,854.8	−3,809.3	+45.5	2/2
lapatinib	2	−4,005.1	−3,997.9	+7.2	2/3
lapatinib	3	−3,986.3	−3,982.2	+4.1	2/3
lapatinib	4	−3,927.9	−3,920.5	+7.4	1/2

Competing risks is ahead in every cell of the grid. On gemcitabine the margin is large and every replicate agrees; on lapatinib it is small and the replicates split, so that is a wash rather than a win.

Three replicates per cell except k=4, which has two — I stopped those runs early.

Number of states selected

Held-out likelihood by number of states, per model:

model	gemcitabine	lapatinib
Gamma/Bernoulli	k=2: −3,916 · k=3: −3,947 · k=4: −3,855	k=2: −4,005 · k=3: −3,986 · k=4: −3,928
competing risks	k=2: −3,847 · k=3: −3,886 · k=4: −3,809	k=2: −3,998 · k=3: −3,982 · k=4: −3,921

The two models select the same k — 4 in both conditions, and they agree on the ordering of k=2 against k=3 as well. Swapping the emission does not change how many states the data supports.

(Within a k the two models see identical masks, so the model comparison is exact. Masks differ between k values, so the comparison across k carries more noise than the numbers suggest.)

Fate and death timing

condition	states	fate LL, Gamma	fate LL, CR	difference	CR death-time LL per death
gemcitabine	2	−284.0	−290.2	−6.2	−4.146
gemcitabine	3	−278.5	−292.0	−13.5	−4.037
gemcitabine	4	−276.5	−279.2	−2.8	−4.086
lapatinib	2	−171.8	−171.5	+0.3	−4.280
lapatinib	3	−182.6	−206.5	−23.9	−4.366
lapatinib	4	−162.3	−179.9	−17.6	−4.319

Fate log-loss is slightly worse under competing risks in most cells, which is expected: its division probability is derived from the two clocks rather than fit freely, so it trades a little fate accuracy for internal consistency and death timing.

The death-time density is around −4.0 to −4.4 nats/death, against −4.56 for a uniform over the 96 h window and −4.10 for an oracle per-phase gamma fit in sample on the same deaths. The death clock is well calibrated, not merely non-degenerate. The Gamma/Bernoulli model cannot produce this number at all.

State assignments

Viterbi assignments on the full unmasked data, against the Gamma/Bernoulli model:

	Rand	adjusted Rand
lapatinib, k=2	0.913	0.794
gemcitabine, k=2	0.865	0.729
lapatinib, k=3	0.791	0.582
gemcitabine, k=3	0.629	0.228

At k=2 both emissions recover recognisably the same two states — the division clocks match closely and the death clocks then separate them further. In lapatinib the death-prone state's G1 death scale is 275 h against 51,000 h for the other; the Bernoulli could only express that as a fate probability.

At k=3 they diverge, badly so for gemcitabine. Death timing is new information and it moves the state boundaries. Anyone switching the emission should expect published state assignments to move with it.

Reading of the result

The competing-risks form wins on gemcitabine and is roughly a wash on lapatinib. That tracks how much death information each dataset carries: gemcitabine has 393 deaths in 4,537 cells (8.7%), lapatinib 187 in 5,290 (3.5%). Where deaths are plentiful the new terms pay for themselves; where they are scarce the extra parameter roughly cancels the gain.

Given that the comparison is stacked against it, and that it also fixes the normalization defect and derives rather than fits the division probability, I read this as worth adopting — but the k=3 assignment shift means it should not be swapped in silently.

Also in this branch

An independent bug fix, in its own commit and separable from the rest. StateDistribution.estimator built filtered copies of every observation array to drop cross-validation-masked cells, then passed the unfiltered arrays to bern_estimator. Masking negates a cell's whole observation, so a hidden dividing cell reads as −1: finite, so counted in the denominator, but never equal to 1, so absent from the numerator. With 25% of cells hidden that pulls a true 0.90 division probability to 0.67. Only the one-condition M step was affected — atonce_estimator already filtered before estimating, which is why the cross validation above (four conditions at once) is unaffected by it.

Limitations

• The G2 death shape is weakly identified for states with few deaths. In the lapatinib k=3 fit one state's G2 death clock came out shape 128 / scale 0.45 — a near-deterministic 58 h death, which is not credible. The shape is currently shared across concentrations but free per state; sharing it across states too, the way the division shape already is, would regularize this. Worth doing before anyone reads G2 death parameters quantitatively.
• Left truncation is still treated as right censoring. Root cells, and cells first seen mid-G2, have partial durations; both models treat those as right-censored at the observed partial time. Conservative and unchanged here, but it is an approximation, and it touches the 64 G2 deaths flagged as censored.
• Fitting is 2–4× slower, from the second gamma_estimator call per state per phase.
• Not evaluated on the MCF10A growth-factor data, only the two drug series.

🤖 Generated with Claude Code