Metadata-Version: 2.4
Name: fram
Version: 1.0.0
Summary: First-passage / extreme-statistics kinetic proofreading model of T cell antigen discrimination (Morgan & Lindsay 2023)
Author-email: Thomas Mourier <mourier.t@northeastern.edu>
License-Expression: MIT
Project-URL: Zenodo, https://doi.org/10.5281/zenodo.22999081
Keywords: kinetic proofreading,T cell activation,antigen discrimination,first-passage time,extreme statistics,Weibull,Gillespie
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.8
Requires-Dist: mpmath>=1.2
Requires-Dist: matplotlib>=3.5
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Dynamic: license-file

# fram

Kinetic proofreading with extreme statistics — reproduction and extension of
Morgan & Lindsay, *Modulation of antigen discrimination by duration of immune
contacts in a kinetic proofreading model of T cell activation with extreme
statistics*, PLoS Comput Biol 19(8): e1011216 (2023).

`fram_weibull.py`, next to the package in the source archive, is the original
prototype and remains the reference implementation; the package reproduces its
results and is pinned against them in the test suite.

## Installation

```
pip install fram
```

Python 3.9 or later, with numpy, scipy, mpmath and matplotlib. The test suite,
the prototype and `validate_new_closed_form.py` are in the source archive (the
PyPI sdist, the same file as the Zenodo record); from its root:

```
python -m pip install -e ".[dev]"
python -m pytest                  # 465 test cases, ~2 min
```

## The model

One TCR/pMHC signalling pair, a continuous-time Markov chain:

```
U --k1--> C_0 --kp--> C_1 --kp--> ... --kp--> C_nKP     (C_nKP absorbing = triggered)
C_i --kd--> U          for i = 0 .. nKP-1
```

`kd = k_minus1` for agonist, `sigma * k_minus1` for self antigen. `T_{1,n}` is
the minimum over `n` independent pairs of the first-passage time to `C_nKP`.

```python
from fram import Chain, quantile_min, survival_min

agonist = Chain.agonist(nKP=3)                    # kd = k_minus1
self_ag = Chain.self_antigen(nKP=3, sigma=100.0)  # kd = sigma * k_minus1
quantile_min(0.5, n=1e4, chain=self_ag)           # 71.16 s
```

`Chain.start` is always explicit — it sets the Weibull shape parameter (`nKP`
for pairs that begin engaged in `C_0`, the paper's assumption; `nKP+1` for pairs
that begin unbound in `U`). With every rate equal to 1.0 that difference is
indistinguishable from an off-by-one in `nKP`, which is why it is never a silent
default.

## Four routes to the same distribution

| route | call | cost | use it for |
|---|---|---|---|
| exact | `chain.survival_min(t, n, chain)` | 0.32 ms | everything, any `n`, including `n > 1e6` |
| reference | `chain.F_expm(t, chain)` | ~30 ms | independent cross-check (mpmath, 60 dps) |
| arbiter | `chain.F_rational(t, chain)` | ~ms | exact rational arithmetic; short times only |
| sampling | `ssa.sample_min_activation(n, chain, n_rep, seed=...)` | 478 ms | validating the exact machinery |
| limit law | `asymptotics.weibull_*(…, n, chain)` | 0.01 ms | Eq (1), the `n -> infinity` Weibull |
| paper's ODE | `ode.ode_cdf_min(t, n, chain)` | 0.9–8 ms | the second-order-binding model, where no product formula exists |
| sampling (fast) | `renewal.sample_min_activation_renewal(n, chain, n_rep, seed=…)` | 1–5 s / 1500 reps | simulating where Gillespie cannot go (large `sigma`, large `nKP`) |
| density | `chain.density_min(t, n, chain)`, `log_time_density_min` | as `survival_min` | the PDF, exactly -- `f = kp * P(X_t = C_{nKP-1})`, never a difference of `F` |
| new closed form | `closed_form.quantile_new(q, n, law)`, `cdf_min_new` | 0.02 ms | a one-liner of Eq (1)'s kind that keeps `koff` (below) |
| reference (spectral) | `chain.log_survival_spectral(t, chain)` | ~40 ms | `O(1)` in `t`, so it checks the default where `Lam*t` is 1e10 |

Costs are for a 200-point grid at `nKP=3, n=1e4`.

