From 329c57e8bc5c0c6650030a8f612459e731fac896 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 18:15:34 +0100 Subject: [PATCH 1/4] Rank get_enzyme_bottlenecks by the re-solved change of the objective The order no longer depends on the shadow price, which depends on the basis at a degenerate optimum. Enzymes at their bound are solved again with the bound raised by fold_change; multiple_active marks rows where another capacity constraint is active. --- CHANGELOG.md | 1 + docs/migrating_from_gecko_matlab.md | 4 +- src/geckopy/utilities/bottlenecks.py | 98 ++++++++++++------ tests/test_bottlenecks.py | 149 ++++++++++++++++++++++++++- tests/test_scipy_solver_support.py | 8 +- 5 files changed, 217 insertions(+), 43 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..2a01a3d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 6c3b465..5ea6eba 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -324,7 +324,7 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `startGECKOproject` | `geckopy init` *(CLI)* | **`[Py]`** Project skeleton generator runs from the `geckopy` CLI (`geckopy init my_project`). | | `updateGECKOdoc` | *(not applicable)* | Python docs are docstring-driven; no MATLAB Toolbox metadata to regenerate. | | — | `pfba_enzymes` *(new)* | **`[Py]`** pFBA variant that minimises total `usage_prot_*` flux instead of all fluxes. Raises `NotImplementedError` on gecko-light. | -| — | `get_enzyme_bottlenecks` *(new)* | **`[Py]`** Per-enzyme dual-based bottleneck ranker; returns a DataFrame. | +| — | `get_enzyme_bottlenecks` *(new)* | **`[Py]`** Per-enzyme bottleneck ranker; returns a DataFrame sorted by the re-solved change of the objective when the bound of one enzyme at its bound is raised by `fold_change`. The shadow price is a reference column. | ### 5.7 Model adapter (`src/geckomat/model_adapter`) @@ -456,7 +456,7 @@ for the MATLAB-side items the port surfaced. ### geckopy-only (no MATLAB counterpart) - `pfba_enzymes` — enzyme-pool-minimising pFBA. -- `get_enzyme_bottlenecks` — dual-based bottleneck ranking. +- `get_enzyme_bottlenecks` — bottleneck ranking by the re-solved change of the objective. - `relax_proteomics_greedy` — greedy alternative to `flexibilize_enz_concs`. - `merge_kcats` — generalised n-source kcat merge (replaces the diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index ffdfdef..0c31ead 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -1,43 +1,42 @@ -"""Find the enzymes that are limiting the current objective. - -When you solve an ecModel, some enzymes are saturated and some -have spare capacity. The saturated ones are the ones whose -removal (or whose increased availability) would shift the -objective — they're the *bottlenecks*. The LP's shadow price for -each ``prot_`` mass-balance constraint quantifies that -sensitivity directly. - -``get_enzyme_bottlenecks`` runs the solve, ranks every enzyme by -absolute shadow price, and returns the top N as a DataFrame. The -top row is the enzyme whose extra availability would help the -objective the most. - -Ported from the legacy geckopy package described in Carrasco et al. -(2023, https://doi.org/10.1128/spectrum.01705-23), file -geckopy/flux_analysis.py:275-281 (get_protein_bottlenecks). +"""Find the enzymes that limit the current objective. + +When an ecModel is solved, some enzymes are used up to their upper bound and some have spare capacity. An enzyme that is not used up cannot limit the objective: raising its bound does not change the optimum. An enzyme that is used up may limit the objective, but it does not have to, because another capacity constraint can be active at the same optimum. + +``get_enzyme_bottlenecks`` ranks the enzymes by the re-solved change of the objective when the upper bound of one enzyme is raised by a factor ``fold_change``. Only the enzymes at their bound are re-solved; every other enzyme has a change of 0 without a solve. The LP shadow price of the ``prot_`` constraint is reported for reference and is not used for the order: at a degenerate optimum, which enzyme-constrained models routinely have, several constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis and not on the model. + +Ported from the legacy geckopy package described in Carrasco et al. (2023, https://doi.org/10.1128/spectrum.01705-23), file geckopy/flux_analysis.py:275-281 (get_protein_bottlenecks). """ from __future__ import annotations +import math from typing import TYPE_CHECKING import pandas as pd +from ..ec_model.constants import POOL_EXCHANGE_ID from ..solver_support import constraint_shadow_prices, optimize_model if TYPE_CHECKING: from ..ec_model.ec_model import EcModel +# An enzyme is at its bound when |flux| / upper_bound is at least 1 - _AT_CAP_TOL. +_AT_CAP_TOL = 1e-6 +# Changes of the objective below this value are set to 0. +_GAIN_TOL = 1e-10 + + +def get_enzyme_bottlenecks( + model: "EcModel", *, top: int = 10, fold_change: float = 2.0 +) -> pd.DataFrame: + """Top-N enzymes by the change of the objective when their bound is raised. -def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: - """Top-N enzymes by absolute shadow price. + Solves the model. An enzyme is at its bound when its upper bound is positive and ``|flux| / upper_bound`` is at least ``1 - 1e-6``. For each such enzyme the upper bound of its usage reaction is raised to ``upper_bound * fold_change``, the model is solved again and ``objective_gain`` is the new objective minus the objective of the first solve (0 when the new solve is not optimal or the change is below 1e-10). The model is restored after each solve. Every other enzyme has ``objective_gain`` 0 without a solve, because raising a bound that does not bind cannot change the optimum. The change is signed in the direction of the objective as the solver returns it, so it is positive for a maximised objective that increases. - Solves the model, ranks enzymes by ``|shadow_price|``, returns a - DataFrame indexed by uniprot id with columns - ``[gene, shadow_price, flux, cap_usage, upper_bound]``. + The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. - Ported from the legacy geckopy package (Carrasco et al., 2023, - https://doi.org/10.1128/spectrum.01705-23), - geckopy/flux_analysis.py:275-281. + ``multiple_active`` is True for an enzyme when a capacity constraint other than its own is active at the same optimum, which is where the dual says little about the enzyme. The capacity constraints counted are the enzymes at their bound and the protein pool exchange (``prot_pool_exchange``) when its flux is within 1e-6 relative of a finite, positive upper bound. + + Ported from the legacy geckopy package (Carrasco et al., 2023, https://doi.org/10.1128/spectrum.01705-23), geckopy/flux_analysis.py:275-281. Parameters ---------- @@ -45,14 +44,13 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: EcModel to solve. top Number of top-ranked enzymes to return. Default 10. + fold_change + Factor by which the upper bound of an enzyme at its bound is raised. Default 2.0. Returns ------- pandas.DataFrame - Indexed by uniprot id, columns - ``[gene, shadow_price, flux, cap_usage, upper_bound]``, - sorted by descending ``|shadow_price|`` and truncated to - ``top`` rows. + Indexed by uniprot id, columns ``[gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active]``, truncated to ``top`` rows. ``df.attrs["n_resolved"]`` is the number of enzymes that were solved again and ``df.attrs["active_constraints"]`` the number of active capacity constraints. Raises ------ @@ -64,6 +62,7 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: raise RuntimeError( f"Solver status {sol.status!r}; cannot rank bottlenecks" ) + base = sol.objective_value enzymes = list(model.enzymes) prices = constraint_shadow_prices(model, [enz.prot_metabolite_id for enz in enzymes]) @@ -72,13 +71,46 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: rows.append({ "uniprot": enz.id, "gene": enz.gene, + "objective_gain": 0.0, "shadow_price": float(price), "flux": enz.flux, "cap_usage": enz.cap_usage, "upper_bound": enz.upper_bound, }) + + at_cap = [ + i for i, row in enumerate(rows) + if row["upper_bound"] > 0 and row["cap_usage"] >= 1 - _AT_CAP_TOL + ] + active = len(at_cap) + if model.reactions.has_id(POOL_EXCHANGE_ID): + pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) + pool_ub = float(pool.upper_bound) + if 0 < pool_ub < math.inf and abs(float(pool.flux)) >= pool_ub * (1 - _AT_CAP_TOL): + active += 1 + + for i in at_cap: + usage = enzymes[i].usage_reaction + with model: + usage.upper_bound = rows[i]["upper_bound"] * fold_change + new = optimize_model(model) + if new.status == "optimal": + gain = new.objective_value - base + rows[i]["objective_gain"] = gain if gain >= _GAIN_TOL else 0.0 + if at_cap: + optimize_model(model) + + at_cap_set = set(at_cap) + for i, row in enumerate(rows): + row["multiple_active"] = (active - (1 if i in at_cap_set else 0)) >= 1 + + rows.sort(key=lambda row: -round(row["objective_gain"], 9)) df = pd.DataFrame(rows).set_index("uniprot") - df = df.reindex( - df["shadow_price"].abs().sort_values(ascending=False).index - ) - return df.head(top) + df = df[[ + "gene", "objective_gain", "shadow_price", "flux", "cap_usage", + "upper_bound", "multiple_active", + ]] + df = df.head(top) + df.attrs["n_resolved"] = len(at_cap) + df.attrs["active_constraints"] = active + return df diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index 8283845..d8a2246 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -2,9 +2,13 @@ from pathlib import Path import cobra +import numpy as np import pytest +from scipy import sparse from geckopy import ModelAdapter, make_ec_model +from geckopy.ec_model import EcModel +from geckopy.ec_model.ec_data import EcData from geckopy.ec_model.pipeline import apply_kcat_constraints from geckopy.utilities import get_enzyme_bottlenecks @@ -43,20 +47,155 @@ def test_returns_at_most_top_when_fewer_enzymes(): assert len(df) == len(model.ec.enzymes) -def test_sorted_by_absolute_shadow_price(): +def test_sorted_by_objective_gain(): model = _ectestgem_solvable() - df = get_enzyme_bottlenecks(model, top=10) - abs_sp = df["shadow_price"].abs().to_list() - assert abs_sp == sorted(abs_sp, reverse=True) + df = get_enzyme_bottlenecks(model, top=999) + rounded = df["objective_gain"].round(9).to_list() + assert rounded == sorted(rounded, reverse=True) + assert (df["objective_gain"] >= 0).all() def test_columns(): model = _ectestgem_solvable() df = get_enzyme_bottlenecks(model, top=3) assert list(df.columns) == [ - "gene", "shadow_price", "flux", "cap_usage", "upper_bound", + "gene", "objective_gain", "shadow_price", "flux", "cap_usage", + "upper_bound", "multiple_active", ] assert df.index.name == "uniprot" + assert "n_resolved" in df.attrs + assert "active_constraints" in df.attrs + + +def _toy(solver="glpk", *, ubs=(4.0, 1.0, 5.0, 3.0), pool_ub=1e6): + """Four enzymes. E1 and E2 are parallel steps A -> B, E3 turns B into C, E4 turns A into D. + The objective is SK_C + SK_D, the substrate is not limiting and R4 is capped at 1000. + + With the default bounds E1, E2 and E3 are at their bound (4 + 1 = 5) and raising any one of + them alone leaves the optimum unchanged, while E4 is at its bound and raising it adds flux. + """ + model = EcModel("toy") + mets = {n: cobra.Metabolite(n, compartment="c") for n in + ["A_e", "A", "B", "C", "D", "pool", "prot_E1", "prot_E2", "prot_E3", "prot_E4"]} + model.add_metabolites(list(mets.values())) + + def rxn(rid, stoich, lb, ub): + r = cobra.Reaction(rid) + r.add_metabolites({mets[k]: v for k, v in stoich.items()}) + r.bounds = (lb, ub) + return r + + reactions = [ + rxn("EX_A", {"A_e": -1.0}, -1e6, 0.0), + rxn("TR_A", {"A_e": -1.0, "A": 1.0}, 0.0, 1e6), + rxn("R1", {"A": -1.0, "B": 1.0, "prot_E1": -1.0}, 0.0, 1000.0), + rxn("R2", {"A": -1.0, "B": 1.0, "prot_E2": -1.0}, 0.0, 1000.0), + rxn("R3", {"B": -1.0, "C": 1.0, "prot_E3": -1.0}, 0.0, 1000.0), + rxn("R4", {"A": -1.0, "D": 1.0, "prot_E4": -1.0}, 0.0, 1000.0), + rxn("SK_C", {"C": -1.0}, 0.0, 1000.0), + rxn("SK_D", {"D": -1.0}, 0.0, 1000.0), + rxn("prot_pool_exchange", {"pool": 1.0}, 0.0, pool_ub), + ] + for i, ub in enumerate(ubs, start=1): + reactions.append(rxn(f"usage_prot_E{i}", {"pool": -1.0, f"prot_E{i}": 1.0}, 0.0, ub)) + model.add_reactions(reactions) + model.objective = {model.reactions.SK_C: 1.0, model.reactions.SK_D: 1.0} + model.ec = EcData( + rxns=["R1", "R2", "R3", "R4"], + kcat=np.ones(4), + source=[""] * 4, + notes=[""] * 4, + eccodes=[""] * 4, + genes=["g1", "g2", "g3", "g4"], + enzymes=["E1", "E2", "E3", "E4"], + mw=np.full(4, 100.0), + sequence=[""] * 4, + concs=np.array(ubs, dtype=float), + rxn_enz_mat=sparse.identity(4, format="csr"), + ) + try: + model.solver = solver + except Exception as exc: + pytest.skip(f"solver {solver} unavailable: {exc}") + return model + + +def test_toy_ranks_by_resolved_gain_not_by_dual(): + model = _toy() + df = get_enzyme_bottlenecks(model, top=10) + assert list(df.index) == ["E4", "E1", "E2", "E3"] + assert df.loc["E4", "objective_gain"] == pytest.approx(3.0) + assert df.loc[["E1", "E2", "E3"], "objective_gain"].to_list() == [0.0, 0.0, 0.0] + assert df.attrs["n_resolved"] == 4 + assert df.attrs["active_constraints"] == 4 + assert df["multiple_active"].all() + assert ((df["cap_usage"] - 1.0).abs() < 1e-6).all() + # A zero gain with a nonzero dual: the dual is reported but does not enter the order. + assert (df.loc[["E1", "E2", "E3"], "shadow_price"].abs() > 1e-9).any() + # The model is restored. + assert model.reactions.usage_prot_E4.upper_bound == 3.0 + assert model.enzymes.get_by_id("E4").flux == pytest.approx(3.0) + + +def test_non_binding_enzymes_are_not_resolved(): + model = _toy(ubs=(4.0, 3.0, 20.0, 5000.0)) + df = get_enzyme_bottlenecks(model, top=10) + # E1 and E2 carry 7 through E3 (ub 20) and are at their bound; E3 is not, and E4 is limited by R4. + assert df.attrs["n_resolved"] == 2 + assert df.attrs["active_constraints"] == 2 + assert df.loc[["E3", "E4"], "objective_gain"].to_list() == [0.0, 0.0] + assert df.loc["E1", "objective_gain"] > 0 + assert df.loc["E1", "multiple_active"] + assert list(df.index[:2]) == ["E1", "E2"] + + +def test_single_active_constraint_is_not_multiple_active(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0)) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["n_resolved"] == 1 + assert df.attrs["active_constraints"] == 1 + assert df.loc["E3", "objective_gain"] == pytest.approx(2.0) + assert not df.loc["E3", "multiple_active"] + assert df.loc[["E1", "E2", "E4"], "multiple_active"].all() + assert list(df.index)[0] == "E3" + + +def test_pool_exchange_at_its_bound_is_an_active_constraint(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0), pool_ub=2.0) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["active_constraints"] == 1 + assert df.attrs["n_resolved"] == 0 + assert df["multiple_active"].all() + assert (df["objective_gain"] == 0).all() + + +def test_pool_exchange_and_an_enzyme_at_their_bounds(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0), pool_ub=1004.0) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["active_constraints"] == 2 + assert df.attrs["n_resolved"] == 1 + assert df.loc["E3", "multiple_active"] + assert df.loc["E3", "objective_gain"] == 0.0 + + +def test_fold_change_is_respected(): + twice = get_enzyme_bottlenecks(_toy(), top=1, fold_change=2.0) + thrice = get_enzyme_bottlenecks(_toy(), top=1, fold_change=3.0) + assert twice.index[0] == thrice.index[0] == "E4" + assert twice["objective_gain"].iloc[0] == pytest.approx(3.0) + assert thrice["objective_gain"].iloc[0] == pytest.approx(6.0) + + +@pytest.mark.parametrize("solver", ["glpk", "scipy", "gurobi"]) +def test_order_does_not_depend_on_the_solver(solver): + reference = get_enzyme_bottlenecks(_toy("glpk"), top=10) + df = get_enzyme_bottlenecks(_toy(solver), top=10) + assert list(df.index) == list(reference.index) + assert df["objective_gain"].to_numpy() == pytest.approx( + reference["objective_gain"].to_numpy(), abs=1e-6 + ) + assert df["multiple_active"].to_list() == reference["multiple_active"].to_list() + assert df.attrs["n_resolved"] == reference.attrs["n_resolved"] def test_raises_on_infeasible(): diff --git a/tests/test_scipy_solver_support.py b/tests/test_scipy_solver_support.py index b0a155f..8e4c102 100644 --- a/tests/test_scipy_solver_support.py +++ b/tests/test_scipy_solver_support.py @@ -161,13 +161,15 @@ def test_get_enzyme_bottlenecks_on_scipy_equals_the_reference(reference, scipy_d scipy_model = _on("scipy", tight=True)[0] actual = get_enzyme_bottlenecks(scipy_model, top=999) assert (expected["shadow_price"].abs() > 1e-9).any() - assert list(actual.index[:1]) == list(expected.index[:1]) - assert set(actual.index) == set(expected.index) + assert list(actual.index) == list(expected.index) + assert actual["objective_gain"].to_numpy() == pytest.approx( + expected["objective_gain"].to_numpy(), abs=1e-6 + ) aligned = actual.loc[expected.index] assert aligned["shadow_price"].to_numpy() == pytest.approx( expected["shadow_price"].to_numpy(), abs=1e-6 ) - assert actual["shadow_price"].abs().is_monotonic_decreasing + assert actual["objective_gain"].round(9).is_monotonic_decreasing for uniprot, row in actual.iterrows(): assert row["flux"] == pytest.approx(scipy_model.enzymes.get_by_id(uniprot).flux) From dfa2172a2b2bfd24ab110fd6a531329239772bbe Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 18:17:11 +0100 Subject: [PATCH 2/4] get_enzyme_bottlenecks reports the gain in the objective's direction and requires a fold_change above 1 --- src/geckopy/utilities/bottlenecks.py | 10 +++++++--- tests/test_bottlenecks.py | 16 ++++++++++++++++ 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 0c31ead..6695438 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -32,7 +32,7 @@ def get_enzyme_bottlenecks( Solves the model. An enzyme is at its bound when its upper bound is positive and ``|flux| / upper_bound`` is at least ``1 - 1e-6``. For each such enzyme the upper bound of its usage reaction is raised to ``upper_bound * fold_change``, the model is solved again and ``objective_gain`` is the new objective minus the objective of the first solve (0 when the new solve is not optimal or the change is below 1e-10). The model is restored after each solve. Every other enzyme has ``objective_gain`` 0 without a solve, because raising a bound that does not bind cannot change the optimum. The change is signed in the direction of the objective as the solver returns it, so it is positive for a maximised objective that increases. - The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. + ``objective_gain`` is the improvement of the objective in its own direction: a rise for a maximised objective and a fall for a minimised one. The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. ``multiple_active`` is True for an enzyme when a capacity constraint other than its own is active at the same optimum, which is where the dual says little about the enzyme. The capacity constraints counted are the enzymes at their bound and the protein pool exchange (``prot_pool_exchange``) when its flux is within 1e-6 relative of a finite, positive upper bound. @@ -45,7 +45,7 @@ def get_enzyme_bottlenecks( top Number of top-ranked enzymes to return. Default 10. fold_change - Factor by which the upper bound of an enzyme at its bound is raised. Default 2.0. + Factor, greater than 1, by which the upper bound of an enzyme at its bound is raised. Default 2.0. Returns ------- @@ -54,9 +54,13 @@ def get_enzyme_bottlenecks( Raises ------ + ValueError + If ``fold_change`` is not greater than 1. RuntimeError If the solve status is not optimal. """ + if not fold_change > 1: + raise ValueError(f"fold_change must be greater than 1, got {fold_change!r}") sol = optimize_model(model) if sol.status != "optimal": raise RuntimeError( @@ -95,7 +99,7 @@ def get_enzyme_bottlenecks( usage.upper_bound = rows[i]["upper_bound"] * fold_change new = optimize_model(model) if new.status == "optimal": - gain = new.objective_value - base + gain = (new.objective_value - base) * (1.0 if model.objective_direction == "max" else -1.0) rows[i]["objective_gain"] = gain if gain >= _GAIN_TOL else 0.0 if at_cap: optimize_model(model) diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index d8a2246..59761c9 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -186,6 +186,22 @@ def test_fold_change_is_respected(): assert thrice["objective_gain"].iloc[0] == pytest.approx(6.0) +@pytest.mark.parametrize("fold_change", [1.0, 0.5, float("nan")]) +def test_fold_change_must_exceed_one(fold_change): + with pytest.raises(ValueError, match="fold_change"): + get_enzyme_bottlenecks(_toy(), fold_change=fold_change) + + +def test_gain_is_the_improvement_of_a_minimised_objective(): + maximised = get_enzyme_bottlenecks(_toy(), top=4) + model = _toy() + model.objective = {model.reactions.get_by_id("SK_C"): -1.0, model.reactions.get_by_id("SK_D"): -1.0} + model.objective_direction = "min" + minimised = get_enzyme_bottlenecks(model, top=4) + assert list(minimised.index) == list(maximised.index) + assert minimised["objective_gain"].to_numpy() == pytest.approx(maximised["objective_gain"].to_numpy()) + + @pytest.mark.parametrize("solver", ["glpk", "scipy", "gurobi"]) def test_order_does_not_depend_on_the_solver(solver): reference = get_enzyme_bottlenecks(_toy("glpk"), top=10) From a2d533795b20da0b40a5d3f6583ca10477d9049c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:13:12 +0100 Subject: [PATCH 3/4] get_enzyme_bottlenecks logs a warning when no enzyme is at its usage bound --- CHANGELOG.md | 2 +- src/geckopy/utilities/bottlenecks.py | 15 +++++++++++++++ tests/test_bottlenecks.py | 25 +++++++++++++++++++++++++ 3 files changed, 41 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a01a3d..6abf96d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ ## Unreleased -* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. +* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. When no enzyme is at its usage bound, as on ecYeastGEM where the protein pool is the active constraint, a warning says so and the rows are in model order. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 6695438..d7d96da 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -8,6 +8,7 @@ """ from __future__ import annotations +import logging import math from typing import TYPE_CHECKING @@ -19,6 +20,8 @@ if TYPE_CHECKING: from ..ec_model.ec_model import EcModel +logger = logging.getLogger(__name__) + # An enzyme is at its bound when |flux| / upper_bound is at least 1 - _AT_CAP_TOL. _AT_CAP_TOL = 1e-6 # Changes of the objective below this value are set to 0. @@ -52,6 +55,10 @@ def get_enzyme_bottlenecks( pandas.DataFrame Indexed by uniprot id, columns ``[gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active]``, truncated to ``top`` rows. ``df.attrs["n_resolved"]`` is the number of enzymes that were solved again and ``df.attrs["active_constraints"]`` the number of active capacity constraints. + Warns + ----- + When no enzyme is at its usage bound a warning is logged: no single enzyme limits the objective, every ``objective_gain`` is 0 and the rows are in model order. + Raises ------ ValueError @@ -87,11 +94,19 @@ def get_enzyme_bottlenecks( if row["upper_bound"] > 0 and row["cap_usage"] >= 1 - _AT_CAP_TOL ] active = len(at_cap) + pool_active = False if model.reactions.has_id(POOL_EXCHANGE_ID): pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) pool_ub = float(pool.upper_bound) if 0 < pool_ub < math.inf and abs(float(pool.flux)) >= pool_ub * (1 - _AT_CAP_TOL): active += 1 + pool_active = True + if not at_cap: + logger.warning( + "get_enzyme_bottlenecks: no enzyme is at its usage bound, so no single enzyme limits " + "the objective; %s. Every objective_gain is 0 and the rows are in model order.", + "the protein pool exchange is at its bound" if pool_active else "no capacity constraint is active", + ) for i in at_cap: usage = enzymes[i].usage_reaction diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index 59761c9..6a4e989 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -221,3 +221,28 @@ def test_raises_on_infeasible(): r3.lower_bound = 1e9 with pytest.raises(RuntimeError, match="Solver status"): get_enzyme_bottlenecks(model, top=3) + + +def test_warns_when_no_enzyme_is_at_its_bound(caplog): + model = _toy(ubs=(1e5, 1e5, 1e5, 1e5)) + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + df = get_enzyme_bottlenecks(model) + assert df.attrs["n_resolved"] == 0 + assert (df["objective_gain"] == 0).all() + assert "no enzyme is at its usage bound" in caplog.text + assert "no capacity constraint is active" in caplog.text + + +def test_warning_names_the_protein_pool_when_it_is_the_active_constraint(caplog): + model = _toy(ubs=(1e5, 1e5, 1e5, 1e5), pool_ub=6.0) + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + df = get_enzyme_bottlenecks(model) + assert df.attrs["n_resolved"] == 0 + assert df.attrs["active_constraints"] == 1 + assert "the protein pool exchange is at its bound" in caplog.text + + +def test_no_warning_when_an_enzyme_is_at_its_bound(caplog): + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + get_enzyme_bottlenecks(_toy()) + assert "no enzyme is at its usage bound" not in caplog.text From d7bfb7f709c76c5b067fabb167ae2d92d6ae5106 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:13:34 +0100 Subject: [PATCH 4/4] Migration notes: no dash in the get_enzyme_bottlenecks bullet --- docs/migrating_from_gecko_matlab.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 5ea6eba..ae8128d 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -456,7 +456,7 @@ for the MATLAB-side items the port surfaced. ### geckopy-only (no MATLAB counterpart) - `pfba_enzymes` — enzyme-pool-minimising pFBA. -- `get_enzyme_bottlenecks` — bottleneck ranking by the re-solved change of the objective. +- `get_enzyme_bottlenecks`: bottleneck ranking by the re-solved change of the objective. - `relax_proteomics_greedy` — greedy alternative to `flexibilize_enz_concs`. - `merge_kcats` — generalised n-source kcat merge (replaces the