enzyme_usage: demand ratio of the usage over the measured concentrations

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+* `enzyme_usage` reports `measured_protein`, `demanded_protein` and `demand_ratio`, the usage of the enzymes with a measured concentration divided by their measured concentrations, and the same restricted to enzymes with a usage above 1e-7 mg/gDCW (`measured_protein_active`, `demand_ratio_active`). `report_enzyme_usage` carries `demand_ratio` and `demand_ratio_active`. GECKO's `getEnzymeUsage` and `reportEnzymeUsage` give the same values.
 * geckopy is licensed under the GNU General Public License v3 (GPL-3.0), as RAVEN and raven-toolbox. `LICENSE`, the `license` field and the classifier of `pyproject.toml`, and the README state it. Releases up to 4.0.0b2 were published under the MIT licence.
 
 ## 4.0.0b2 — 2026-11-02

src/geckopy/utilities/enzyme_usage.py

@@ -20,6 +20,8 @@
 
 logger = logging.getLogger(__name__)
 
+_ACTIVE_USAGE = 1e-7
+
 
 @dataclass
 class EnzymeUsageResult:
@@ -41,6 +43,30 @@ class EnzymeUsageResult:
     fluxes
         The original input flux distribution (carried through for
         downstream use, e.g. by ``report_enzyme_usage``).
+    measured_protein
+        Sum of the measured concentrations in ``model.ec.concs`` (mg/gDCW)
+        over the enzymes with a non-NaN entry.
+    demanded_protein
+        Sum of ``abs_usage`` over the same enzymes (mg/gDCW).
+    demand_ratio
+        ``demanded_protein / measured_protein``. NaN when no concentration
+        is measured; Inf when the measured concentrations sum to 0 and the
+        demand is positive. Above 1, the solution uses more of the measured
+        enzymes than was measured, which happens when the upper bounds of
+        the usage reactions are above the measurements (for instance after
+        ``flexibilize_enz_concs``). A ratio at or below 1 does not exclude
+        that one enzyme exceeds its measurement, since the sum is over all
+        measured enzymes; ``cap_usage`` gives the ratio per enzyme. Enzymes
+        without a measured concentration draw from the protein pool and are
+        not counted.
+    measured_protein_active, demand_ratio_active
+        The same two values restricted to the enzymes with an
+        ``abs_usage`` above 1e-7 mg/gDCW. Measured enzymes that are inactive
+        in the solution raise ``measured_protein`` and can hide a conflict
+        among the active ones.
+
+    The five demand fields are computed over all enzymes with a usage
+    reaction, whatever the value of ``include_zero``.
     """
 
     prot_id: list[str] = field(default_factory=list)
@@ -54,6 +80,11 @@ class EnzymeUsageResult:
         default_factory=lambda: np.empty(0, dtype=float)
     )
     fluxes: object = None
+    measured_protein: float = 0.0
+    measured_protein_active: float = 0.0
+    demanded_protein: float = 0.0
+    demand_ratio: float = float("nan")
+    demand_ratio_active: float = float("nan")
 
 
 def enzyme_usage(
@@ -103,9 +134,11 @@ def enzyme_usage(
     prot_ids: list[str] = []
     abs_usage_list: list[float] = []
     ub_list: list[float] = []
+    concs_list: list[float] = []
+    concs = np.asarray(model.ec.concs, dtype=float).ravel()
 
     left_out: list[str] = []
-    for enzyme in model.ec.enzymes:
+    for index, enzyme in enumerate(model.ec.enzymes):
         rxn_id = f"{USAGE_PREFIX}{enzyme}"
         if rxn_id not in cobra_rxn_ids:
             left_out.append(enzyme)
@@ -115,6 +148,7 @@ def enzyme_usage(
         prot_ids.append(enzyme)
         abs_usage_list.append(abs(flux))
         ub_list.append(rxn_ub)
+        concs_list.append(float(concs[index]) if index < concs.size else np.nan)
 
     if left_out:
         logger.warning(
@@ -129,6 +163,14 @@ def enzyme_usage(
             ub_arr != 0, abs_usage_arr / ub_arr, np.nan,
         )
 
+    concs_arr = np.array(concs_list, dtype=float)
+    measured = ~np.isnan(concs_arr)
+    active = abs_usage_arr > _ACTIVE_USAGE
+    measured_protein = float(concs_arr[measured].sum())
+    measured_active = float(concs_arr[measured & active].sum())
+    demanded = float(abs_usage_arr[measured].sum())
+    demanded_active = float(abs_usage_arr[measured & active].sum())
+
     if not include_zero:
         keep = abs_usage_arr > 0
         prot_ids = [p for p, k in zip(prot_ids, keep) if k]
@@ -142,9 +184,23 @@ def enzyme_usage(
         cap_usage=cap_arr,
         ub=ub_arr,
         fluxes=fluxes,
+        measured_protein=measured_protein,
+        measured_protein_active=measured_active,
+        demanded_protein=demanded,
+        demand_ratio=_ratio(demanded, measured_protein),
+        demand_ratio_active=_ratio(demanded_active, measured_active),
     )
 
 
+def _ratio(numerator: float, denominator: float) -> float:
+    """``numerator / denominator``: NaN for 0/0 and Inf for x/0, as in ``getEnzymeUsage``."""
+    if denominator == 0:
+        if numerator == 0 or np.isnan(numerator):
+            return float("nan")
+        return float("inf")
+    return numerator / denominator
+
+
 def _lookup_flux(
     fluxes: Union[pd.Series, Mapping[str, float]], rxn_id: str,
 ) -> float:

src/geckopy/utilities/report_enzyme_usage.py

@@ -56,6 +56,14 @@ class EnzymeUsageReport:
     total_usage_flux
         The protein pool exchange capacity (mg/gDCW) used as the
         denominator for ``perc_usage``.
+    demand_ratio
+        ``demand_ratio`` of the usage data: the protein used of the enzymes
+        with a measured concentration divided by their measured
+        concentrations. Above 1, the solution needs more of those enzymes
+        than was measured.
+    demand_ratio_active
+        ``demand_ratio_active`` of the usage data: the same ratio over the
+        enzymes that carry flux in the solution.
     """
 
     high_cap_usage: pd.DataFrame = field(
@@ -65,6 +73,8 @@ class EnzymeUsageReport:
         default_factory=lambda: pd.DataFrame(columns=_TOP_ABS_COLUMNS)
     )
     total_usage_flux: float = 0.0
+    demand_ratio: float = float("nan")
+    demand_ratio_active: float = float("nan")
 
 
 def report_enzyme_usage(
@@ -229,6 +239,8 @@ def report_enzyme_usage(
         high_cap_usage=pd.DataFrame(high_rows, columns=_HIGH_CAP_COLUMNS),
         top_abs_usage=pd.DataFrame(top_rows, columns=_TOP_ABS_COLUMNS),
         total_usage_flux=total_pool,
+        demand_ratio=usage_data.demand_ratio,
+        demand_ratio_active=usage_data.demand_ratio_active,
     )
 
 

tests/test_enzyme_usage.py

@@ -202,3 +202,126 @@ def test_enzyme_without_usage_reaction_is_left_out_with_a_warning(caplog):
         result = enzyme_usage(model, {"usage_prot_E": 3.0})
     assert result.prot_id == ["E"]
     assert "no usage reaction for the enzymes F" in caplog.text
+
+
+# --------------------------------------------------------------------------- #
+# Demand ratio
+# --------------------------------------------------------------------------- #
+
+def test_demand_ratio_at_the_measurement():
+    """The enzyme runs at its measured concentration: ratio 1."""
+    model = _build_toy(enzyme_ub=5.0)
+    result = enzyme_usage(model, model.optimize().fluxes)
+    assert result.measured_protein == pytest.approx(5.0)
+    assert result.demanded_protein == pytest.approx(5.0)
+    assert result.demand_ratio == pytest.approx(1.0)
+    assert result.demand_ratio_active == pytest.approx(1.0)
+
+
+def test_demand_ratio_above_one_when_the_bound_is_relaxed():
+    """A usage bound above the measurement lets the solution use more of the
+    enzyme than was measured."""
+    model = _build_toy(enzyme_ub=5.0)
+    model.reactions.get_by_id("usage_prot_E").upper_bound = 10.0
+    result = enzyme_usage(model, model.optimize().fluxes)
+    assert result.abs_usage[0] == pytest.approx(10.0)
+    assert result.measured_protein == pytest.approx(5.0)
+    assert result.demand_ratio == pytest.approx(2.0)
+    assert result.demand_ratio_active == pytest.approx(2.0)
+
+
+def test_demand_ratio_of_an_inactive_measured_enzyme():
+    """Without flux the demand is 0, and no active enzyme is measured."""
+    model = _build_toy(enzyme_ub=5.0)
+    model.reactions.get_by_id("EX_A").lower_bound = 0.0
+    result = enzyme_usage(model, model.optimize().fluxes)
+    assert result.measured_protein == pytest.approx(5.0)
+    assert result.demand_ratio == pytest.approx(0.0, abs=1e-9)
+    assert result.measured_protein_active == 0.0
+    assert np.isnan(result.demand_ratio_active)
+
+
+def test_demand_ratio_is_nan_without_a_measured_concentration():
+    """An enzyme with a NaN concentration draws from the protein pool and is
+    not counted."""
+    model = _build_toy(enzyme_ub=5.0)
+    model.ec.concs[:] = np.nan
+    result = enzyme_usage(model, model.optimize().fluxes)
+    assert result.measured_protein == 0.0
+    assert result.demanded_protein == 0.0
+    assert np.isnan(result.demand_ratio)
+    assert np.isnan(result.demand_ratio_active)
+
+
+def test_demand_ratio_does_not_depend_on_include_zero():
+    model, fluxes = _build_mixed()
+    fluxes["usage_prot_F"] = 0.0
+    with_zero = enzyme_usage(model, fluxes)
+    without_zero = enzyme_usage(model, fluxes, include_zero=False)
+    assert with_zero.prot_id == ["E", "F", "G"]
+    assert without_zero.prot_id == ["E", "G"]
+    for name in (
+        "measured_protein", "measured_protein_active", "demanded_protein",
+        "demand_ratio", "demand_ratio_active",
+    ):
+        assert getattr(without_zero, name) == getattr(with_zero, name)
+    assert without_zero.measured_protein == pytest.approx(8.0)
+    assert without_zero.demand_ratio == pytest.approx(0.25)
+
+
+def test_report_carries_the_demand_ratios():
+    from geckopy.utilities import report_enzyme_usage
+
+    model = _build_toy(enzyme_ub=5.0)
+    model.reactions.get_by_id("usage_prot_E").upper_bound = 10.0
+    result = enzyme_usage(model, model.optimize().fluxes)
+    report = report_enzyme_usage(model, result)
+    assert report.demand_ratio == pytest.approx(2.0)
+    assert report.demand_ratio_active == pytest.approx(2.0)
+
+
+def _build_mixed():
+    """Enzymes E (measured 5), H (measured 4, no usage reaction), F (measured
+    3) and G (not measured), with usage fluxes 2, 1e-8 and 7."""
+    model = _build_toy(enzyme_ub=5.0)
+    for name in ("F", "G"):
+        usage = cobra.Reaction(f"usage_prot_{name}")
+        usage.upper_bound = 10.0
+        model.add_reactions([usage])
+    model.ec.enzymes[:] = ["E", "H", "F", "G"]
+    model.ec.concs = np.array([5.0, 4.0, 3.0, np.nan])
+    fluxes = {"usage_prot_E": 2.0, "usage_prot_F": 1e-8, "usage_prot_G": 7.0}
+    return model, fluxes
+
+
+def test_demand_ratio_over_measured_and_active_enzymes():
+    model, fluxes = _build_mixed()
+    result = enzyme_usage(model, fluxes)
+    assert result.prot_id == ["E", "F", "G"]
+    assert result.measured_protein == pytest.approx(8.0)
+    assert result.demanded_protein == pytest.approx(2.0 + 1e-8)
+    assert result.demand_ratio == pytest.approx((2.0 + 1e-8) / 8.0)
+    assert result.measured_protein_active == pytest.approx(5.0)
+    assert result.demand_ratio_active == pytest.approx(0.4)
+
+
+def test_demand_ratio_with_a_short_concentration_vector():
+    model, fluxes = _build_mixed()
+    model.ec.concs = np.array([5.0])
+    result = enzyme_usage(model, fluxes)
+    assert result.measured_protein == pytest.approx(5.0)
+    assert result.demand_ratio_active == pytest.approx(0.4)
+    model.ec.concs = np.array([])
+    result = enzyme_usage(model, fluxes)
+    assert np.isnan(result.demand_ratio)
+
+
+def test_report_carries_nan_demand_ratios():
+    from geckopy.utilities import report_enzyme_usage
+
+    model = _build_toy(enzyme_ub=5.0)
+    model.ec.concs[:] = np.nan
+    result = enzyme_usage(model, model.optimize().fluxes)
+    report = report_enzyme_usage(model, result)
+    assert np.isnan(report.demand_ratio)
+    assert np.isnan(report.demand_ratio_active)