`chain` also takes `backend="auto"`, which uses uniformization while `Lam*t` is
small and `mp.expm` beyond — needed because large `sigma` pushes median
activation times to `1e6`–`1e35` s, where an `O(Lam*t)` sweep is hopeless.

## How good is the paper's ODE approximation?

Measured against exact `S(t)^n` (see `ode.compare_to_exact`, and
`tests/test_ode.py` which pins all of this):

| | ODE error on the median | Eq (1) error on the median |
|---|---|---|
| `nKP=3, sigma=1` | 1.2e-5 | 3.8e-2 |
| `nKP=10, sigma=1` | 6.5e-6 | 3.9e-1 |
| `nKP=3, sigma=100` | 3.6e-8 | 1.0 (three orders of magnitude) |
| `nKP=10, sigma=100` | 2.2e-15 | 1.0 |

Two things make this sharp rather than anecdotal.

**Its worst case is analytic.** At small `sigma` the recycled activation flux
has no time to redistribute, so the population never depletes and the ODE's
cumulative hazard is `n F(t)` — the Poisson clock — instead of the exact
`-n log(1-F(t))`. The resulting sup-norm error is therefore
`sup_t |exp(-nF) - (1-F)^n| = 2/(e^2 n)`, i.e. `0.2707/n`, independent of
`nKP`; measured `0.27067` at `n = 1e3, 1e4, 1e5`.

**It converges like `1/n`, not `n^(-1/nKP)`.** So the ODE gains a decade of
accuracy per decade of `n` where Eq (1) gains ten per cent — and it is best
exactly where Eq (1) is worst, becoming exact to machine precision at large
`sigma` (the activation time then far exceeds the chain relaxation time, so the
recycled flux redistributes into the conditional quasi-steady state before it
can matter, which is what the exact closure does).

Two caveats worth keeping in view:

* For the *independent-pair* model there is no reason to prefer the ODE: exact
  `S(t)^n` is both more accurate and ~3-25x faster. The ODE earns its place on
  the paper's **second-order** binding, `k1 R (L_T - (R_T - R))`, where the
  pairs are coupled and no product formula exists at all.
* Second-order binding is a genuinely different model, not a reparametrisation.
  At `nKP=10, sigma=10` its median is 10x below the independent-pair one. It
  does, though, explain the paper's initial condition: Eq (14) starts every
  pair unbound, yet Eq (1) needs shape `nKP`, and at `R_T = L_T = 1e4` the
  binding propensity `k1 n^2 = 1e8 /s` makes engagement instantaneous, so the
  unbound start leaves no trace (measured ratio 1.002 against the engaged
  exact median).

`closure="proportional"` implements the exact redistribution instead of
Eq (14a)'s recycling; it reproduces `S(t)^n` to solver tolerance and exists so
that the difference between the two closures is attributable to the closure
rather than to the integrator.

The default engine is **uniformization**, not a matrix exponential: with
`Lam = max` exit rate, `P = I + Q/Lam` is nonnegative, so
`exp(Qt) = e^{-x} sum_m x^m/m! P^m` is a sum of strictly positive terms and `F`
keeps full relative accuracy at any magnitude — verified against `F_rational`
down to `F ~ 1e-67`, where `1 - S(t)` in float64 is pure roundoff. It is also
~500x faster than `mp.expm` and, because the absorbed-mass sequence `a_m` does
not depend on `t`, one cached sweep serves a whole grid of contact durations.

The SSA runs on the **population** process (aggregated occupancy counts across
all `n` pairs as a single CTMC), which is exact and costs one trajectory rather
than `n`. Its event count grows like `n (kp+kd) T_{1,n}`, so it is the
validation layer only — any claim about large `n` or large `sigma` should come
from `survival_min`.

## Established facts, and where they are pinned

| fact | test |
|---|---|
| `F(t) ~ A t^nKP`, `A = kp^nKP/nKP!` | `test_chain.py::test_short_time_prefactor_*` |
| initial condition sets the shape (`nKP` vs `nKP+1`) | `test_chain.py::test_shape_is_nKP_when_engaged_and_nKP_plus_one_when_unbound` |
| `F = A t^nKP [1 - nKP(kp+kd)t/(nKP+1) + O(t^2)]` | `test_chain.py::test_first_correction_coefficient` |
| `err * n^(1/nKP) -> 0.804` (nKP=3), `0.794` (nKP=10) | `test_asymptotics.py::test_*convergence_constant*` |
| Eq (1) median / exact = 0.962 (nKP=3), 0.609 (nKP=10) at `n=1e4` | `test_asymptotics.py::test_eq1_accuracy_at_the_papers_n_self` |
| sigma absent from Eq (1); exact median 71.16 s at sigma=100 | `test_asymptotics.py::test_eq1_fails_badly_at_large_sigma` |
| SSA agrees with `S(t)^n` | `test_ssa.py::test_ssa_samples_the_exact_distribution` (KS on the probability-integral transform) |
| renewal sampler agrees with both `S(t)^n` and the SSA | `test_renewal.py::test_renewal_min_samples_the_exact_distribution`, `::test_renewal_agrees_with_the_jump_by_jump_ssa` |
| `sigma` multiplies the Eq (1) error by `(1+sigma)/2` at every `nKP` | `test_renewal.py::test_sigma_multiplies_the_eq1_error_by_one_plus_sigma_over_two` |
| the error collapses onto `Lambda = (kp+kd) t_Eq1/(nKP+1)` | `test_renewal.py::test_eq1_error_collapses_onto_the_validity_group` |
| exact `1 - S(t)^n` matches the simulation at the `R^2` ceiling in all 25 cells | `test_figures.py::test_exact_law_matches_the_simulation_in_every_cell` |
| `f(t)` is Erlang at `kd = 0`, integrates to 1, and matches `dF/dt` | `test_chain.py::test_single_pair_density_is_erlang_when_it_cannot_dissociate` and neighbours |
| `t f(t)` peaks at `1/e` for an exponential (so peak height reads the shape) | `test_chain.py::test_log_time_density_peaks_at_one_over_e_for_an_exponential` |
| the new closed form's mean is `-Q_T^{-1} 1` exactly (so `tau_c`, `kappa` are derived) | `test_closed_form.py::test_mean_first_passage_matches_the_generator` |
| it reproduces Eq (1) *and* Eq (1)'s exact first correction at short time | `test_closed_form.py::test_working_form_reproduces_eq1_and_its_exact_first_correction` |
| its median is within 6.2% of exact over the 5x5 grid, and beats Eq (1) in every cell | `test_closed_form.py::test_the_new_form_beats_eq1_at_every_quantile_of_the_grid` |
| the working form is monotone and stable where the bounded parent form is not | `test_closed_form.py::test_the_working_form_is_the_numerically_stable_one` |
| uniformization agrees with a spectral evaluation of `Q_T` out to `t = 1e9` | `test_closed_form.py::test_spectral_survival_matches_uniformization` |

`asymptotics.convergence_constant` derives the 0.804 / 0.794 constants in closed
form, `(c/k) (A^{-1} log 2)^{1/k}`, so they are predictions rather than fitted
numbers.

## What breaks Eq (1), and how many knobs there are

Eq (1) contains `nKP`, `kp` and `n`, and nothing else. The model has three rates and
`nKP`, `n`, and a start state; rescaling time removes one rate, so beyond
Eq (1)'s own variables the entire remaining space is two ratios,
`sigma = kd/kp` and `k1/kp`. Only one of them matters, and it matters a lot:

| knob | what it does to the Eq (1) median error | coupled to `nKP`? |
|---|---|---|
| `nKP` | grows it like `t_q/(nKP+1)`, `t_q = (nKP!/n)^(1/nKP)` | -- |
| `n` | shrinks it like `n^(-1/nKP)` | **yes**, the exponent is `1/nKP` |
| `sigma = kd/kp` | multiplies it by `(1+sigma)/2` | **no** |
| `k1/kp` | 1e-2 -> 1e4 moves it from 0.0379 to 0.0314 (17%) | second order |
| all rates x s | nothing: the error is dimensionless | -- |

The reason `sigma` is the clean one is that the correction term factorises.
From `chain.short_time_law`,

```
F(t) = A t^nKP [1 - nKP (kp + kd) t/(nKP+1) + O(t^2)],   A = kp^nKP/nKP!
```

so the relative error of the Eq (1) quantile is

```
Lambda = (kp + kd) t_q/(nKP+1) = (1 + sigma) * [t_q/(nKP+1)]
```

-- a `sigma`-only factor times an `(nKP, n)`-only factor. Equivalently,
`F(t)/(A t^nKP)` depends on `sigma` **only** through the product `(kp+kd) t`:
across `sigma` from 1 to 200 that ratio moves by under 2% at fixed `(kp+kd) t`.
`sigma` does not change the shape of Eq (1)'s validity window, it shrinks the
window by `(1+sigma)`. `nKP` and `n` decide where the activation time sits
inside it.

`figures.figure_eq1_separability` shows both halves: all 63 `(nKP, sigma)`
pairs collapse onto `error = min(Lambda, 1)` (within 2% for `Lambda < 0.05`,
and the shortfall is itself `O(Lambda)`), and dividing out the `sigma = 1`
column leaves `(1+sigma)/2` at every `nKP` until the error saturates against
its own ceiling of 1.

The mechanism is a change of regime, not a loss of precision. Eq (1) is the
first-attempt law: it counts the paths that climb all `nKP` steps without ever
falling off, which requires `t << 1/(kp+kd)`. Once `(1+sigma) t_q` exceeds 1
the typical pair has dissociated and rebound many times, the first-passage time
becomes a geometric number of attempts and therefore approximately
*exponential*, and the extreme statistic inherits shape 1 instead of shape
`nKP`. `figures.effective_weibull_shape` measures it: at `n = 1e4` the slope of
`log(-log S)` falls from 5.4 to 1.0 across `sigma = 1 -> 5` at `nKP = 10`.

### The 5x5 grid

```python
from fram.figures import figure_eq1_grid
figure_eq1_grid("fram_eq1_grid.png")     # ~4.5 min (or pass `cells=`)
```

Rows `nKP = 2, 4, 6, 8, 10`, columns `sigma = 1, 2, 3, 4, 5`, `n = 1e4`, 6000
simulated contacts per cell, each panel simulation vs Eq (1) vs exact on a
shared axis of `t / median_Eq(1)`. Selected cells:

| | `sigma = 1` | `sigma = 3` | `sigma = 5` |
|---|---|---|---|
| `nKP = 2` | x1.01, shape 2.0 | x1.01, shape 2.0 | x1, shape 2.0 |
| `nKP = 6` | x1.23, shape 4.8 | x1.76, shape 2.9 | x6.16, shape 1.0 |
| `nKP = 10` | x1.64, shape 5.4 | x58.5, shape 1.0 | x2.92e+03, shape 1.0 |

(`x` = simulated median / Eq (1) median; `shape` = measured Weibull shape, which
Eq (1) says must be `nKP`.) The top row is the control: 6000 replicates put the
Monte Carlo floor on the median at ~0.8% and on the sup-norm gap at 0.018, and
every cell of that row sits on the floor.

### The same grid without Eq (1), and the R^2 page

```python
from fram.figures import (sim_vs_exact_cells, figure_pdf_grid,
                          figure_sim_vs_exact_grid, figure_r2_page,
                          figure_eq1_grid)
cells = sim_vs_exact_cells((2, 4, 6, 8, 10), (1., 2., 3., 4., 5.))   # ~4.5 min
figure_pdf_grid("fram_pdf_grid.png", cells=cells)
figure_sim_vs_exact_grid("fram_sim_vs_exact_grid.png", cells=cells)
figure_r2_page("fram_r2_page.png", cells=cells)
figure_eq1_grid("fram_eq1_grid.png", cells=cells)
```

All four figures take the same `cells`, so they describe one set of 150 000
simulated contacts rather than four runs that happen to agree.

`figure_sim_vs_exact_grid` drops Eq (1) and gives each panel its own axis in
seconds, because with the approximation gone there is nothing left to hold
fixed across cells -- and the absolute axis says what the normalised one hides:
the median activation time runs from 12 ms at `(2, 1)` to 5.1e3 s at `(10, 5)`,
five orders of magnitude on a grid whose two axes span a factor of 5 each.

`figure_r2_page` is its caption: the generator, the exact law, the definition of
`R^2`, and the 25 values. Both take the *same* `cells`, so the page describes
the panels rather than an independent run that happens to agree.

The page states its own ceiling, which is the part that is easy to get wrong.
Nothing is fitted -- `1 - S(t)^n` is a prediction -- so the residual is pure
sampling noise of known size: the total sum of squares is exactly `(N^2-1)/(12N)`
for midpoint plotting positions and the expected residual sum is `N/(6(N+1))`,
so a model that is *exactly right* scores

```
E[R^2] = 1 - 2N/(N^2 - 1)   =   0.999667   at N = 6000
```

