Skip to content

assign_compartments gap-fill: reliable flux-based fill, not cobra's MILP - #82

Open
edkerk wants to merge 2 commits into
developfrom
feat/flux-based-gapfill
Open

assign_compartments gap-fill: reliable flux-based fill, not cobra's MILP#82
edkerk wants to merge 2 commits into
developfrom
feat/flux-based-gapfill

Conversation

@edkerk

@edkerk edkerk commented Jul 19, 2026

Copy link
Copy Markdown
Member

The universal-DB gap-fill in assign_compartments used cobra.flux_analysis.gapfill, whose indicator MILP is unreliable at genome scale. This replaces it with a plain flux-based fill.

The problem (measured)

cobra's GapFiller MILP fails to find a valid gap-fill in the majority of genome-scale cases even when the exact reaction that restores growth is present in the universal:

  • Knockout-recovery (remove an essential reaction, gap-fill from a universal that contains it): cobra recovers only 45%. On failures the MILP doesn't select the restoring reaction, and its own validate() correctly rejects the broken incumbent and raises.
  • It is not a validation-tolerance issue (skipping validate() gives identical results) and not a universal-size issue (a minimal 1-reaction universal fails identically) — the MILP itself is the ceiling.
  • On realistic incomplete drafts (12% of internal reactions dropped): cobra restores growth 0/12.

The fix

On a working copy of the model: add the universal candidates, hold biomass at the growth floor, run pFBA, and keep the added reactions that carry flux. A plain LP — it cannot have the MILP's failure mode, so a returned set always actually restores growth. The caller still re-certifies with a real FBA, so no false certificate is possible.

The set is flux-parsimonious (pFBA minimises total flux) and returned sorted, so it does not depend on which co-optimal vertex the solver lands on. It is not guaranteed reaction-count-minimal; the caller re-certifies regardless.

Result (measured, fixed reproducible knockout sample)

knockout-recovery realistic drafts
cobra gapfill (before) 45% 0 / 12
pFBA fill (this PR) 60 / 60 (exact reaction each time) 12 / 12

Namespace safety

Candidates are matched to the model by metabolite id (as cobra's gapfill required) — the universal must share the draft's metabolite namespace. A candidate whose ids don't resolve becomes a dead-end that can't carry flux and is left out. Previously that produced a silent []; now, when most candidates fail to resolve, _gapfill warns, so "found nothing" is distinguishable from "wrong namespace".

Tests

Existing gap-fill tests pass unchanged; new direct _gapfill tests cover the growth-restoring set, the empty-universal case, and the namespace-mismatch warning. Full suite green; ruff and mypy clean.

Note

The assignment-ablation study in #75 documents the old 45% ceiling. This fix supersedes that framing — its gap-fill section should be updated to the reliable pFBA fill when #75 merges.

edkerk added 2 commits July 19, 2026 09:02
…atch

Harden the flux-based _gapfill (universal-DB fill in localization.certify):
operate on a model copy so the caller's model is never mutated; add candidates
in one batch (per-reaction adds are super-linear at scale); return the
flux-carrying set sorted, so the result does not depend on which co-optimal
vertex the solver picked; and warn when most universal candidates share no
metabolite id with the model, so a namespace mismatch is distinguishable from
a genuine empty fill instead of returning [] silently.

Knockout-recovery 60/60 (exact reaction each time) and 12/12 on realistic
incomplete drafts, vs cobra's 45% and 0/12.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant