Close the proton leak when fitting the measured rates

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

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)

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.
 

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)

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

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