get_enzyme_bottlenecks ranks by the re-solved change of the objective

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_<id>` 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_conc_control_coeffs` takes the control coefficient by finite difference by default, as GECKO's `getConcControlCoeffs` does: it raises one usage bound by `fold_change` (2), re-solves and divides the rise of the optimum by the rise of the bound. The reduced-cost path, which reads the dual of the usage reaction and is unreliable at a degenerate vertex, is `method="shadow_price"`. The finite-difference path needs one solve per protein, so each iteration of `flexibilize_enz_concs`, which calls it, takes minutes on a genome-scale ecModel.
 * `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.

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

src/geckopy/utilities/bottlenecks.py

@@ -1,69 +1,79 @@
-"""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_<id>`` 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_<id>`` 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 logging
+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
 
+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.
+_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.
+
+    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.
 
-def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame:
-    """Top-N enzymes by absolute shadow price.
+    ``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_<id>`` 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.
 
-    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]``.
+    ``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.
+    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
     ----------
     model
         EcModel to solve.
     top
         Number of top-ranked enzymes to return. Default 10.
+    fold_change
+        Factor, greater than 1, 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.
+
+    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
+        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(
             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 +82,54 @@ 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)
+    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
+        with model:
+            usage.upper_bound = rows[i]["upper_bound"] * fold_change
+            new = optimize_model(model)
+            if new.status == "optimal":
+                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)
+
+    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