Measured over the grid: min 0.99874, mean 0.99962, max 0.99995, with 16 of 25
cells at or above the ceiling. Two-sided KS per cell: smallest `p = 0.025`,
2 of 25 below 0.05 against 1.25 expected by chance. So the answer the page
gives is "indistinguishable at this sample size", not "0.04% wrong" -- and
without the ceiling written down, 0.9996 reads as a small discrepancy instead
of the best attainable score.

The page also carries the *method*: the generator, the exact CDF and density,
and why the density is drawn as `t f(t)` against `log t`. None of that is on
the grids. A small-multiples grid earns its keep by being scannable, and four
paragraphs of method above one is four paragraphs the eye has to get past
twenty-five times; the legend sits under the axis for the same reason.

### Densities on one scale

The medians across the grid differ by five orders of magnitude, so `f(t)`
against `t` cannot share an axis at all: the `(2, 1)` density peaks at 58 /s and
the `(10, 5)` density at 7e-5 /s, and on one linear scale 24 of the 25 curves
are a flat line on the floor. `t f(t)` against `log t` can:

* its integral over `d log t` is 1 for every parameter set, so all 25 curves
  enclose the same area and their heights compare directly;
* for a Weibull of shape `k` it peaks at exactly `k/e`, so **the height of a
  curve is a reading of its shape parameter**. Measured across the grid:
  peak 0.729 -> 1.98 at `(2, 1)`, 1.728 -> 4.70 at `(6, 1)`, 1.889 -> 5.13 at
  `(10, 1)`, 0.368 -> 1.00 at `(10, 5)`, against `effective_weibull_shape` of
  1.98, 4.80, 5.38 and 1.00. The `(10, 1)` cell is the loose one, at 5%: the
  exact law there is not quite a Weibull, so peak height and interquartile
  slope are measuring slightly different things about the same curve.

So the figure says in one look what the CDF grid says in numbers: down the rows
the peak grows *taller and narrower* -- proofreading sharpens the timer -- and
across the columns it collapses to a low broad memoryless hump three to four
decades to the right. The simulated histogram uses a fixed bin width in
decades, so bin choice does no work either.

The density is exact, not a difference quotient. `C_{nKP-1} --kp--> C_nKP` is
the only way into the absorbing state, so `f(t) = kp * P(X_t = C_{nKP-1})`, and
that occupancy is another positive Poisson sum out of the same uniformization
sweep as `F` and `S` (`c` in `_uniformization_series`, carried at no extra
cost). Differencing `F` would be the obvious route and the wrong one: `f` is
wanted exactly where `F` changes fastest *and* where `F` may be 1e-40. Pinned
against closed-form Erlang at `kd = 0` to 1e-12, against `d/dt cdf_min` to
1e-5, and by `int f dt = 1` to 1e-6.

### Simulating where Gillespie cannot

`ssa.sample_min_activation` costs `n (kp+kd) T_{1,n}` events, which is ~3e8 per
replicate in the bottom-right cell and ~1e20 at the paper's own `sigma = 100`.
`renewal.sample_min_activation_renewal` draws the same law attempt-by-attempt
instead of jump-by-jump: the number of failed attempts is exactly
`Geometric(r^nKP)`, `r = kp/(kp+kd)`, and the holding times inside them merge
into three Gamma draws, so one pair costs O(1) draws however many million
attempts it makes. It touches neither the generator nor mpmath, so it is an
independent check rather than a restatement -- pinned in `tests/test_renewal.py`
against `-Q_T^{-1} 1` for the mean, against `cdf_min` by the probability-integral
transform out to `(nKP=10, sigma=5)`, and against the Gillespie SSA by two-sample
KS wherever Gillespie is still affordable.

## A closed form that keeps `koff`

Eq (1) fails because it is the first-attempt law. The repair is to keep the
attempts: from `C_0` a pair succeeds with probability `sigma = r^nKP`,
`r = kp/mu`, `mu = kp + koff`, in a time `Gamma(nKP, mu)`, and otherwise burns
a failed cycle of mean length `tau_c` and starts again. Replacing the
compound-geometric sum of those cycles by a single exponential of the same
mean -- its exact `sigma -> 0` limit -- gives

