From 9a1656f84e0a5cb89332cd0d0acaedd9a87a6ef7 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 01:36:36 +0200
Subject: [PATCH 1/7] feat: random sampling of the conventional model

Two passes per condition, one with formate left out of the measured set
so the sampler can say which byproducts the model makes unprompted, and
the ATP and redox budget over the sampled means.

yeast-GEM has merged glucokinase (r_4235) into hexokinase, so the
glycolytic ATP expenditure is now read from r_0534 alone; the accounting
takes the reaction sets as arguments, which is what lets it be checked
against the published summary on that summary's own reaction set.
---
 README.md                    |  15 +++
 src/overflow/atp.py          | 140 ++++++++++++++++++++++
 src/overflow/run_sampling.py | 149 +++++++++++++++++++++++
 src/overflow/sampling.py     | 221 +++++++++++++++++++++++++++++++++++
 tests/test_sampling.py       | 156 +++++++++++++++++++++++++
 5 files changed, 681 insertions(+)
 create mode 100644 src/overflow/atp.py
 create mode 100644 src/overflow/run_sampling.py
 create mode 100644 src/overflow/sampling.py
 create mode 100644 tests/test_sampling.py

diff --git a/README.md b/README.md
index 49a72d8..16d0894 100644
--- a/README.md
+++ b/README.md
@@ -110,6 +110,21 @@ Writes, per enzyme and condition, how much of it the model uses and what
 fraction of what was available that is, plus the capacity usage of the
 annotated systems and the two figures over them.
 
