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"])