```
mu    = kp + koff = kp (1 + omega)                omega = koff/kp
sigma = (kp/mu)^nKP                               per-attempt success prob.
tau_c = 1/k1 + ((1-sigma)/koff - nKP sigma/mu)/(1-sigma)
kappa = sigma/tau_c                               effective trigger rate

F1(t) = sigma [ P(nKP, mu t)
                + ( t P(nKP, mu t) - (nKP/mu) P(nKP+1, mu t) ) / tau_c ]
P(T_1n <= t) = 1 - (1 - F1(t))^n
```

with `P` the regularized lower incomplete gamma. Nothing in it is fitted:
`tau_c` is the conditional mean number of bound-state holds in a failed
attempt, `kappa` is fixed by making the mean exact, and
`nKP/mu + tau_c (1-sigma)/sigma` reproduces `-Q_T^{-1} 1` from the generator to
1e-12 (`test_closed_form.py::test_mean_first_passage_matches_the_generator`).

```python
from fram import RetryLaw, quantile_new, cdf_min_new
law = RetryLaw(n_kp=10, kp=1.0, koff=5.0, k1=1.0)
quantile_new(0.5, n=1e4, law=law)      # 5029.7 s   (exact: 5030.1, Eq (1): 1.74)
```

Measured against the exact law on the note's 5x5 grid (`n = 1e4`, `kp = k1 = 1`),
relative error of the median, `t_new/t_exact - 1` (the note's `E` with the sign
flipped):

| | `omega=1` | `omega=2` | `omega=3` | `omega=4` | `omega=5` |
|---|---|---|---|---|---|
| `nKP=2` | -0.12% | -0.14% | -0.15% | -0.16% | -0.17% |
| `nKP=4` | -0.62% | -0.85% | -1.05% | -1.24% | -1.47% |
| `nKP=6` | -1.06% | -1.78% | -2.96% | **-6.12%** | -4.33% |
| `nKP=8` | -1.46% | -3.58% | -3.59% | -0.58% | -0.12% |
| `nKP=10` | -1.92% | -4.01% | -0.26% | -0.03% | -0.01% |

Eq (1) over the same cells runs -0.78% to -99.97%, i.e. -0.78% in the top-left
and a factor of 2900 in the bottom-right. The new form is worst at
`(nKP=6, omega=4)`, the crossover where neither the Erlang first attempt nor
the memoryless tail dominates, and best at the two ends. Over a 1080-cell sweep in
`(nKP, omega, k1/kp, n)` -- `omega` to 28, the self-ligand edge of the note's
Table IV, `k1/kp` over four decades, `n` from 1e2 to 1e5 -- the worst error at
any of the 5%/25%/50%/75%/95% quantiles is 12.5%, at `n = 1e2` where the
extreme-value target `-log(1-q)/n` is largest; Eq (1) reaches 1.000 there.

Three things are worth knowing about the form itself.

**It contains Eq (1) and Eq (1)'s exact first correction.** Its short-time
expansion is `A t^nKP [1 - nKP mu t/(nKP+1) + t/((nKP+1) tau_c) + O(t^2)]` with
`A = kp^nKP/nKP!`. The first two terms are `chain.short_time_law` term for
term. The third is spurious -- a retry cannot contribute before
`O(t^(nKP+2))` in the exact law -- and is smaller than the genuine correction
by exactly `mu tau_c >= 1 + mu/k1 > 1`, so it can never overturn it. That
ratio is `1/(mu tau_c)` identically, verified to 2e-16 at 270 parameter
combinations.

**Its far field is memoryless, and exactly so.** At large `omega` the median of
`T_{1,n}` is `nKP/mu + (target - sigma)/kappa` with
`target = 1 - (1-q)^(1/n)`, to machine precision. Dropping `sigma` and
linearising `target` as `ln2/n` gives the familiar
`ln2/(n kappa) + nKP/mu`, which is what the residual `1 - ln2/(2n)` in the
summary table is. The effective Weibull shape goes to 1, matching the exact
law to three decimals.

**The bounded form is the numerically bad one.** `F1_parent` (the unexpanded
Erlang-plus-exponential convolution, `closed_form.F1_parent`) is bounded and
`-> 1`, but it is a difference of two quantities both close to
`P(nKP, mu t)`, so in float64 it keeps only `log10(F1/P)` digits: 4e-6
relative at `nKP=10, omega=28`, where the working form holds 1e-15. The
working form is `sigma` times a sum of two nonnegative terms -- its derivative
is `sigma Erlang(nKP, mu; t) + kappa P(nKP, mu t)`, so monotonicity is a
theorem rather than a measurement -- and it is unbounded, crossing 1 at the
mean first-passage time. Since the extreme statistic lives at
`F1 ~ -log(1-q)/n`, decades below that, the unbounded one is the one to
compute with.