+## Random sampling
+
+```bash
+python -m overflow.run_sampling --procs 16
+```
+
+Samples the conventional model under the measured rates, with and without
+formate in the measured set, and writes the sampled means and standard
+deviations, the byproducts the model secretes when it is not told to, and the
+ATP and redox budget.
+
+Sampling draws vertices by maximising small random objectives, as the published
+analysis did; means over such draws are vertex means rather than an average over
+the interior of the flux space.
+
 ## Tests
 
 ```bash
diff --git a/src/overflow/atp.py b/src/overflow/atp.py
new file mode 100644
index 0000000..49b6487
--- /dev/null
+++ b/src/overflow/atp.py
@@ -0,0 +1,140 @@
+"""Where the ATP comes from and what turns over the redox cofactors.
+
+The sampled flux distributions are summarised into a budget: how much
+ATP respiration makes, how much substrate-level phosphorylation makes,
+and how much growth and maintenance spend. Expressed per unit of glucose
+and per unit of growth, so conditions with different uptake rates can be
+compared.
+"""
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Mapping, Optional, Sequence
+
+import numpy as np
+import pandas as pd
+
+from overflow.config import BIO_RXN, NGAM_RXN
+
+if TYPE_CHECKING:  # pragma: no cover - typing only
+    import cobra
+
+#: ATP made by the respiratory chain and by substrate-level
+#: phosphorylation in the TCA cycle.
+RESPIRATION = ("r_0226", "r_1022")
+
+#: ATP made in glycolysis, and the steps that spend it. yeast-GEM merged
+#: glucokinase (r_4235 in 8.3.4) into hexokinase, so r_0534 now carries
+#: both and is counted once.
+GLYCOLYSIS_PRODUCING = ("r_0892", "r_0962")
+GLYCOLYSIS_CONSUMING = ("r_0886", "r_0534")
+
+GLUCOSE_TRANSPORT = "r_1166"
+
+#: Reactions whose flux the published summary reports directly.
+TURNOVER = {
+    "rPDH": "r_0961",
+    "rIDH": "r_0658",
+    "rMDHc": "r_0714",
+    "rMDHm": "r_0713",
+    "rNDE": "r_0770",
+}
+
+
+def _sum(means: Mapping[str, float], reactions: Sequence[str]) -> float:
+    return float(sum(float(means.get(r, 0.0)) for r in reactions))
+
+
+def atp_budget(
+    means: Mapping[str, float],
+    gam: float,
+    growth_rate: float,
+    glucose_rate: Optional[float] = None,
+    respiration: Sequence[str] = RESPIRATION,
+    glycolysis_producing: Sequence[str] = GLYCOLYSIS_PRODUCING,
+    glycolysis_consuming: Sequence[str] = GLYCOLYSIS_CONSUMING,
+    glucose_transport: str = GLUCOSE_TRANSPORT,
+    turnover: Mapping[str, str] = TURNOVER,
+) -> pd.Series:
+    """The ATP budget implied by one flux distribution.
+
+    ``gam`` is the growth-associated ATP cost written into the biomass
+    reaction. Rates are mmol ATP/gDW/h; yields are per mmol glucose and
+    per unit of growth.
+    """
+    if glucose_rate is None:
+        glucose_rate = float(means.get(glucose_transport, 0.0))
+
+    respiration_rate = _sum(means, respiration)
+    glycolysis = _sum(means, glycolysis_producing) - _sum(means, glycolysis_consuming)
+    growth_cost = gam * float(means.get(BIO_RXN, 0.0))
+    maintenance = float(means.get(NGAM_RXN, 0.0))
+    metabolism = respiration_rate + glycolysis - growth_cost - maintenance
+
+    rows: dict[str, float] = {"rGlu": glucose_rate}
+    for name, rate in (
+        ("ETC", respiration_rate),
+        ("glycolysis", glycolysis),
+        ("GAEC", growth_cost),
+        ("NGAM", maintenance),
+        ("Metabolism", metabolism),
+    ):
+        rows[f"{name}_rATP"] = rate
+        rows[f"{name}_YATP_glu"] = rate / glucose_rate if glucose_rate else np.nan
+        rows[f"{name}_YATP_mu"] = rate / growth_rate if growth_rate else np.nan
+
+    spent = growth_cost + maintenance + metabolism
+    rows["GAEC+NGAM+Metabolism_rATP"] = spent
+    rows["GAEC+NGAM+Metabolism_YATP_glu"] = spent / glucose_rate if glucose_rate else np.nan
+    rows["GAEC+NGAM+Metabolism_YATP_mu"] = spent / growth_rate if growth_rate else np.nan
+
+    for name, reaction in turnover.items():
+        rows[name] = float(means.get(reaction, 0.0))
+
+    return pd.Series(rows)
+
+
+def cofactor_turnover(
+    model: "cobra.Model",
+    means: Mapping[str, float],
+    metabolite_name: str,
+    label: str,
+) -> pd.Series:
+    """Production of one cofactor per compartment.
+
+    Only production counts: consumption is its mirror image, and summing
+    both would report zero everywhere.
+    """
+    metabolites = [m for m in model.metabolites if m.name == metabolite_name]
+    rows: dict[str, float] = {}
+    total = 0.0
+    for metabolite in metabolites:
+        produced = 0.0
+        for reaction in metabolite.reactions:
+            flux = float(means.get(reaction.id, 0.0))
+            rate = reaction.metabolites[metabolite] * flux
+            if rate > 0:
+                produced += rate
+        if produced == 0:
+            continue
+        rows[f"{label}[{metabolite.compartment}]"] = produced
+        total += produced
+    rows[f"{label}[tot]"] = total
+    return pd.Series(rows)
+
+
+def selected_fluxes(
+    model: "cobra.Model",
+    means: Mapping[str, float],
+    gam: float,
+    growth_rate: float,
+    polymerization: Optional[Mapping[str, float]] = None,
+) -> pd.Series:
+    """The published summary of one condition's sampled fluxes."""
+    parts = [
+        atp_budget(means, gam, growth_rate),
+        cofactor_turnover(model, means, "NAD", "rNAD"),
+        cofactor_turnover(model, means, "NADP(+)", "rNADP"),
+    ]
+    if polymerization:
+        parts.append(pd.Series({f"GAECpol_{k}": v for k, v in polymerization.items()}))
+    return pd.concat(parts)
diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
new file mode 100644
index 0000000..a90482d
--- /dev/null
+++ b/src/overflow/run_sampling.py
@@ -0,0 +1,149 @@
+"""Random sampling of every condition.
+
+Two passes per condition: one with formate excluded from the measured
+set, which asks the model which byproducts it would make if it were not
+told, and one with every measured rate applied.
+
+``python -m overflow.run_sampling``.
+"""
+from __future__ import annotations
+
+import argparse
+import time
+from pathlib import Path
+from typing import Optional
+
+import pandas as pd
+
+from overflow.adapter import load_gem
+from overflow.atp import selected_fluxes
+from overflow.biomass import GAM_NO_POLYMERIZATION, current_gam
+from overflow.config import (
+    BIO_RXN,
+    CONDITION_ORDER,
+    RESULTS_DIR,
+    load_conditions,
+)
+from overflow.sampling import (
+    polymerization_breakdown,
+    prepare_model,
+    sample_condition,
+)
+
+#: Mean flux above which a byproduct counts as predicted [mmol/gDW/h].
+DETECTION = 0.001
+
+
+def run(
+    conditions: Optional[list[str]] = None,
+    n_samples: int = 5000,
+    seed: Optional[int] = 0,
+    n_proc: Optional[int] = None,
+    verbose: bool = True,
+) -> dict[str, dict]:
+    """Sample each condition twice and summarise."""
+    all_conditions = load_conditions()
+    results: dict[str, dict] = {}
+    good_free: Optional[list[str]] = None
+    good_full: Optional[list[str]] = None
+
+    for name in conditions or CONDITION_ORDER:
+        condition = all_conditions[name]
+        started = time.time()
+
+        model = load_gem()
+        prepare_model(model, condition, gam_base=GAM_NO_POLYMERIZATION)
+        gam = current_gam(model)
+        polymerization = polymerization_breakdown(model)
+
+        free, good_free = sample_condition(
+            model.copy(), condition, n_samples=n_samples, include_formate=False,
+            seed=seed, good_reactions=good_free, n_proc=n_proc,
+        )
+        full, good_full = sample_condition(
+            model, condition, n_samples=n_samples, include_formate=True,
+            seed=seed, good_reactions=good_full, n_proc=n_proc,
+        )
+
+        summary = selected_fluxes(
+            model, full.means, gam=gam, growth_rate=condition.d_rate,
+            polymerization=polymerization,
+        )
+        results[name] = {
+            "free": free, "full": full, "summary": summary,
+            "gam": gam, "model": model,
+        }
+        if verbose:
+            print(
+                f"[{name}] {full.n_samples} samples, GAM {gam:.2f}, "
+                f"maintenance up to {full.max_maintenance:.2f}; "
+                f"growth {full.means.get(BIO_RXN, float('nan')):.5f} "
+                f"[{time.time() - started:.0f}s]",
+                flush=True,
+            )
+    return results
+
+
+def alternative_exchanges(results: dict[str, dict], model, threshold: float = DETECTION) -> pd.DataFrame:
+    """Byproducts the model secretes when it is not told to.
+
+    Taken from the pass with formate left out of the measured set.
+    """
+    exchanges = [r for r in model.reactions if len(r.metabolites) == 1]
+    rows = []
+    for reaction in exchanges:
+        means = {
+            name: float(result["free"].means.get(reaction.id, 0.0))
+            for name, result in results.items()
+        }
+        if max(means.values()) <= threshold:
+            continue
+        rows.append({"rxnID": reaction.id, "rxnName": reaction.name, **means})
+    return pd.DataFrame(rows).sort_values("rxnID")
+
+
+def main(argv: Optional[list[str]] = None) -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("conditions", nargs="*", default=list(CONDITION_ORDER))
+    parser.add_argument("--results-dir", type=Path, default=RESULTS_DIR)
+    parser.add_argument("--samples", type=int, default=5000)
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--procs", type=int, default=None)
+    parser.add_argument("--solver")
+    args = parser.parse_args(argv)
+
+    if args.solver:
+        import cobra
+
+        cobra.Configuration().solver = args.solver
+
+    results = run(args.conditions, n_samples=args.samples, seed=args.seed, n_proc=args.procs)
+
+    out = args.results_dir / "randomSampling"
+    out.mkdir(parents=True, exist_ok=True)
+
+    any_model = next(iter(results.values()))["model"]
+    alternative_exchanges(results, any_model).to_csv(
+        out / "altExchangeFlux.tsv", sep="\t", index=False
+    )
+
+    wide = pd.DataFrame({"rxnID": next(iter(results.values()))["full"].reactions})
+    wide["rxnName"] = [any_model.reactions.get_by_id(r).name for r in wide["rxnID"]]
+    for name, result in results.items():
+        frame = result["full"].to_frame().set_index("rxnID")
+        wide[f"{name}_AVERAGE"] = wide["rxnID"].map(frame["mean"])
+        wide[f"{name}_STDEV"] = wide["rxnID"].map(frame["sd"])
+        wide[f"{name}_STDEV_percent"] = wide["rxnID"].map(frame["sd_percent"])
+    wide.to_csv(out / "allFluxes.tsv", sep="\t", index=False)
+
+    summary = pd.DataFrame({name: r["summary"] for name, r in results.items()})
+    summary.index.name = "Row"
+    summary.to_csv(out / "selectedFluxes.tsv", sep="\t")
+
+    print()
+    print(summary.head(24).round(4).to_string())
+    return 0
+
+
+if __name__ == "__main__":  # pragma: no cover
+    raise SystemExit(main())
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
new file mode 100644
index 0000000..304f21e
--- /dev/null
+++ b/src/overflow/sampling.py
@@ -0,0 +1,221 @@
+"""Random sampling of the conventional model under the measured rates.
+
+No enzyme constraints here: the question is what the measured exchange
+rates alone imply about the flux distribution, so the model is the plain
+yeast-GEM with its biomass rescaled to the condition and every measured
+rate held inside a narrow band around its measurement.
+
+Sampling draws vertices by maximising small random objectives, which is
+what the published analysis did. It is not a uniform draw from the
+interior; means over such draws are vertex means.
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING, Mapping, Optional
+
+import numpy as np
+import pandas as pd
+from raven_toolbox.analysis.sampling import random_sampling
+
+from overflow.biomass import (
+    GAM_NO_POLYMERIZATION,
+    POLYMERIZATION_COST,
+    YEAST_BIOMASS,
+    apply_protein_content,
+)
+from overflow.config import BIO_RXN, GROWTH_RXN, NGAM_RXN, Condition
+
+if TYPE_CHECKING:  # pragma: no cover - typing only
+    import cobra
+
+#: Exchange reactions held at their measurement, and the sign convention
+#: each is measured with: -1 for an uptake, +1 for a secretion.
+MEASURED_EXCHANGES = {
+    "r_1714": ("Glucose", -1),
+    "r_1992": ("Oxygen", -1),
+    "r_1672": ("CO2", +1),
+    "r_1808": ("Glycerol", +1),
+    "r_1634": ("Acetate", +1),
+    "r_1761": ("Ethanol", +1),
+    "r_1793": ("Formate", +1),
+    GROWTH_RXN: ("D", +1),
+}
+
+#: Formate dehydrogenase, blocked so that formate leaves the cell rather
+#: than being oxidised to CO2 inside it.
+FORMATE_DEHYDROGENASE = "r_0445"
+
+FORMATE_EXCHANGE = "r_1793"
+
+
+def measured_rate(condition: Condition, field_name: str) -> float:
+    """One measured rate, as a positive magnitude."""
+    if field_name == "D":
+        return condition.d_rate
+    if field_name == "Glucose":
+        return condition.glucose
+    if field_name == "Oxygen":
+        return condition.oxygen
+    if field_name == "CO2":
+        return condition.co2
+    return condition.byproducts.get(field_name, 0.0)
+
+
+def sampling_bounds(
+    condition: Condition,
+    tolerance: float = 0.05,
+    minimum_secretion: float = 0.01,
+    include_formate: bool = True,
+) -> dict[str, tuple[float, float]]:
+    """Bounds holding every measured rate within ``tolerance`` of itself.
+
+    A byproduct that was not detected gets a small allowance rather than
+    being blocked, so that the sampler can say whether the model wants to
+    make it.
+    """
+    bounds: dict[str, tuple[float, float]] = {}
+    for reaction_id, (field_name, sign) in MEASURED_EXCHANGES.items():
+        if reaction_id == FORMATE_EXCHANGE and not include_formate:
+            continue
+        value = measured_rate(condition, field_name)
+        low = (1 - tolerance) * value
+        high = (1 + tolerance) * value
+        if sign < 0:
+            low, high = -high, -low
+        else:
+            high = high if high != 0 else minimum_secretion
+        bounds[reaction_id] = (low, high)
+    return bounds
+
+
+def prepare_model(
+    model: "cobra.Model",
+    condition: Condition,
+    gam_base: float = GAM_NO_POLYMERIZATION,
+) -> dict[str, float]:
+    """Set the condition's biomass and open maintenance.
+
+    Biomass is rescaled to the measured protein content and maintenance
+    rebuilt from it, which is what this analysis fixes the ATP budget
+    against. Unlike the enzyme-constrained models there is no protein
+    budget here to make that unaffordable.
+    """
+    composition = apply_protein_content(
+        model, condition.p_tot, gam="polymerization", gam_base=gam_base,
+        scale_protein=True,
+    )
+    model.reactions.get_by_id(FORMATE_DEHYDROGENASE).bounds = (0.0, 0.0)
+    model.reactions.get_by_id(NGAM_RXN).bounds = (0.0, 1000.0)
+    return composition
+
+
+def apply_bounds(model: "cobra.Model", bounds: Mapping[str, tuple[float, float]]) -> None:
+    for reaction_id, (low, high) in bounds.items():
+        model.reactions.get_by_id(reaction_id).bounds = (low, high)
+
+
+def constrain_maintenance(model: "cobra.Model", fraction: float = 0.95) -> float:
+    """Hold maintenance near the most the condition can support.
+
+    Everything else is fixed by the measured rates, so whatever ATP the
+    model does not spend on growth it spends here; leaving maintenance
+    free would let it dissipate that instead through whichever internal
+    cycle happens to be cheapest.
+    """
+    maintenance = model.reactions.get_by_id(NGAM_RXN)
+    maintenance.bounds = (0.0, 1000.0)
+    model.objective = NGAM_RXN
+    model.objective.direction = "max"
+    solution = model.optimize()
+    if solution.status != "optimal":
+        raise RuntimeError("no feasible solution under the measured rates")
+    highest = float(solution.fluxes[NGAM_RXN])
+    maintenance.lower_bound = fraction * highest
+    return highest
+
+
+def polymerization_breakdown(model: "cobra.Model", gam_base: float = GAM_NO_POLYMERIZATION) -> dict[str, float]:
+    """ATP spent polymerising each component of one gram of biomass."""
+    from raven_toolbox.biomass.scale import sum_biomass
+
+    fractions = sum_biomass(model, YEAST_BIOMASS)
+    costs = {
+        component: cost * fractions.get(component, 0.0)
+        for component, cost in POLYMERIZATION_COST.items()
+    }
+    costs["total"] = sum(costs.values())
+    costs["GAEC_noPol"] = gam_base
+    costs["GAEC_total"] = costs["total"] + gam_base
+    return costs
+
+
+@dataclass
+class SamplingResult:
+    """Sampled fluxes for one condition and one formate setting."""
+
+    condition: str
+    include_formate: bool
+    reactions: list[str] = field(default_factory=list)
+    mean: np.ndarray = field(default_factory=lambda: np.empty(0))
+    sd: np.ndarray = field(default_factory=lambda: np.empty(0))
+    n_samples: int = 0
+    max_maintenance: float = 0.0
+
+    def to_frame(self) -> pd.DataFrame:
+        with np.errstate(divide="ignore", invalid="ignore"):
+            relative = np.where(self.mean != 0, self.sd / self.mean * 100, np.nan)
+        return pd.DataFrame(
+            {"rxnID": self.reactions, "mean": self.mean, "sd": self.sd,
+             "sd_percent": relative}
+        )
+
+    @property
+    def means(self) -> pd.Series:
+        return pd.Series(self.mean, index=self.reactions)
+
+
+def sample_condition(
+    model: "cobra.Model",
+    condition: Condition,
+    n_samples: int = 5000,
+    include_formate: bool = True,
+    tolerance: float = 0.05,
+    seed: Optional[int] = None,
+    good_reactions: Optional[list[str]] = None,
+    n_proc: Optional[int] = None,
+    min_flux: Optional[bool] = None,
+) -> tuple[SamplingResult, list[str]]:
+    """Sample one condition, returning the summary and the good reactions.
+
+    ``good_reactions`` is the loop-free set the sampler may use as random
+    objectives; it depends only on the network, so it is computed once
+    and handed back for reuse.
+    """
+    apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate))
+    highest = constrain_maintenance(model)
+
+    sampled = random_sampling(
+        model,
+        n_samples,
+        method="random_objective",
+        n_objectives=2,
+        replace_max_bound=True,
+        suppress_errors=True,
+        good_reactions=good_reactions,
+        min_flux=include_formate if min_flux is None else min_flux,
+        seed=seed,
+        n_proc=n_proc,
+    )
+    samples = sampled.samples
+    result = SamplingResult(
+        condition=condition.name,
+        include_formate=include_formate,
+        reactions=list(samples.columns),
+        mean=samples.mean(axis=0).to_numpy(),
+        sd=samples.std(axis=0, ddof=1).to_numpy(),
+        n_samples=len(samples),
+        max_maintenance=highest,
+    )
+    returned = sampled.good_reactions
+    return result, list(returned) if returned is not None else (good_reactions or [])
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
new file mode 100644
index 0000000..b05ae54
--- /dev/null
+++ b/tests/test_sampling.py
@@ -0,0 +1,156 @@
+"""Sampling the conventional model under the measured rates.
+
+The ATP accounting is arithmetic on a flux vector, so it is checked
+against the published summary: fed the committed sampling means it must
+reproduce the committed `selectedFluxes.txt`.
+"""
+import numpy as np
+import pandas as pd
+import pytest
+
+from overflow.atp import atp_budget, cofactor_turnover
+from overflow.config import CONDITION_ORDER, ROOT, load_conditions
+from overflow.sampling import (
+    FORMATE_DEHYDROGENASE,
+    FORMATE_EXCHANGE,
+    MEASURED_EXCHANGES,
+    measured_rate,
+    sampling_bounds,
+)
+
+LEGACY = ROOT / "legacy_matlab" / "results" / "randomSampling"
+
+#: yeast-GEM 8.3.4 split glucose phosphorylation over two reactions.
+LEGACY_GLYCOLYSIS_CONSUMING = ("r_4235", "r_0886", "r_0534")
+
+
+@pytest.fixture(scope="module")
+def legacy_means():
+    table = pd.read_csv(LEGACY / "allFluxes.txt", sep="\t")
+    return {
+        condition: dict(zip(table["ID"], table[f"{condition}_AVERAGE"]))
+        for condition in CONDITION_ORDER
+    }
+
+
+@pytest.fixture(scope="module")
+def legacy_selected():
+    return pd.read_csv(LEGACY / "selectedFluxes.txt", sep="\t").set_index("Row")
+
+
+# --- the bounds -------------------------------------------------------
+
+def test_uptakes_are_negative_and_secretions_positive(conditions):
+    bounds = sampling_bounds(conditions["CN75"])
+    low, high = bounds["r_1714"]
+    assert low < high < 0, "glucose is taken up"
+    low, high = bounds["r_1761"]
+    assert 0 < low < high, "ethanol is secreted"
+
+
+def test_every_measured_rate_is_held_within_five_percent(conditions):
+    condition = conditions["CN75"]
+    for reaction_id, (field_name, sign) in MEASURED_EXCHANGES.items():
+        value = measured_rate(condition, field_name)
+        if value == 0:
+            continue
+        low, high = bounds = sampling_bounds(condition)[reaction_id]
+        assert abs(low) == pytest.approx(0.95 * value if sign > 0 else 1.05 * value)
+        assert abs(high) == pytest.approx(1.05 * value if sign > 0 else 0.95 * value)
+
+
+def test_an_undetected_byproduct_gets_a_small_allowance(conditions):
+    """Blocking it would prevent the sampler from saying whether the
+    model wants to make it, which is the question the run asks."""
+    low, high = sampling_bounds(conditions["CN4"])["r_1761"]
+    assert low == 0.0
+    assert high == pytest.approx(0.01)
+
+
+def test_formate_can_be_left_out_of_the_measurement_set(conditions):
+    with_formate = sampling_bounds(conditions["CN22"], include_formate=True)
+    without = sampling_bounds(conditions["CN22"], include_formate=False)
+    assert FORMATE_EXCHANGE in with_formate
+    assert FORMATE_EXCHANGE not in without
+    assert set(with_formate) - {FORMATE_EXCHANGE} == set(without)
+
+
+def test_growth_is_held_at_the_dilution_rate(conditions):
+    low, high = sampling_bounds(conditions["hGR"])["r_2111"]
+    assert low == pytest.approx(0.95 * 0.29)
+    assert high == pytest.approx(1.05 * 0.29)
+
+
+# --- the ATP accounting, against the published summary ---------------
+
+@pytest.mark.parametrize("condition", CONDITION_ORDER)
+def test_atp_budget_reproduces_the_published_summary(
+    legacy_means, legacy_selected, condition
+):
+    means = legacy_means[condition]
+    published = legacy_selected[condition]
+    gam = float(published["GAEC_rATP"]) / means["r_4041"]
+
+    ours = atp_budget(
+        means,
+        gam=gam,
+        growth_rate=load_conditions()[condition].d_rate,
+        glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING,
+    )
+
+    # The published table is written to four significant figures, so a
+    # quantity derived from several of them carries their rounding.
+    for row in (
+        "rGlu", "ETC_rATP", "ETC_YATP_glu", "ETC_YATP_mu",
+        "glycolysis_rATP", "glycolysis_YATP_glu", "glycolysis_YATP_mu",
+        "GAEC_rATP", "GAEC_YATP_glu", "NGAM_rATP", "NGAM_YATP_glu",
+        "Metabolism_rATP", "GAEC+NGAM+Metabolism_rATP",
+        "rPDH", "rIDH", "rMDHc", "rMDHm", "rNDE",
+    ):
+        assert float(ours[row]) == pytest.approx(float(published[row]), rel=2e-3), row
+
+
+def test_the_budget_balances(legacy_means, legacy_selected):
+    """What respiration and glycolysis make is what growth, maintenance
+    and the rest of metabolism spend."""
+    means = legacy_means["CN4"]
+    gam = float(legacy_selected.loc["GAEC_rATP", "CN4"]) / means["r_4041"]
+    budget = atp_budget(means, gam=gam, growth_rate=0.1,
+                        glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING)
+    made = budget["ETC_rATP"] + budget["glycolysis_rATP"]
+    spent = budget["GAEC+NGAM+Metabolism_rATP"]
+    assert made == pytest.approx(spent)
+
+
+def test_yields_are_rates_divided_by_uptake_and_growth(legacy_means, legacy_selected):
+    means = legacy_means["CN38"]
+    gam = float(legacy_selected.loc["GAEC_rATP", "CN38"]) / means["r_4041"]
+    budget = atp_budget(means, gam=gam, growth_rate=0.1,
+                        glycolysis_consuming=LEGACY_GLYCOLYSIS_CONSUMING)
+    assert budget["ETC_YATP_glu"] == pytest.approx(budget["ETC_rATP"] / budget["rGlu"])
+    assert budget["ETC_YATP_mu"] == pytest.approx(budget["ETC_rATP"] / 0.1)
+
+
+def test_a_condition_with_no_glucose_does_not_divide_by_zero():
+    budget = atp_budget({"r_1166": 0.0, "r_4041": 0.1}, gam=50.0, growth_rate=0.1)
+    assert np.isnan(budget["ETC_YATP_glu"])
+    assert np.isfinite(budget["ETC_YATP_mu"])
+
+
+# --- cofactor turnover ------------------------------------------------
+
+@pytest.mark.slow
+def test_cofactor_turnover_counts_production_only(conv_model, legacy_means):
+    turnover = cofactor_turnover(conv_model, legacy_means["CN4"], "NAD", "rNAD")
+    assert (turnover >= 0).all()
+    assert turnover["rNAD[tot]"] == pytest.approx(
+        turnover.drop("rNAD[tot]").sum()
+    )
+
+
+@pytest.mark.slow
+def test_cofactor_turnover_is_reported_per_compartment(conv_model, legacy_means):
+    turnover = cofactor_turnover(conv_model, legacy_means["hGR"], "NAD", "rNAD")
+    compartments = {k for k in turnover.index if k != "rNAD[tot]"}
+    assert len(compartments) >= 2, "NAD turns over in more than one compartment"
+    assert all(k.startswith("rNAD[") for k in compartments)

From 18f4df70f86de912f6a48fa5ab103bc97ce3f500 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 01:54:52 +0200
Subject: [PATCH 2/7] feat: cache the loop-free screening and make bound
 replacement optional

Screening every reaction for loop involvement is the dominant cost of a
sampling run and depends only on the network, so it is kept across runs.

Opening the arbitrary 1000 bounds to infinity, as the published analysis
did, is now off by default: the screening that picks the random
objectives runs on the finite bounds, so once they are opened an
objective can turn out unbounded, and this sampler raises on that where
RAVEN's returned no solution and moved on.
---
 src/overflow/run_sampling.py | 57 +++++++++++++++++++++++++++++++++---
 src/overflow/sampling.py     | 10 ++++++-
 tests/test_sampling.py       | 23 +++++++++++++++
 3 files changed, 85 insertions(+), 5 deletions(-)

diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index a90482d..fb47339 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -9,6 +9,7 @@
 from __future__ import annotations
 
 import argparse
+import json
 import time
 from pathlib import Path
 from typing import Optional
@@ -34,18 +35,48 @@
 DETECTION = 0.001
 
 
+def load_good_reactions(cache: Optional[Path], name: str) -> Optional[list[str]]:
+    """The cached loop-free reaction set, if one was kept."""
+    if cache is None:
+        return None
+    path = cache / f"goodReactions_{name}.json"
+    if not path.is_file():
+        return None
+    return json.loads(path.read_text())
+
+
+def save_good_reactions(cache: Optional[Path], name: str, reactions: list[str]) -> None:
+    if cache is None or not reactions:
+        return
+    cache.mkdir(parents=True, exist_ok=True)
+    (cache / f"goodReactions_{name}.json").write_text(json.dumps(sorted(reactions)))
+
+
 def run(
     conditions: Optional[list[str]] = None,
     n_samples: int = 5000,
     seed: Optional[int] = 0,
     n_proc: Optional[int] = None,
+    cache: Optional[Path] = None,
+    replace_max_bound: bool = False,
     verbose: bool = True,
 ) -> dict[str, dict]:
-    """Sample each condition twice and summarise."""
+    """Sample each condition twice and summarise.
+
+    Which reactions are free of loops depends on the network and not on
+    the condition, so that screening is done once and reused; ``cache``
+    keeps it across runs, where it is otherwise the dominant cost.
+    """
     all_conditions = load_conditions()
     results: dict[str, dict] = {}
-    good_free: Optional[list[str]] = None
-    good_full: Optional[list[str]] = None
+    good_free = load_good_reactions(cache, "free")
+    good_full = load_good_reactions(cache, "full")
+    if verbose and (good_free or good_full):
+        print(
+            f"reusing cached loop-free sets: "
+            f"{len(good_free or [])} / {len(good_full or [])} reactions",
+            flush=True,
+        )
 
     for name in conditions or CONDITION_ORDER:
         condition = all_conditions[name]
@@ -59,16 +90,20 @@ def run(
         free, good_free = sample_condition(
             model.copy(), condition, n_samples=n_samples, include_formate=False,
             seed=seed, good_reactions=good_free, n_proc=n_proc,
+            replace_max_bound=replace_max_bound,
         )
         full, good_full = sample_condition(
             model, condition, n_samples=n_samples, include_formate=True,
             seed=seed, good_reactions=good_full, n_proc=n_proc,
+            replace_max_bound=replace_max_bound,
         )
 
         summary = selected_fluxes(
             model, full.means, gam=gam, growth_rate=condition.d_rate,
             polymerization=polymerization,
         )
+        save_good_reactions(cache, "free", good_free)
+        save_good_reactions(cache, "full", good_full)
         results[name] = {
             "free": free, "full": full, "summary": summary,
             "gam": gam, "model": model,
@@ -109,6 +144,16 @@ def main(argv: Optional[list[str]] = None) -> int:
     parser.add_argument("--samples", type=int, default=5000)
     parser.add_argument("--seed", type=int, default=0)
     parser.add_argument("--procs", type=int, default=None)
+    parser.add_argument(
+        "--cache", type=Path, default=None,
+        help="directory to keep the loop-free reaction screening in",
+    )
+    parser.add_argument(
+        "--replace-max-bound", action="store_true",
+        help="open the arbitrary 1000 bounds to infinity, as the published "
+             "analysis did; the sampler raises if a random objective then "
+             "turns out unbounded",
+    )
     parser.add_argument("--solver")
     args = parser.parse_args(argv)
 
@@ -117,7 +162,11 @@ def main(argv: Optional[list[str]] = None) -> int:
 
         cobra.Configuration().solver = args.solver
 
-    results = run(args.conditions, n_samples=args.samples, seed=args.seed, n_proc=args.procs)
+    results = run(
+        args.conditions, n_samples=args.samples, seed=args.seed,
+        n_proc=args.procs, cache=args.cache,
+        replace_max_bound=args.replace_max_bound,
+    )
 
     out = args.results_dir / "randomSampling"
     out.mkdir(parents=True, exist_ok=True)
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
index 304f21e..408e60e 100644
--- a/src/overflow/sampling.py
+++ b/src/overflow/sampling.py
@@ -185,12 +185,20 @@ def sample_condition(
     good_reactions: Optional[list[str]] = None,
     n_proc: Optional[int] = None,
     min_flux: Optional[bool] = None,
+    replace_max_bound: bool = False,
 ) -> tuple[SamplingResult, list[str]]:
     """Sample one condition, returning the summary and the good reactions.
 
     ``good_reactions`` is the loop-free set the sampler may use as random
     objectives; it depends only on the network, so it is computed once
     and handed back for reuse.
+
+    ``replace_max_bound`` opens the arbitrary 1000 bounds to infinity, as
+    the published analysis did, so that a reaction in a loop cannot sit
+    at 1000 and be mistaken for flux. It is off here: the screening that
+    picks the random objectives runs on the finite bounds, so once they
+    are opened an objective can turn out unbounded, and this sampler
+    raises on that where RAVEN's returned no solution and moved on.
     """
     apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate))
     highest = constrain_maintenance(model)
@@ -200,7 +208,7 @@ def sample_condition(
         n_samples,
         method="random_objective",
         n_objectives=2,
-        replace_max_bound=True,
+        replace_max_bound=replace_max_bound,
         suppress_errors=True,
         good_reactions=good_reactions,
         min_flux=include_formate if min_flux is None else min_flux,
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index b05ae54..f7a3de4 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -154,3 +154,26 @@ def test_cofactor_turnover_is_reported_per_compartment(conv_model, legacy_means)
     compartments = {k for k in turnover.index if k != "rNAD[tot]"}
     assert len(compartments) >= 2, "NAD turns over in more than one compartment"
     assert all(k.startswith("rNAD[") for k in compartments)
+
+
+# --- the loop-free screening cache -----------------------------------
+
+def test_the_loop_free_screening_round_trips(tmp_path):
+    """Screening every reaction for loop involvement is the dominant
+    cost of a sampling run and depends only on the network, so it is
+    kept rather than repeated."""
+    from overflow.run_sampling import load_good_reactions, save_good_reactions
+
+    assert load_good_reactions(tmp_path, "free") is None
+    save_good_reactions(tmp_path, "free", ["r_2", "r_1"])
+    assert load_good_reactions(tmp_path, "free") == ["r_1", "r_2"]
+    assert load_good_reactions(tmp_path, "full") is None
+
+
+def test_no_cache_directory_means_no_caching(tmp_path):
+    from overflow.run_sampling import load_good_reactions, save_good_reactions
+
+    save_good_reactions(None, "free", ["r_1"])
+    assert load_good_reactions(None, "free") is None
+    save_good_reactions(tmp_path, "free", [])
+    assert load_good_reactions(tmp_path, "free") is None

From 020dbbc344ba9305f18510de100617f105dcbaf0 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 02:31:23 +0200
Subject: [PATCH 3/7] perf: give the loop-free screening the process count too

The screening is a flux variability analysis, which reads its process
count from cobra's configuration rather than from the sampler argument.
Left unset it ran on one core and took most of the run.
---
 src/overflow/run_sampling.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index fb47339..66ff820 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -157,10 +157,17 @@ def main(argv: Optional[list[str]] = None) -> int:
     parser.add_argument("--solver")
     args = parser.parse_args(argv)
 
-    if args.solver:
+    if args.solver or args.procs:
         import cobra
 
-        cobra.Configuration().solver = args.solver
+        if args.solver:
+            cobra.Configuration().solver = args.solver
+        if args.procs:
+            # The loop-free screening is a flux variability analysis, which
+            # reads its process count from here rather than from the
+            # sampler's argument. Left unset it runs on one core and
+            # dominates the run.
+            cobra.Configuration().processes = args.procs
 
     results = run(
         args.conditions, n_samples=args.samples, seed=args.seed,

From 4c72d46737103ff244f1de7222536c5413e98356 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 11:03:22 +0200
Subject: [PATCH 4/7] feat: make the parsimonious step within a draw optional

The published second pass minimised total flux within each draw to
suppress loops. Here one numerically awkward draw in several thousand
raises out of the parsimonious solve and aborts the whole run, where
RAVEN's solveLP returned no solution and the loop simply retried, so it
is off unless asked for.
---
 src/overflow/run_sampling.py | 12 +++++++++---
 src/overflow/sampling.py     |  5 +++++
 2 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index 66ff820..a3620e1 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -59,6 +59,7 @@ def run(
     n_proc: Optional[int] = None,
     cache: Optional[Path] = None,
     replace_max_bound: bool = False,
+    min_flux: bool = False,
     verbose: bool = True,
 ) -> dict[str, dict]:
     """Sample each condition twice and summarise.
@@ -90,12 +91,12 @@ def run(
         free, good_free = sample_condition(
             model.copy(), condition, n_samples=n_samples, include_formate=False,
             seed=seed, good_reactions=good_free, n_proc=n_proc,
-            replace_max_bound=replace_max_bound,
+            replace_max_bound=replace_max_bound, min_flux=min_flux,
         )
         full, good_full = sample_condition(
             model, condition, n_samples=n_samples, include_formate=True,
             seed=seed, good_reactions=good_full, n_proc=n_proc,
-            replace_max_bound=replace_max_bound,
+            replace_max_bound=replace_max_bound, min_flux=min_flux,
         )
 
         summary = selected_fluxes(
@@ -154,6 +155,11 @@ def main(argv: Optional[list[str]] = None) -> int:
              "analysis did; the sampler raises if a random objective then "
              "turns out unbounded",
     )
+    parser.add_argument(
+        "--min-flux", action="store_true",
+        help="minimise total flux within each draw, as the published second "
+             "pass did; one numerically awkward draw then aborts the run",
+    )
     parser.add_argument("--solver")
     args = parser.parse_args(argv)
 
@@ -172,7 +178,7 @@ def main(argv: Optional[list[str]] = None) -> int:
     results = run(
         args.conditions, n_samples=args.samples, seed=args.seed,
         n_proc=args.procs, cache=args.cache,
-        replace_max_bound=args.replace_max_bound,
+        replace_max_bound=args.replace_max_bound, min_flux=args.min_flux,
     )
 
     out = args.results_dir / "randomSampling"
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
index 408e60e..b431a54 100644
--- a/src/overflow/sampling.py
+++ b/src/overflow/sampling.py
@@ -199,6 +199,11 @@ def sample_condition(
     picks the random objectives runs on the finite bounds, so once they
     are opened an objective can turn out unbounded, and this sampler
     raises on that where RAVEN's returned no solution and moved on.
+
+    ``min_flux`` minimises total flux within each draw. The published
+    second pass did this to suppress loops; here one awkward draw in
+    several thousand raises out of the parsimonious solve and takes the
+    whole run with it, so it is off unless asked for.
     """
     apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate))
     highest = constrain_maintenance(model)

From b8e4ec02c1ca912c752f5e0d4be81d47d86f2a0a Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 13:10:07 +0200
Subject: [PATCH 5/7] fix: hold reactions to their loop-free range before
 sampling

The malate dehydrogenases cycled against each other at several hundred
mmol/gDW/h, against 0.1 and 0.9 in the published result: ordinary
variability analysis lets a reaction in a thermodynamically infeasible
cycle reach the model's arbitrary 1000 bound, and sampling then draws
from that room. Loop-free variability analysis puts the same pair at
-2.1 to 2.9 and -4.6 to 4.5.

Every flux distribution free of closed cycles already lies inside those
bounds, so nothing sensible is excluded. It costs one loop-free
variability analysis per condition, about a minute on six cores.

Also: have multiprocessing fork rather than use a forkserver. Python
3.14 made forkserver the default and its helper process does not survive
here, which killed runs partway through as a broken pipe.
---
 src/overflow/run_sampling.py | 13 +++++++
 src/overflow/sampling.py     | 73 ++++++++++++++++++++++++++++++++++++
 tests/test_sampling.py       | 61 ++++++++++++++++++++++++++++++
 3 files changed, 147 insertions(+)

diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index a3620e1..fa69b9f 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -29,6 +29,7 @@
     polymerization_breakdown,
     prepare_model,
     sample_condition,
+    use_fork_start_method,
 )
 
 #: Mean flux above which a byproduct counts as predicted [mmol/gDW/h].
@@ -60,6 +61,7 @@ def run(
     cache: Optional[Path] = None,
     replace_max_bound: bool = False,
     min_flux: bool = False,
+    loopless: bool = True,
     verbose: bool = True,
 ) -> dict[str, dict]:
     """Sample each condition twice and summarise.
@@ -92,11 +94,13 @@ def run(
             model.copy(), condition, n_samples=n_samples, include_formate=False,
             seed=seed, good_reactions=good_free, n_proc=n_proc,
             replace_max_bound=replace_max_bound, min_flux=min_flux,
+            loopless=loopless,
         )
         full, good_full = sample_condition(
             model, condition, n_samples=n_samples, include_formate=True,
             seed=seed, good_reactions=good_full, n_proc=n_proc,
             replace_max_bound=replace_max_bound, min_flux=min_flux,
+            loopless=loopless,
         )
 
         summary = selected_fluxes(
@@ -113,6 +117,7 @@ def run(
             print(
                 f"[{name}] {full.n_samples} samples, GAM {gam:.2f}, "
                 f"maintenance up to {full.max_maintenance:.2f}; "
+                f"{full.loopless_tightened} reactions held to their loop-free range; "
                 f"growth {full.means.get(BIO_RXN, float('nan')):.5f} "
                 f"[{time.time() - started:.0f}s]",
                 flush=True,
@@ -160,9 +165,16 @@ def main(argv: Optional[list[str]] = None) -> int:
         help="minimise total flux within each draw, as the published second "
              "pass did; one numerically awkward draw then aborts the run",
     )
+    parser.add_argument(
+        "--keep-loops", action="store_true",
+        help="do not tighten reactions to their loop-free range first",
+    )
     parser.add_argument("--solver")
     args = parser.parse_args(argv)
 
+    if args.procs:
+        use_fork_start_method()
+
     if args.solver or args.procs:
         import cobra
 
@@ -179,6 +191,7 @@ def main(argv: Optional[list[str]] = None) -> int:
         args.conditions, n_samples=args.samples, seed=args.seed,
         n_proc=args.procs, cache=args.cache,
         replace_max_bound=args.replace_max_bound, min_flux=args.min_flux,
+        loopless=not args.keep_loops,
     )
 
     out = args.results_dir / "randomSampling"
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
index b431a54..c857758 100644
--- a/src/overflow/sampling.py
+++ b/src/overflow/sampling.py
@@ -49,6 +49,27 @@
 FORMATE_EXCHANGE = "r_1793"
 
 
+def use_fork_start_method() -> None:
+    """Have multiprocessing fork rather than use a forkserver.
+
+    Python 3.14 made forkserver the default on Linux. Its helper process
+    does not survive in every environment -- on this cluster it dies and
+    takes the run with it, as a broken pipe or a reset connection partway
+    through sampling. Fork is what cobrapy's process pool was written
+    against and works here.
+    """
+    import multiprocessing
+
+    if "fork" not in multiprocessing.get_all_start_methods():
+        return
+    if multiprocessing.get_start_method(allow_none=True) == "fork":
+        return
+    try:
+        multiprocessing.set_start_method("fork", force=True)
+    except RuntimeError:  # already started a pool; nothing to do
+        pass
+
+
 def measured_rate(condition: Condition, field_name: str) -> float:
     """One measured rate, as a positive magnitude."""
     if field_name == "D":
@@ -150,6 +171,48 @@ def polymerization_breakdown(model: "cobra.Model", gam_base: float = GAM_NO_POLY
     return costs
 
 
+def loopless_bounds(model: "cobra.Model", processes: Optional[int] = None) -> pd.DataFrame:
+    """Flux range of every reaction that does not need a closed cycle.
+
+    Ordinary variability analysis lets a reaction in a thermodynamically
+    infeasible cycle run to whatever arbitrary bound the model carries;
+    the loop-free range is what the reaction can do while the rest of
+    the network stays consistent.
+    """
+    from cobra.flux_analysis import flux_variability_analysis
+
+    if processes is not None:
+        import cobra
+
+        cobra.Configuration().processes = processes
+    return flux_variability_analysis(
+        model, fraction_of_optimum=0.0, loopless="cycleFreeFlux"
+    )
+
+
+def apply_loopless_bounds(
+    model: "cobra.Model", ranges: pd.DataFrame, tolerance: float = 1e-9
+) -> int:
+    """Tighten every reaction to its loop-free range.
+
+    Every flux distribution free of closed cycles already lies inside
+    these bounds, so nothing thermodynamically sensible is excluded --
+    what goes is the room a cycle needs to inflate a flux. Returns how
+    many reactions were tightened.
+    """
+    tightened = 0
+    for reaction_id, row in ranges.iterrows():
+        reaction = model.reactions.get_by_id(reaction_id)
+        low = max(reaction.lower_bound, float(row["minimum"]) - tolerance)
+        high = min(reaction.upper_bound, float(row["maximum"]) + tolerance)
+        if low > high:  # numerical crossing; leave the reaction alone
+            continue
+        if low > reaction.lower_bound or high < reaction.upper_bound:
+            reaction.bounds = (low, high)
+            tightened += 1
+    return tightened
+
+
 @dataclass
 class SamplingResult:
     """Sampled fluxes for one condition and one formate setting."""
@@ -161,6 +224,7 @@ class SamplingResult:
     sd: np.ndarray = field(default_factory=lambda: np.empty(0))
     n_samples: int = 0
     max_maintenance: float = 0.0
+    loopless_tightened: int = 0
 
     def to_frame(self) -> pd.DataFrame:
         with np.errstate(divide="ignore", invalid="ignore"):
@@ -186,6 +250,7 @@ def sample_condition(
     n_proc: Optional[int] = None,
     min_flux: Optional[bool] = None,
     replace_max_bound: bool = False,
+    loopless: bool = True,
 ) -> tuple[SamplingResult, list[str]]:
     """Sample one condition, returning the summary and the good reactions.
 
@@ -208,6 +273,13 @@ def sample_condition(
     apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate))
     highest = constrain_maintenance(model)
 
+    # Without this the malate dehydrogenases cycle against each other at
+    # several hundred mmol/gDW/h, which is the arbitrary 1000 bound
+    # showing through rather than anything the cell does.
+    tightened = 0
+    if loopless:
+        tightened = apply_loopless_bounds(model, loopless_bounds(model, n_proc))
+
     sampled = random_sampling(
         model,
         n_samples,
@@ -229,6 +301,7 @@ def sample_condition(
         sd=samples.std(axis=0, ddof=1).to_numpy(),
         n_samples=len(samples),
         max_maintenance=highest,
+        loopless_tightened=tightened,
     )
     returned = sampled.good_reactions
     return result, list(returned) if returned is not None else (good_reactions or [])
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index f7a3de4..9fa6e12 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -177,3 +177,64 @@ def test_no_cache_directory_means_no_caching(tmp_path):
     assert load_good_reactions(None, "free") is None
     save_good_reactions(tmp_path, "free", [])
     assert load_good_reactions(tmp_path, "free") is None
+
+
+# --- holding reactions to their loop-free range ----------------------
+
+def test_loopless_bounds_tighten_only_what_a_cycle_inflates():
+    """A reaction whose ordinary range reaches the model's arbitrary
+    bound, but whose loop-free range is a few units, is held to the
+    latter; a reaction already inside its loop-free range is left be."""
+    import cobra
+    import pandas as pd
+
+    from overflow.sampling import apply_loopless_bounds
+
+    model = cobra.Model("m")
+    a = cobra.Reaction("a", lower_bound=-1000.0, upper_bound=1000.0)
+    b = cobra.Reaction("b", lower_bound=0.0, upper_bound=3.0)
+    model.add_reactions([a, b])
+
+    ranges = pd.DataFrame(
+        {"minimum": [-2.079, 0.0], "maximum": [2.908, 5.0]}, index=["a", "b"]
+    )
+    tightened = apply_loopless_bounds(model, ranges)
+
+    assert tightened == 1
+    assert a.lower_bound == pytest.approx(-2.079, abs=1e-6)
+    assert a.upper_bound == pytest.approx(2.908, abs=1e-6)
+    assert b.bounds == (0.0, 3.0), "already inside its loop-free range"
+
+
+def test_loopless_bounds_never_widen_a_reaction():
+    """The loop-free range of a reaction can exceed a bound the
+    condition imposes; the condition wins."""
+    import cobra
+    import pandas as pd
+
+    from overflow.sampling import apply_loopless_bounds
+
+    model = cobra.Model("m")
+    r = cobra.Reaction("r", lower_bound=-1.0, upper_bound=1.0)
+    model.add_reactions([r])
+    apply_loopless_bounds(
+        model, pd.DataFrame({"minimum": [-50.0], "maximum": [50.0]}, index=["r"])
+    )
+    assert r.bounds == (-1.0, 1.0)
+
+
+def test_a_crossed_range_leaves_the_reaction_alone():
+    """Numerical noise can put a loop-free minimum above its maximum;
+    that must not produce an infeasible reaction."""
+    import cobra
+    import pandas as pd
+
+    from overflow.sampling import apply_loopless_bounds
+
+    model = cobra.Model("m")
+    r = cobra.Reaction("r", lower_bound=0.0, upper_bound=1.0)
+    model.add_reactions([r])
+    apply_loopless_bounds(
+        model, pd.DataFrame({"minimum": [5.0], "maximum": [2.0]}, index=["r"])
+    )
+    assert r.bounds == (0.0, 1.0)

From e7933d7607bee82b46b99a47541dd0680c879520 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Sun, 20 Sep 2026 13:51:45 +0200
Subject: [PATCH 6/7] fix: take the objective set from the same loop-free
 analysis

The reactions worth using as random objectives are the ones that can
carry flux without a closed cycle, which is what the loop-free
variability analysis already reports. Deriving both from it costs
nothing extra and, more to the point, makes the objective set belong to
the condition being sampled: reusing one condition's set for another
leaves reactions in it that are pinned here, and each time one is picked
the sampler retries a hundred times before giving up on that draw. hGR
spent thirty minutes on what takes the other conditions two.
---
 src/overflow/run_sampling.py | 11 +++++++----
 src/overflow/sampling.py     | 23 ++++++++++++++++++++++-
 2 files changed, 29 insertions(+), 5 deletions(-)

diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index fa69b9f..6b194b0 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -72,8 +72,10 @@ def run(
     """
     all_conditions = load_conditions()
     results: dict[str, dict] = {}
-    good_free = load_good_reactions(cache, "free")
-    good_full = load_good_reactions(cache, "full")
+    # With loop-free bounding the objective set is derived per condition
+    # from that condition's own analysis, so a cached one would be ignored.
+    good_free = None if loopless else load_good_reactions(cache, "free")
+    good_full = None if loopless else load_good_reactions(cache, "full")
     if verbose and (good_free or good_full):
         print(
             f"reusing cached loop-free sets: "
@@ -107,8 +109,9 @@ def run(
             model, full.means, gam=gam, growth_rate=condition.d_rate,
             polymerization=polymerization,
         )
-        save_good_reactions(cache, "free", good_free)
-        save_good_reactions(cache, "full", good_full)
+        if not loopless:
+            save_good_reactions(cache, "free", good_free)
+            save_good_reactions(cache, "full", good_full)
         results[name] = {
             "free": free, "full": full, "summary": summary,
             "gam": gam, "model": model,
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
index c857758..d17404a 100644
--- a/src/overflow/sampling.py
+++ b/src/overflow/sampling.py
@@ -190,6 +190,21 @@ def loopless_bounds(model: "cobra.Model", processes: Optional[int] = None) -> pd
     )
 
 
+def good_reactions_from(ranges: pd.DataFrame, tolerance: float = 1e-9) -> list[str]:
+    """Reactions that can carry flux without a closed cycle.
+
+    These are the ones worth using as random objectives; a reaction that
+    only moves inside a cycle gives a degenerate solve every time it is
+    picked, and the sampler retries up to a hundred times before giving
+    up on that draw.
+    """
+    return [
+        reaction_id
+        for reaction_id, row in ranges.iterrows()
+        if max(abs(float(row["minimum"])), abs(float(row["maximum"]))) > tolerance
+    ]
+
+
 def apply_loopless_bounds(
     model: "cobra.Model", ranges: pd.DataFrame, tolerance: float = 1e-9
 ) -> int:
@@ -278,7 +293,13 @@ def sample_condition(
     # showing through rather than anything the cell does.
     tightened = 0
     if loopless:
-        tightened = apply_loopless_bounds(model, loopless_bounds(model, n_proc))
+        # One loop-free variability analysis gives both the bounds and
+        # the objective set, and the objective set has to come from this
+        # condition: a reaction that moves freely in one condition can be
+        # pinned in another, and picking it then wastes a whole draw.
+        ranges = loopless_bounds(model, n_proc)
+        tightened = apply_loopless_bounds(model, ranges)
+        good_reactions = good_reactions_from(ranges)
 
     sampled = random_sampling(
         model,

From 4c90c018cba7b7a46a7b6927b5176aaaefc54846 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Fri, 25 Sep 2026 12:55:18 +0200
Subject: [PATCH 7/7] fix: minimise total flux within each draw, and bound the
 loop search

The sampled means drift into high-flux routes without the second step the
published analysis had: the pentose phosphate pathway and a reversed NADP
isocitrate dehydrogenase carry several times the published flux. The step
now tolerates a draw whose optimum cannot be pinned exactly.

cobra picks the cyclic reactions with unseeded random LP weights, and one
seed in a while runs for hours. The search now has a time limit and retries
with a new seed.
---
 README.md                    |  8 ++--
 src/overflow/run_sampling.py |  8 ++--
 src/overflow/sampling.py     | 77 +++++++++++++++++++++++++++++++-----
 tests/test_sampling.py       | 60 ++++++++++++++++++++++++++++
 4 files changed, 136 insertions(+), 17 deletions(-)

diff --git a/README.md b/README.md
index 16d0894..0fa5161 100644
--- a/README.md
+++ b/README.md
@@ -121,9 +121,11 @@ formate in the measured set, and writes the sampled means and standard
 deviations, the byproducts the model secretes when it is not told to, and the
 ATP and redox budget.
 
-Sampling draws vertices by maximising small random objectives, as the published
-analysis did; means over such draws are vertex means rather than an average over
-the interior of the flux space.
+Sampling draws vertices by maximising small random objectives and then minimising
+total flux at each, as the published analysis did; means over such draws are vertex
+means rather than an average over the interior of the flux space. `--no-min-flux`
+skips the second step, which lets the sampler wander into high-flux routes: the
+pentose phosphate pathway then carries several times the published flux.
 
 ## Tests
 
diff --git a/src/overflow/run_sampling.py b/src/overflow/run_sampling.py
index 6b194b0..ad5fa62 100644
--- a/src/overflow/run_sampling.py
+++ b/src/overflow/run_sampling.py
@@ -60,7 +60,7 @@ def run(
     n_proc: Optional[int] = None,
     cache: Optional[Path] = None,
     replace_max_bound: bool = False,
-    min_flux: bool = False,
+    min_flux: bool = True,
     loopless: bool = True,
     verbose: bool = True,
 ) -> dict[str, dict]:
@@ -164,9 +164,9 @@ def main(argv: Optional[list[str]] = None) -> int:
              "turns out unbounded",
     )
     parser.add_argument(
-        "--min-flux", action="store_true",
-        help="minimise total flux within each draw, as the published second "
-             "pass did; one numerically awkward draw then aborts the run",
+        "--no-min-flux", dest="min_flux", action="store_false",
+        help="skip the minimisation of total flux within each draw that the "
+             "published second pass did",
     )
     parser.add_argument(
         "--keep-loops", action="store_true",
diff --git a/src/overflow/sampling.py b/src/overflow/sampling.py
index d17404a..8ff280a 100644
--- a/src/overflow/sampling.py
+++ b/src/overflow/sampling.py
@@ -11,11 +11,16 @@
 """
 from __future__ import annotations
 
+import signal
+from contextlib import contextmanager
 from dataclasses import dataclass, field
 from typing import TYPE_CHECKING, Mapping, Optional
 
 import numpy as np
 import pandas as pd
+import raven_toolbox.analysis.sampling as raven_sampling
+from cobra.exceptions import Infeasible
+from cobra.flux_analysis import pfba
 from raven_toolbox.analysis.sampling import random_sampling
 
 from overflow.biomass import (
@@ -49,6 +54,22 @@
 FORMATE_EXCHANGE = "r_1793"
 
 
+@contextmanager
+def time_limit(seconds: float):
+    """Raise ``TimeoutError`` in the main thread once ``seconds`` have passed."""
+
+    def expire(signum, frame):
+        raise TimeoutError
+
+    previous = signal.signal(signal.SIGALRM, expire)
+    signal.setitimer(signal.ITIMER_REAL, seconds)
+    try:
+        yield
+    finally:
+        signal.setitimer(signal.ITIMER_REAL, 0)
+        signal.signal(signal.SIGALRM, previous)
+
+
 def use_fork_start_method() -> None:
     """Have multiprocessing fork rather than use a forkserver.
 
@@ -70,6 +91,19 @@ def use_fork_start_method() -> None:
         pass
 
 
+def tolerant_pfba(model: "cobra.Model", fraction_of_optimum: float = 1 - 1e-6):
+    """Minimise total flux at the optimum, keeping the plain solution if that is infeasible.
+
+    Pinning the objective exactly at its optimum can be numerically
+    infeasible for a single draw; a small slack avoids most of those and
+    the plain solution covers the rest.
+    """
+    try:
+        return pfba(model, fraction_of_optimum=fraction_of_optimum)
+    except Infeasible:
+        return model.optimize()
+
+
 def measured_rate(condition: Condition, field_name: str) -> float:
     """One measured rate, as a positive magnitude."""
     if field_name == "D":
@@ -171,7 +205,13 @@ def polymerization_breakdown(model: "cobra.Model", gam_base: float = GAM_NO_POLY
     return costs
 
 
-def loopless_bounds(model: "cobra.Model", processes: Optional[int] = None) -> pd.DataFrame:
+def loopless_bounds(
+    model: "cobra.Model",
+    processes: Optional[int] = None,
+    seed: int = 0,
+    attempts: int = 4,
+    seconds: float = 900.0,
+) -> pd.DataFrame:
     """Flux range of every reaction that does not need a closed cycle.
 
     Ordinary variability analysis lets a reaction in a thermodynamically
@@ -185,8 +225,22 @@ def loopless_bounds(model: "cobra.Model", processes: Optional[int] = None) -> pd
         import cobra
 
         cobra.Configuration().processes = processes
-    return flux_variability_analysis(
-        model, fraction_of_optimum=0.0, loopless="cycleFreeFlux"
+    for attempt in range(attempts):
+        # cobra picks the cyclic reactions with random LP weights, and how
+        # long that takes depends on them: most seeds finish in a minute or
+        # two, an occasional one runs for hours. A fresh seed after a time
+        # limit gets past it.
+        np.random.seed(seed + attempt)
+        try:
+            with time_limit(seconds):
+                return flux_variability_analysis(
+                    model, fraction_of_optimum=0.0, loopless="cycleFreeFlux"
+                )
+        except TimeoutError:
+            continue
+    raise TimeoutError(
+        f"the loop-free variability analysis did not finish in {seconds:.0f} s "
+        f"in any of {attempts} attempts"
     )
 
 
@@ -263,7 +317,7 @@ def sample_condition(
     seed: Optional[int] = None,
     good_reactions: Optional[list[str]] = None,
     n_proc: Optional[int] = None,
-    min_flux: Optional[bool] = None,
+    min_flux: bool = True,
     replace_max_bound: bool = False,
     loopless: bool = True,
 ) -> tuple[SamplingResult, list[str]]:
@@ -280,12 +334,15 @@ def sample_condition(
     are opened an objective can turn out unbounded, and this sampler
     raises on that where RAVEN's returned no solution and moved on.
 
-    ``min_flux`` minimises total flux within each draw. The published
-    second pass did this to suppress loops; here one awkward draw in
-    several thousand raises out of the parsimonious solve and takes the
-    whole run with it, so it is off unless asked for.
+    ``min_flux`` minimises total flux within each draw, as the published
+    second pass did. Without it the sampled means drift into
+    high-flux routes: the pentose phosphate pathway and a reversed
+    NADP isocitrate dehydrogenase carry several times the published
+    flux. The step is tolerant of a draw whose optimum cannot be pinned
+    exactly.
     """
     apply_bounds(model, sampling_bounds(condition, tolerance, include_formate=include_formate))
+    raven_sampling.pfba = tolerant_pfba
     highest = constrain_maintenance(model)
 
     # Without this the malate dehydrogenases cycle against each other at
@@ -297,7 +354,7 @@ def sample_condition(
         # the objective set, and the objective set has to come from this
         # condition: a reaction that moves freely in one condition can be
         # pinned in another, and picking it then wastes a whole draw.
-        ranges = loopless_bounds(model, n_proc)
+        ranges = loopless_bounds(model, n_proc, seed=seed or 0)
         tightened = apply_loopless_bounds(model, ranges)
         good_reactions = good_reactions_from(ranges)
 
@@ -309,7 +366,7 @@ def sample_condition(
         replace_max_bound=replace_max_bound,
         suppress_errors=True,
         good_reactions=good_reactions,
-        min_flux=include_formate if min_flux is None else min_flux,
+        min_flux=min_flux,
         seed=seed,
         n_proc=n_proc,
     )
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index 9fa6e12..34a5d8e 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -17,6 +17,7 @@
     measured_rate,
     sampling_bounds,
 )
+from routes import two_route_model
 
 LEGACY = ROOT / "legacy_matlab" / "results" / "randomSampling"
 
@@ -238,3 +239,62 @@ def test_a_crossed_range_leaves_the_reaction_alone():
         model, pd.DataFrame({"minimum": [5.0], "maximum": [2.0]}, index=["r"])
     )
     assert r.bounds == (0.0, 1.0)
+
+
+def test_a_variability_analysis_that_overruns_is_retried_with_a_new_seed(monkeypatch):
+    import time
+
+    import cobra.flux_analysis as flux_analysis
+    import numpy as np
+
+    from overflow.sampling import loopless_bounds
+
+    seeds = []
+
+    def analysis(model, **kwargs):
+        seeds.append(int(np.random.get_state()[1][0]))
+        if len(seeds) == 1:
+            time.sleep(5)
+        return "ranges"
+
+    monkeypatch.setattr(flux_analysis, "flux_variability_analysis", analysis)
+    assert loopless_bounds(object(), seed=7, attempts=3, seconds=0.2) == "ranges"
+    assert len(seeds) == 2 and seeds[0] != seeds[1]
+
+
+def test_a_variability_analysis_that_never_finishes_is_reported(monkeypatch):
+    import time
+
+    import cobra.flux_analysis as flux_analysis
+
+    from overflow.sampling import loopless_bounds
+
+    monkeypatch.setattr(
+        flux_analysis, "flux_variability_analysis", lambda model, **kw: time.sleep(5)
+    )
+    with pytest.raises(TimeoutError, match="2 attempts"):
+        loopless_bounds(object(), attempts=2, seconds=0.1)
+
+
+def test_a_draw_that_cannot_be_pinned_keeps_the_plain_solution(monkeypatch):
+    from cobra.exceptions import Infeasible
+
+    import overflow.sampling as sampling
+
+    model = two_route_model()
+
+    def refuse(*args, **kwargs):
+        raise Infeasible("numerically infeasible")
+
+    monkeypatch.setattr(sampling, "pfba", refuse)
+    solution = sampling.tolerant_pfba(model)
+    assert solution.status == "optimal"
+    assert solution.fluxes["demand"] == pytest.approx(1.0)
+
+
+def test_the_tolerant_step_minimises_total_flux_when_it_can():
+    import overflow.sampling as sampling
+
+    solution = sampling.tolerant_pfba(two_route_model())
+    assert solution.fluxes["direct"] == pytest.approx(0.0, abs=1e-5)
+    assert solution.fluxes["step1"] == pytest.approx(1.0, abs=1e-5)