fix: fraction_of_optimum of pfba_enzymes for a negative or minimised objective

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+* `pfba_enzymes` with `fraction_of_optimum` below 1 lets the objective fall by that fraction of the magnitude of its optimum, for a negative optimum and for a minimised objective as for a maximised positive one, as MATLAB's `getPfbaEnzymes` does; it raised `Infeasible` for these objectives.
 * `ec_fseof` ranks the targets by the score `|slope| * |correlation|` rounded to 9 significant digits, with scores that are equal after rounding in model order, as GECKO's `ecFSEOF` does; scores that differ by solver noise used to decide the order.
 * `geckopy.kcat_tuning.evotune` exports `best_parsimonious`, `n_changed`, `parsimony_frontier` and `isozyme_tie_map`, the helpers behind the parsimony sweep and the tying of isozyme copies, which GECKO has as `bestParsimonious`, `nChanged`, `parsimonyFrontier` and `isozymeTieMap`.
 * `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""<TAB>""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number.

src/geckopy/utilities/pfba_enzymes.py

@@ -48,9 +48,9 @@ def pfba_enzymes(
 ) -> cobra.Solution:
     """Parsimonious FBA minimising the L1 norm of ``usage_prot_*`` fluxes.
 
-    Sets the objective (if given), fixes it at
-    ``fraction_of_optimum * 0.999999 * optimal_value``, then minimises
-    ``sum(usage_prot_*.forward_variable)``. Usage reactions are
+    Sets the objective (if given), lets it fall from its optimum by
+    ``1 - fraction_of_optimum * 0.999999`` times the magnitude of the
+    optimum, then minimises ``sum(usage_prot_*.forward_variable)``. Usage reactions are
     forward-only in the GECKO 3 layout (``lb=0``); if a user has
     flipped one to allow reverse flux, its ``reverse_variable`` is
     included automatically.
@@ -69,13 +69,16 @@ def pfba_enzymes(
         enzyme usage. If ``None``, the model's current objective is
         used as-is.
     fraction_of_optimum
-        Fraction of the objective's optimal value to fix as a
+        Fraction of the objective's optimal value to keep as a
         constraint before minimising enzyme usage (default 1.0). The
         fraction is multiplied by 0.999999, the margin of RAVEN's
-        ``solveLP`` that MATLAB's ``getPfbaEnzymes`` applies, so the
-        default fixes the objective at 0.999999 of its optimum and
-        0.9 fixes it at 0.8999991. A lower value trades some
-        objective value for potentially lower enzyme usage.
+        ``solveLP`` that MATLAB's ``getPfbaEnzymes`` applies. The
+        objective may fall from its optimum by ``1 - fraction * 0.999999``
+        times the magnitude of the optimum, so the default keeps it within
+        1e-6 of the optimum and 0.9 within 10 %, for a negative optimum
+        (the exchange flux of a nutrient that is maximised) as for a
+        positive one. A lower value trades some objective value for
+        potentially lower enzyme usage.
 
     Returns
     -------
@@ -108,9 +111,16 @@ def pfba_enzymes(
                 f"Objective is already {_PFBA_OBJ_NAME}; "
                 "pfba_enzymes appears to be nested."
             )
-        cobra.util.fix_objective_as_constraint(
-            m, fraction=fraction_of_optimum * _PFBA_MARGIN,
-        )
+        optimum = m.slim_optimize(error_value=None)
+        if optimum is None:
+            raise cobra.exceptions.OptimizationError(
+                "pfba_enzymes: the model is infeasible for its objective."
+            )
+        # The objective may fall by (1 - fraction * margin) of the magnitude
+        # of its optimum, for a positive and a negative optimum alike.
+        slack = (1.0 - fraction_of_optimum * _PFBA_MARGIN) * abs(optimum)
+        bound = optimum - slack if m.objective.direction == "max" else optimum + slack
+        cobra.util.fix_objective_as_constraint(m, bound=bound)
 
         usage_rxns = [r for r in m.reactions if r.id.startswith(USAGE_PREFIX)]
         variables = []

tests/test_pfba_enzymes.py

@@ -112,3 +112,26 @@ def test_gecko_light_raises():
     model.ec.gecko_light = True
     with pytest.raises(NotImplementedError):
         pfba_enzymes(model)
+
+
+def test_fraction_of_a_negative_optimum_lets_the_objective_fall():
+    model = _ectestgem_ec_model_with_kcats()
+    # Maximising -R3 with R3 >= 1 gives an optimum of -1.
+    model.reactions.get_by_id("R3").lower_bound = 1.0
+    model.objective = {model.reactions.get_by_id("R3"): -1.0}
+    optimum = model.slim_optimize()
+    assert optimum == pytest.approx(-1.0)
+    sol = pfba_enzymes(model, fraction_of_optimum=0.9)
+    assert sol.status == "optimal"
+    assert -sol.fluxes["R3"] >= optimum * 1.1 - 1e-6
+
+
+def test_fraction_of_a_minimised_objective_lets_it_rise():
+    model = _ectestgem_ec_model_with_kcats()
+    model.reactions.get_by_id("R3").lower_bound = 1.0
+    model.objective = "R3"
+    model.objective_direction = "min"
+    optimum = model.slim_optimize()
+    sol = pfba_enzymes(model, fraction_of_optimum=0.9)
+    assert sol.status == "optimal"
+    assert sol.fluxes["R3"] <= optimum * 1.1 + 1e-6