### The brackets, and the one that does not bracket

`quantile_new` root-finds on `[t_W, t_U]`:

* `t_W = (nKP! L/(n kp^nKP))^(1/nKP)`, `L = -ln(1-q)` -- Eq (1)'s own
  quantile, and a **lower** bound, because the net `O(t^(nKP+1))` coefficient
  of `F1` is negative;
* `t_U = nKP/mu + tau_c * target/sigma` -- an **upper** bound, because
  `F1(t) >= kappa (t - nKP/mu)` for every `t`.

The natural-looking `t_E = nKP/mu + tau_c (L/(n sigma) - 1)`, the quantile of
the far field with the `sigma` dropped, is **not** a bound: it undershoots by
about `tau_c` and goes negative outright whenever the first attempt dominates
(`L/n < sigma`, e.g. `nKP=2, omega=1`, where `t_E = -0.67`). It overtakes even
`t_U` once `sigma < L^2/(2n^2)`, so `[min(t_W,t_E), max(t_W,t_E)]` is a valid
bracket deep in the retry regime and not in between -- 22 of 45 grid
cell/quantile combinations at `n = 1e4`. `quantile_new(..., literal_bracket=True)`
uses it anyway and raises where it fails; `bracket_straddles` reports which
regime a cell is in.

### Validating it

```
python validate_new_closed_form.py            # ~50 min
python validate_new_closed_form.py --quick    # ~1 min, coarser
```

writes `results/new_closed_form.md` (the summary, with the worst-case tables),
`results/new_closed_form.csv` and four companion CSVs, and
`results/new_closed_form_grid.png` -- the Fig 1 layout with the new curve
overlaid on the simulation, the exact law and Eq (1), from the same 150 000
simulated contacts. Six checks, in the order they have to be believed: the
exact reference against a spectral evaluation of the sub-generator (agreeing
to 9e-11 from `t = 1e-3` to `1e9`), the 5x5 grid, 6000 attempt-level contacts
per cell, the 1080-cell stress sweep, both limits, and the size of the one
spurious term.

## The classifier (Figs 2 and 3)

```python
from fram.classifier import Environment, optimal_contact_duration, auc

env = Environment.from_paper(nKP=10, sigma=15.0, rho_ag=0.5)
opt = optimal_contact_duration(env)      # tau* = 1.7e4 s, Gamma(tau*) = 0.999
auc(env)                                 # 0.9997
```

`Gamma = P(TP) + P(TN)` with the four outcomes joint over the APC condition and
the activation time, so `Gamma(0) = 1 - rho_ag`, `Gamma(inf) = rho_ag`, and the
baseline (the best constant strategy) is `max(rho_ag, 1 - rho_ag)` -- the
paper's Eq (4). Where no contact duration beats it, `Optimum.at_baseline` is
set and `tau*` is reported as 0 or infinity, which is what makes the published
curves stop short.

### Reproducing the published figures

`figures.figure2()` and `figures.figure3()` redraw both, and
`figures.diff_published_fig2()` compares against landmarks recovered from the
published PNGs by colour matching:

| nKP | sigma | published Gamma(tau*) | ours | published tau* | ours |
|---|---|---|---|---|---|
| 3 | 1e4 | 0.9992 | 0.9971 | 87 | 102 |
| 3 | 1e3 | 0.8555 | 0.8339 | 33.1 | 31.5 |
| 3 | 1e2 | 0.5076 | 0.5045 | 2.13 | 2.05 |
| 10 | 15 | 0.9992 | 0.9987 | 1.36e4 | 1.67e4 |
| 10 | 10 | 0.9789 | 0.9775 | 9913 | 9917 |
| 10 | 5 | 0.5875 | 0.5846 | 921 | 955 |

Two parameter readings had to be pinned down to get there, and both are
ambiguities in the paper rather than choices:

* **Two different binding rates.** Table 1 lists one `k1 = 1.0 s^-1`, but the
  self population is driven by the ODE's bimolecular term `k1 R L`, whose
  per-pair rate at full receptor availability is `k1 R_T = 1e4 /s`, while the
  single agonist is driven by the Sec 1.1 formula at `k1 = 1.0 /s`. Using
  either value for both misses the published `tau*` by 1.5-2x; using both
  reproduces it. `Environment.from_paper` sets `k1_self = n_self`, `k1_ag = 1`.
