From 23b8604b4e2966ad7e0574565cecfd7bf6f87305 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 26 Sep 2026 00:16:46 +0200 Subject: [PATCH] feat: close named reactions before building, and close the proton leak --block sets the bounds of the named reactions to zero. The proton leak (r_2129) has no enzyme and no bound, so with the measured oxygen fitted it absorbed most of the respiration at high C/N (P/O 0.19 at CN75); the built models can carry no maintenance ATP because the rate fit raises each cap only as far as the rates need. With the leak closed the respiration is coupled (P/O 1.1-1.2) and the maintenance appears as non-growth ATP use. The README recipe closes it. --- README.md | 12 ++++++++++-- src/overflow/build.py | 14 ++++++++++++-- src/overflow/constraints.py | 13 +++++++++++++ tests/test_build.py | 17 +++++++++++++++++ tests/test_constraints.py | 14 ++++++++++++++ tests/test_objective.py | 13 +++++++++++++ 6 files changed, 79 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 0634b57..86c5c96 100644 --- a/README.md +++ b/README.md @@ -85,6 +85,13 @@ Options change what the models are asked to account for: pool (`protein`, the default), or the smallest total flux at the dilution rate over every reaction (`flux`) or over the metabolic reactions only (`metabolic-flux`). The flux objectives need `--fit-rates`. +- `--block RXN ...` closes the named reactions before anything else is applied. The + analysis blocks `r_2129`, the proton leak: it has no enzyme and no bound, so it is a + free sink for oxygen. Fitting the measured oxygen with it open sent most of the + respiration through it at high C/N (P/O 0.19 at CN75, against 1.1-1.2 with it closed), + because the fit raises each cap only as far as the rates need and the built models + can carry no maintenance ATP. Closed, the respiration is coupled and the maintenance + shows up as non-growth ATP use. - `--scale-protein` rescales biomass to the measured protein content. Off by default: the model has no carbon to spare at the measured glucose uptake, and CN4 then falls short of its dilution rate. @@ -146,10 +153,11 @@ measurements and raises the measured enzyme abundances by the least that makes that feasible. Without it those rates are free, and on the earlier tutorial model the enzyme-constrained solution disposed of surplus carbon through unmeasured exits instead of respiring it, so the flux distribution did not describe the -measured physiology. `--uptake-flex 1.08` gives CN4 the glucose it needs. +measured physiology. `--uptake-flex 1.08` gives CN4 the glucose it needs, and `--block r_2129` closes the +proton leak (see the build options). ```bash -python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 +python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 --block r_2129 python -m overflow.build_ribosome python -m overflow.analyze_usage python -m overflow.run_sampling --procs 12 diff --git a/src/overflow/build.py b/src/overflow/build.py index 33e4bcf..76d5cee 100644 --- a/src/overflow/build.py +++ b/src/overflow/build.py @@ -13,7 +13,7 @@ import time from dataclasses import dataclass, field from pathlib import Path -from typing import Optional +from typing import Optional, Sequence import numpy as np import pandas as pd @@ -39,6 +39,7 @@ InfeasibleCondition, constrain_byproducts, constrain_measured_rates, + block_reactions, constrain_uptake, free_ngam, set_chemostat_constraints, @@ -102,6 +103,7 @@ def build_condition( rate_tolerance: float = 0.05, uptake_flex: float = 1.05, objective: str = "protein", + block: Sequence[str] = (), verbose: bool = True, ) -> BuildResult: """Build one proteome-constrained ecModel. @@ -123,6 +125,8 @@ def build_condition( protein pool (``"protein"``), or the smallest total flux at the dilution rate over every reaction (``"flux"``) or over the metabolic reactions only (``"metabolic-flux"``); both need ``fit_rates``. + + ``block`` lists reactions to close before anything else is applied. """ if objective not in OBJECTIVES: raise ValueError(f"objective must be one of {OBJECTIVES}, not {objective!r}") @@ -142,6 +146,7 @@ def build_condition( model = load_model(build_adapter(condition)) free_ngam(model) + block_reactions(model, block) # Biomass is rescaled to the protein content as measured, not to the # value the proteomics coverage implies. @@ -358,6 +363,11 @@ def main(argv: Optional[list[str]] = None) -> int: "raise abundances by the least that makes that feasible", ) parser.add_argument("--rate-tolerance", type=float, default=0.05) + parser.add_argument( + "--block", nargs="+", default=[], metavar="RXN", + help="reaction IDs to close before the model is built, e.g. r_2129 (the " + "proton leak)", + ) parser.add_argument( "--objective", choices=OBJECTIVES, @@ -403,7 +413,7 @@ def main(argv: Optional[list[str]] = None) -> int: fit_rates=args.fit_rates, rate_tolerance=args.rate_tolerance, uptake_flex=args.uptake_flex, - objective=args.objective, + objective=args.objective, block=args.block, ) write_outputs(result, args.models_dir, args.results_dir) results.append(result) diff --git a/src/overflow/constraints.py b/src/overflow/constraints.py index 4d04921..ce25cf8 100644 --- a/src/overflow/constraints.py +++ b/src/overflow/constraints.py @@ -45,6 +45,19 @@ def constrain_uptake(model: cobra.Model, condition: Condition, flex: float = 1.0 model.reactions.get_by_id(C_SOURCE).bounds = (-flex * condition.glucose, 0.0) +def block_reactions(model: cobra.Model, reaction_ids) -> None: + """Set the bounds of each named reaction to zero. + + An unknown ID is an error rather than a reaction that quietly stays open. + """ + for reaction_id in reaction_ids: + try: + reaction = model.reactions.get_by_id(reaction_id) + except KeyError: + raise KeyError(f"cannot block {reaction_id}: the model has no such reaction") from None + reaction.bounds = (0.0, 0.0) + + def free_ngam(model: cobra.Model, lower: float = 0.0) -> None: """Release the maintenance reaction so it can be fitted. diff --git a/tests/test_build.py b/tests/test_build.py index 6acc746..fc4d482 100644 --- a/tests/test_build.py +++ b/tests/test_build.py @@ -54,3 +54,20 @@ def test_the_summary_reports_the_condition(built): table = summary_table([built]) assert list(table["condition"]) == ["CN4"] assert table.loc[0, "growth"] >= 0.099 + + +@pytest.fixture(scope="module") +def built_without_the_leak(conditions): + return build_condition( + conditions["CN22"], ngam_steps=5, fit_rates=True, rate_tolerance=0.08, + uptake_flex=1.08, block=["r_2129"], verbose=False, + ) + + +def test_a_blocked_leak_carries_no_flux_and_respiration_stays_coupled(built_without_the_leak): + fluxes = built_without_the_leak.fluxes + assert fluxes["r_2129"] == pytest.approx(0.0, abs=1e-9) + atp_synthase = fluxes["r_0226"] - fluxes.get("r_0226_REV", 0.0) + oxygen = -fluxes["r_1992"] + assert atp_synthase / (2 * oxygen) > 0.9 + assert built_without_the_leak.growth == pytest.approx(0.1, rel=0.01) diff --git a/tests/test_constraints.py b/tests/test_constraints.py index 058b462..2d06c6e 100644 --- a/tests/test_constraints.py +++ b/tests/test_constraints.py @@ -2,6 +2,7 @@ from overflow.config import BIO_RXN, BYPRODUCT_RXNS, C_SOURCE, NGAM_RXN, POOL_RXN from overflow.constraints import ( + block_reactions, constrain_byproducts, constrain_uptake, free_ngam, @@ -96,3 +97,16 @@ def test_the_chemostat_objective_survives_being_written_and_read(ec_model, condi after = reloaded.optimize().fluxes[POOL_RXN] assert after == pytest.approx(before, rel=1e-4) # YAML rounds the bounds assert after < reloaded.reactions.get_by_id(POOL_RXN).upper_bound + + +def test_blocked_reactions_carry_no_flux(ec_model): + with ec_model as model: + target = next(r for r in model.reactions if not r.id.startswith(("usage_prot_", "prot_"))) + block_reactions(model, [target.id]) + assert target.bounds == (0.0, 0.0) + + +def test_blocking_an_unknown_reaction_is_an_error(ec_model): + with ec_model as model: + with pytest.raises(KeyError, match="no such reaction"): + block_reactions(model, ["r_9999999"]) diff --git a/tests/test_objective.py b/tests/test_objective.py index 3fd8936..b873a8e 100644 --- a/tests/test_objective.py +++ b/tests/test_objective.py @@ -2,6 +2,7 @@ from overflow.build import OBJECTIVES, build_condition, metabolic_reactions from overflow.config import POOL_RXN +from overflow.constraints import block_reactions from routes import two_route_model from tiny import tiny_ec_model @@ -43,3 +44,15 @@ def test_parsimony_over_a_subset_ignores_the_others_but_returns_every_flux(): assert subset.fluxes["usage_prot_E"] == pytest.approx(10.0) assert set(subset.fluxes.index) == {r.id for r in model.reactions} + +def test_a_blocked_reaction_stays_shut_on_a_small_model(): + model = two_route_model() + block_reactions(model, ["direct"]) + assert model.reactions.direct.bounds == (0.0, 0.0) + assert model.slim_optimize() == pytest.approx(1.0) + assert model.optimize().fluxes["step1"] == pytest.approx(1.0) + + +def test_blocking_an_unknown_reaction_is_refused_on_a_small_model(): + with pytest.raises(KeyError, match="cannot block r_missing"): + block_reactions(two_route_model(), ["r_missing"])