From b69f7777ee25d814e0f9bc36e247ec436c591427 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 2 Nov 2026 23:42:17 +0100 Subject: [PATCH] enzyme_usage: demand ratio of the usage over the measured concentrations --- CHANGELOG.md | 1 + src/geckopy/utilities/enzyme_usage.py | 58 ++++++++- src/geckopy/utilities/report_enzyme_usage.py | 12 ++ tests/test_enzyme_usage.py | 123 +++++++++++++++++++ 4 files changed, 193 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fbcd3dc..be62228 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/geckopy/utilities/enzyme_usage.py b/src/geckopy/utilities/enzyme_usage.py index ac780dc..7080988 100644 --- a/src/geckopy/utilities/enzyme_usage.py +++ b/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: diff --git a/src/geckopy/utilities/report_enzyme_usage.py b/src/geckopy/utilities/report_enzyme_usage.py index 65af600..3b9b639 100644 --- a/src/geckopy/utilities/report_enzyme_usage.py +++ b/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, ) diff --git a/tests/test_enzyme_usage.py b/tests/test_enzyme_usage.py index 4b3869d..119071d 100644 --- a/tests/test_enzyme_usage.py +++ b/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)