* **The initial condition stops mattering.** Eq (1)'s shape of `nKP` needs
  engaged pairs, while the ODE starts them unbound. At `k1_self = 1e4` binding
  takes 1e-4 s and the two initial conditions agree to 0.1%, so the tension is
  moot at the paper's own parameters. At `k1 = 1` they differ by ~3x.

### The independence assumption

The paper factorises `P(no activation | xi=1) = P(T_ag > tau) P(T_self > tau)`,
noting it ignores "interactions, or interference, between self and agonist
ligands". In the independent-pair model that factorisation is **exact** --
distinct pairs, distinct first-passage times -- so the assumption costs nothing
there. It costs something only through the mechanism the paper names: with
`R_T = L_T = n_self` the self antigen occupy most of the receptors, so an
agonist binds at `k1 R(t)/R_T`. `independent=False` restores that competition,
taking the free-receptor fraction from the Sec 1.2 ODE. Measured fractions are
0.10 (sigma=1e2), 0.27 (1e3), 0.62 (1e4) at nKP=3 -- so the interference is
worst exactly where discrimination is hardest, and it can only lower accuracy.

## Futile reactions (Fig 4)

```python
from fram import Chain, mean_futile_reactions
mean_futile_reactions(tau=100.0, n=1e4, chain=Chain(3))    # 4.00e5
```

S1 Sec 5 Eq (15) augments the conditioned ODE with

```
dn_e/dt = k_1 (C_1 + 2 C_2 + ... + (nKP-2) C_{nKP-2}) + (k_1 + kp)(nKP-1) C_{nKP-1}
```

counting `i` wasted phosphorylations whenever a pair leaves `C_i`. It is linear
in the same states as Eq (14), so it rides along in the matrix exponential at
no extra cost, and `terminate_on_activation=True` selects Eq (16), which gates
the count on `1 - A` so counting stops at activation (the paper's figures use
Eq (15); Eq (16) is necessarily smaller and saturates).

Two checks worth knowing about, both pinned in `tests/test_futile.py`:

* At `nKP=3` with all rates 1, the quasi-steady occupancies are
  `C_0, C_1, C_2 = 2667, 1333, 667`, so Eq (15) gives
  `1333 + 4(667) = 4000` reactions per second — which is what the code returns
  to 2%.
* `kd * r = kd kp/(kp+kd) -> kp` as sigma grows, so the whole population wastes
  roughly one phosphorylation per `1/kp` **regardless of sigma**. That is why
  Fig 4's `n_e` axis barely shifts across its three columns.

Past relaxation `n_e` is exactly affine in tau, and `figures.figure4` inverts
that to pick tau from the `n_e` axis rather than from the accuracy curves —
which span tau to 1e30 s where the axis stops at ~1e6 s. Without it every point
falls back to 60-digit mpmath and one panel takes 187 s instead of 44 s.

`figure4` reproduces the published staircase: the leading edge of the
high-accuracy band moves from `n_e ~ 5e5` at `nKP=5` to `~5e8` at `nKP=15`
(published: ~2e5 to ~5e8). It uses a **sequential single-hue** ramp rather than
the published `jet`; a rainbow invents category boundaries in a magnitude and
is unreadable under colour-vision deficiency, so the data is reproduced and the
encoding is not.

## Figures

```python
from fram.figures import figure_weibull_check
figure_weibull_check("fram_weibull_check.png")      # + .csv table view
```

Every figure function returns the records it plotted and writes them as CSV
next to the image.

## Citation

If you use fram, please cite the archived release:

> T. Mourier, *fram: kinetic proofreading with extreme statistics*, version
> 1.0.0, Zenodo (2026), doi:[10.5281/zenodo.22999081](https://doi.org/10.5281/zenodo.22999081).

```bibtex
@misc{mourier_fram_2026,
  author       = {Mourier, Thomas},
  title        = {{fram}: kinetic proofreading with extreme statistics},
  howpublished = {Zenodo},
  note         = {Version 1.0.0},
  year         = {2026},
  doi          = {10.5281/zenodo.22999081},
  url          = {https://doi.org/10.5281/zenodo.22999081}
}
```

`CITATION.cff` in the source archive carries the same metadata.

## License

MIT; see `LICENSE`.
