read_yaml_model: flip usage and pool reactions by the sign of the pool coefficient

CHANGELOG.md

@@ -6,6 +6,7 @@ Milestones in the raven-toolbox port. For function-level status see
 
 ## Unreleased
 
+* `read_yaml_model` turns a `usage_prot_*` reaction forward when its `prot_pool` coefficient is positive and `prot_pool_exchange` when it is negative, as GECKO's `upgradeEcModel` does, and changes the sign of the objective coefficient of a reaction it turns (was: a negative lower bound decided, and the coefficient stayed, which reversed the objective of the model).
 * `write_yaml_model` writes an ec reaction without enzymes with no `enzymes` entry (was: `enzymes: !!omap []`). `read_yaml_model` reads a missing `enzymes` entry, a content-less `enzymes: !!omap` and `enzymes: !!omap []` as an empty enzyme map, so files written by earlier versions of GECKO, RAVEN and raven-toolbox still load.
 * `init.fill_tasks` with `resolve_ties=True` on Gurobi restores `MIPGap` and `MIPGapAbs` of `reference_model` after each gap-fill (was: 0.0 and 0.4 stayed set).
 * `tasks.find_task_essential_reactions` writes its `cache_path` as JSON (was: pickle) and ignores, with a warning, an existing file that is not a JSON cache of the current version, such as a pickle. `checkpoint_every` below 1 raises `ValueError` (was: a write after every task).

src/raven_toolbox/io/yaml.py

@@ -34,8 +34,8 @@
 * very old RAVEN files written as a bare ``-`` sequence of single-key mappings
   rather than one big mapping;
 * MATLAB GECKO ecModels whose ``usage_prot_*`` and ``prot_pool_exchange``
-  reactions use the older reverse-direction convention (negative lower bound,
-  swapped stoichiometry signs); these are flipped on load so consumers always
+  reactions use the older reverse-direction convention (swapped
+  stoichiometry signs, negative flux); these are flipped on load so consumers always
   see the forward convention.
 """
 from __future__ import annotations
@@ -597,19 +597,27 @@ def _flip_legacy_prot_direction(model: cobra.Model) -> None:
     The current convention (in both geckopy and recent MATLAB GECKO)
     treats them as ordinary forward reactions with positive flux. When a
     loaded model still uses the older convention we flip the affected
-    reactions in place so consumers never have to handle two shapes.
+    reactions in place so consumers never have to handle two shapes. The
+    objective coefficient of a flipped reaction changes sign with its flux.
 
-    The signature we look for is any ``usage_prot_*`` or
-    ``prot_pool_exchange`` reaction whose lower bound is negative.
+    The signature we look for is the sign of the ``prot_pool`` coefficient:
+    positive in a ``usage_prot_*`` reaction, negative in
+    ``prot_pool_exchange``. The bounds do not decide, so a reaction with
+    the current stoichiometry and a negative lower bound is left as it is.
     """
+    try:
+        pool = model.metabolites.get_by_id("prot_pool")
+    except KeyError:
+        return
     flipped: list[str] = []
     for rxn in model.reactions:
-        if not (
-            rxn.id.startswith("usage_prot_")
-            or rxn.id == "prot_pool_exchange"
-        ):
+        if rxn.id == "prot_pool_exchange":
+            legacy = rxn.metabolites.get(pool, 0.0) < 0
+        elif rxn.id.startswith("usage_prot_"):
+            legacy = rxn.metabolites.get(pool, 0.0) > 0
+        else:
             continue
-        if rxn.lower_bound >= -1e-9:
+        if not legacy:
             continue
         rxn.add_metabolites(
             {met: -2.0 * coef for met, coef in rxn.metabolites.items()},
@@ -618,6 +626,7 @@ def _flip_legacy_prot_direction(model: cobra.Model) -> None:
         lb, ub = rxn.lower_bound, rxn.upper_bound
         rxn.lower_bound = -ub
         rxn.upper_bound = -lb
+        rxn.objective_coefficient = -rxn.objective_coefficient
         flipped.append(rxn.id)
     if flipped:
         warnings.warn(

tests/test_io_yaml_ec_data.py

@@ -341,6 +341,7 @@ def test_legacy_reverse_direction_prot_flipped(tmp_path):
             "metabolites": {"prot_P1": -1, "prot_pool": 1},  # swapped signs
             "lower_bound": -1000.0,
             "upper_bound": 0.0,
+            "objective_coefficient": 2.0,
         }
     )
     with pytest.warns(UserWarning, match="reverse-direction"):
@@ -352,6 +353,27 @@ def test_legacy_reverse_direction_prot_flipped(tmp_path):
     coefs = {m.id: c for m, c in r.metabolites.items()}
     assert coefs["prot_P1"] == 1.0
     assert coefs["prot_pool"] == -1.0
+    # The objective coefficient changes sign with the flux.
+    assert r.objective_coefficient == -2.0
+
+
+def test_forward_prot_reaction_with_negative_lower_bound_not_flipped(tmp_path):
+    """The sign of the ``prot_pool`` coefficient decides, not the bounds."""
+    doc = _minimal_model_doc()
+    doc["metabolites"].append({"id": "prot_P1", "name": "prot", "compartment": "c"})
+    doc["metabolites"].append({"id": "prot_pool", "name": "pool", "compartment": "c"})
+    doc["reactions"].append(
+        {
+            "id": "usage_prot_P1",
+            "metabolites": {"prot_P1": 1, "prot_pool": -1},
+            "lower_bound": -1000.0,
+            "upper_bound": 1000.0,
+        }
+    )
+    model = read_yaml_model(_write_yaml(doc, tmp_path / "m.yml"))
+    r = model.reactions.get_by_id("usage_prot_P1")
+    assert (r.lower_bound, r.upper_bound) == (-1000.0, 1000.0)
+    assert {m.id: c for m, c in r.metabolites.items()}["prot_pool"] == -1.0
 
 
 def test_legacy_bare_sequence_root_merged(tmp_path):