From 0c8fb29ba85235f5531c2bb9e84814cf6abc1c00 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 14 Oct 2026 23:13:24 +0200 Subject: [PATCH 1/5] feat: level complex abundances and flexibilize per complex level_complex_concs (and level_complex_abundances for molar data) gives every enzyme complex one abundance, the median per copy over its measured subunits, so a single badly quantified subunit no longer sets the complex's capacity. complex_families groups isozyme variants into one complex and leaves subunits shared between complexes as measured. flexibilize_enz_concs(per_complex=True) raises all subunits of a limiting complex by the same factor. See raven-gecko-parity#131. --- src/geckopy/__init__.py | 10 + src/geckopy/limit_proteins/__init__.py | 12 + .../limit_proteins/flexibilize_enz_concs.py | 55 +++- .../limit_proteins/level_complex_concs.py | 243 ++++++++++++++++++ tests/test_level_complex_concs.py | 171 ++++++++++++ 5 files changed, 486 insertions(+), 5 deletions(-) create mode 100644 src/geckopy/limit_proteins/level_complex_concs.py create mode 100644 tests/test_level_complex_concs.py diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index 90d3e8a..fd1b166 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -111,16 +111,21 @@ # limit_proteins from geckopy.limit_proteins import ( + ComplexFamilies, FlexEnzResult, GreedyRelaxResult, + LevelComplexResult, RelaxationStep, apply_flux_data_constraints, calculate_f_factor, + complex_families, constrain_enz_concs, constrain_flux_data, # deprecated; alias of apply_flux_data_constraints fill_enz_concs, flexibilize_enz_concs, get_conc_control_coeffs, + level_complex_abundances, + level_complex_concs, relax_proteomics_greedy, ) @@ -206,15 +211,20 @@ "sensitivity_tuning", "truncate_values", # limit_proteins + "ComplexFamilies", "FlexEnzResult", "GreedyRelaxResult", + "LevelComplexResult", "RelaxationStep", "apply_flux_data_constraints", "calculate_f_factor", + "complex_families", "constrain_enz_concs", "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "level_complex_abundances", + "level_complex_concs", "relax_proteomics_greedy", # utilities "AddNewRxnsResult", diff --git a/src/geckopy/limit_proteins/__init__.py b/src/geckopy/limit_proteins/__init__.py index 79f5f1a..f7c5374 100644 --- a/src/geckopy/limit_proteins/__init__.py +++ b/src/geckopy/limit_proteins/__init__.py @@ -11,6 +11,13 @@ from .fill_enz_concs import fill_enz_concs from .flexibilize_enz_concs import FlexEnzResult, flexibilize_enz_concs from .get_conc_control_coeffs import get_conc_control_coeffs +from .level_complex_concs import ( + ComplexFamilies, + LevelComplexResult, + complex_families, + level_complex_abundances, + level_complex_concs, +) from .relax_proteomics_greedy import ( GreedyRelaxResult, RelaxationStep, @@ -18,15 +25,20 @@ ) __all__ = [ + "ComplexFamilies", "FlexEnzResult", + "LevelComplexResult", "GreedyRelaxResult", "RelaxationStep", "apply_flux_data_constraints", "calculate_f_factor", + "complex_families", "constrain_enz_concs", "constrain_flux_data", # deprecated; alias of apply_flux_data_constraints "fill_enz_concs", "flexibilize_enz_concs", "get_conc_control_coeffs", + "level_complex_abundances", + "level_complex_concs", "relax_proteomics_greedy", ] diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 04e5426..797a60b 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -37,6 +37,7 @@ USAGE_PREFIX, ) from .get_conc_control_coeffs import get_conc_control_coeffs +from .level_complex_concs import ComplexFamilies, complex_families if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -103,6 +104,7 @@ def flexibilize_enz_concs( iter_per_enzyme: int = 5, bio_rxn: Optional[str] = None, verbose: bool = True, + per_complex: bool = False, ) -> FlexEnzResult: """Iteratively relax enzyme concentration constraints until the model reaches ``exp_growth``. @@ -153,6 +155,12 @@ def flexibilize_enz_concs( adapter; pass it (with ``exp_growth``) to run without an adapter. verbose Whether per-iteration progress is logged at INFO. + per_complex + When the most limiting enzyme is a non-shared subunit of a + complex (see ``complex_families``), raise every measured subunit + of that complex by the same factor, and keep one factor per + complex in the refinement pass, so that concentrations levelled + by ``level_complex_concs`` stay levelled. Off by default. Returns ------- @@ -203,6 +211,9 @@ def flexibilize_enz_concs( proteins = [model.ec.enzymes[i] for i in measured_idx] n_proteins = len(proteins) frequence = np.zeros(n_proteins, dtype=int) + touched = np.zeros(n_proteins, dtype=bool) + position = {p: i for i, p in enumerate(proteins)} + families: Optional[ComplexFamilies] = complex_families(model) if per_complex else None initial_sol = model.optimize() if ( @@ -265,11 +276,14 @@ def flexibilize_enz_concs( break increase = fold_change * frequence[max_idx] - ec_idx = measured_idx[max_idx] - usage_rxn = model.reactions.get_by_id( - f"{USAGE_PREFIX}{proteins[max_idx]}" - ) - usage_rxn.upper_bound = float(model.ec.concs[ec_idx]) * (1.0 + increase) + members = families.members(proteins[max_idx]) if families else [proteins[max_idx]] + for member in members: + if member not in position: + continue + i = position[member] + touched[i] = True + usage_rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{member}") + usage_rxn.upper_bound = float(model.ec.concs[measured_idx[i]]) * (1.0 + increase) new_sol = model.optimize() pred_growth = float(new_sol.objective_value or 0.0) @@ -282,6 +296,7 @@ def flexibilize_enz_concs( return _build_result( model, proteins, measured_idx, frequence, flex_break, exp_growth, pool_old, pool_new, bio_rxn_id, + touched=touched, families=families, ) @@ -345,9 +360,13 @@ def _build_result( pool_old: Optional[float], pool_new: Optional[float], bio_rxn_id: str, + touched: Optional[np.ndarray] = None, + families: Optional[ComplexFamilies] = None, ) -> FlexEnzResult: """Build the FlexEnzResult and apply the post-loop UB refinement.""" flex_mask = frequence > 0 + if touched is not None: + flex_mask = flex_mask | touched flex_indices = np.where(flex_mask)[0] if len(flex_indices) == 0: return _maybe_add_pool(FlexEnzResult(), pool_old, pool_new) @@ -377,6 +396,8 @@ def _build_result( dtype=float, ) + if families is not None: + new_concs = _per_complex(flex_proteins, old_concs_arr, new_concs, families) keep_old_mask = new_concs < old_concs_arr for i, rxn_id in enumerate(usage_rxn_ids): ub = ( @@ -412,6 +433,30 @@ def _build_result( return _finalize_and_sort(_maybe_add_pool(result, pool_old, pool_new)) +def _per_complex( + proteins: list[str], + old: np.ndarray, + new: np.ndarray, + families: ComplexFamilies, +) -> np.ndarray: + """One factor per complex: every flexibilized subunit takes the largest new/old of its complex.""" + out = new.copy() + index = {p: i for i, p in enumerate(proteins)} + seen: set[str] = set() + for protein in proteins: + if protein in seen: + continue + group = [m for m in families.members(protein) if m in index] + seen.update(group) + if len(group) < 2: + continue + ratios = [new[index[m]] / old[index[m]] for m in group if old[index[m]] > 0] + factor = max(ratios) if ratios else 1.0 + for m in group: + out[index[m]] = old[index[m]] * factor + return out + + def _maybe_add_pool( result: FlexEnzResult, pool_old: Optional[float], diff --git a/src/geckopy/limit_proteins/level_complex_concs.py b/src/geckopy/limit_proteins/level_complex_concs.py new file mode 100644 index 0000000..ae59f95 --- /dev/null +++ b/src/geckopy/limit_proteins/level_complex_concs.py @@ -0,0 +1,243 @@ +"""One abundance per enzyme complex. + +Absolute proteomics recovers the subunits of a complex unevenly, so +after ``fill_enz_concs`` the subunits of one complex can disagree several +fold about how much of it there is. Capped one by one, the complex can +then carry no more flux than its least abundant subunit allows. This +module groups the enzymes of every multi-subunit reaction into complexes +and gives each complex one abundance, scaled to every subunit by its +stoichiometry. + +A complex is a set of enzymes that catalyses a reaction together, taken +from ``ec.rxn_enz_mat``. Enzyme sets that share at least half of the +smaller set's subunits are one complex (isozyme variants that swap a +subunit, the same complex on different reactions). A subunit that ends +up in more than one complex, or that also catalyses a reaction on its +own, is *shared*: it is left as measured and does not take part in the +summary. +""" +from __future__ import annotations + +import logging +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Optional, Sequence + +import numpy as np + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + +_STATISTICS = {"median": np.median, "mean": np.mean} + + +@dataclass +class ComplexFamilies: + """The complexes of an ecModel and their shared subunits. + + Attributes + ---------- + complexes + One dict per complex, ``{enzyme: copies}`` over every subunit + that appears in any of its enzyme sets. + shared + Enzymes that belong to more than one complex, or to a complex + and also catalyse a reaction alone. + """ + + complexes: list[dict[str, float]] = field(default_factory=list) + shared: set[str] = field(default_factory=set) + + def members(self, enzyme: str) -> list[str]: + """The non-shared subunits of the complex ``enzyme`` belongs to, or ``[enzyme]``.""" + if enzyme in self.shared: + return [enzyme] + for complex_ in self.complexes: + if enzyme in complex_: + return [e for e in complex_ if e not in self.shared] + return [enzyme] + + +@dataclass +class LevelComplexResult: + """What ``level_complex_concs`` changed. + + Attributes + ---------- + adjusted + Measured subunits whose concentration was replaced by the complex's. + added + Unmeasured subunits given the complex's concentration. + shared + Subunits left as measured because they belong to more than one complex. + skipped + Complexes (as sorted subunit tuples) without any measured, non-shared subunit. + """ + + adjusted: list[str] = field(default_factory=list) + added: list[str] = field(default_factory=list) + shared: list[str] = field(default_factory=list) + skipped: list[tuple[str, ...]] = field(default_factory=list) + + +def complex_families(model: "EcModel", min_overlap: float = 0.5) -> ComplexFamilies: + """Group the enzymes of multi-subunit reactions into complexes. + + Parameters + ---------- + model + An EcModel with ``ec.rxn_enz_mat`` (copies per subunit). + min_overlap + Two enzyme sets are one complex when they share at least this + fraction of the smaller set's subunits. + """ + matrix = model.ec.rxn_enz_mat.tocsr() + enzymes = list(model.ec.enzymes) + sets: dict[frozenset, dict[str, float]] = {} + alone: set[str] = set() + for row in range(matrix.shape[0]): + start, end = matrix.indptr[row], matrix.indptr[row + 1] + present = {enzymes[j]: float(v) for j, v in zip(matrix.indices[start:end], matrix.data[start:end]) if v} + if len(present) == 1: + alone.update(present) + elif present: + sets.setdefault(frozenset(present), present) + + keys = list(sets) + parent = list(range(len(keys))) + + def find(i: int) -> int: + while parent[i] != i: + parent[i] = parent[parent[i]] + i = parent[i] + return i + + by_enzyme: dict[str, list[int]] = {} + for i, key in enumerate(keys): + for enzyme in key: + by_enzyme.setdefault(enzyme, []).append(i) + for indices in by_enzyme.values(): + for x in range(len(indices)): + for y in range(x + 1, len(indices)): + i, j = indices[x], indices[y] + if len(keys[i] & keys[j]) >= min_overlap * min(len(keys[i]), len(keys[j])): + parent[find(i)] = find(j) + + grouped: dict[int, dict[str, float]] = {} + for i, key in enumerate(keys): + complex_ = grouped.setdefault(find(i), {}) + for enzyme, copies in sets[key].items(): + complex_.setdefault(enzyme, copies) + complexes = list(grouped.values()) + + count: dict[str, int] = {} + for complex_ in complexes: + for enzyme in complex_: + count[enzyme] = count.get(enzyme, 0) + 1 + shared = {e for e, n in count.items() if n > 1} | (alone & set(count)) + return ComplexFamilies(complexes=complexes, shared=shared) + + +def level_complex_abundances( + model: "EcModel", + uniprot_ids: Sequence[str], + abundances: Sequence[float], + *, + statistic: str = "median", + fill_missing: bool = True, + families: Optional[ComplexFamilies] = None, +) -> tuple[list[str], np.ndarray, LevelComplexResult]: + """Level molar abundances to one abundance per complex. + + Parameters + ---------- + model + The EcModel whose complexes are used. + uniprot_ids, abundances + Measured molar abundances (any molar unit; per subunit, not per complex). + statistic + ``"median"`` or ``"mean"`` over the measured, non-shared subunits' + abundance divided by their copies. + fill_missing + Give unmeasured, non-shared subunits the complex's abundance. + families + Precomputed ``complex_families(model)``. + + Returns + ------- + ids, abundances, result + The levelled data (unmeasured subunits appended when filled) and + what changed. + """ + if statistic not in _STATISTICS: + raise ValueError(f"statistic must be one of {sorted(_STATISTICS)}, not {statistic!r}") + summarise = _STATISTICS[statistic] + families = families or complex_families(model) + values = {p: float(a) for p, a in zip(uniprot_ids, abundances)} + order = list(uniprot_ids) + result = LevelComplexResult(shared=sorted(families.shared & set(values))) + + for complex_ in families.complexes: + members = {e: c for e, c in complex_.items() if e not in families.shared} + measured = [values[e] / c for e, c in members.items() if np.isfinite(values.get(e, np.nan)) and values[e] > 0] + if not measured: + result.skipped.append(tuple(sorted(complex_))) + continue + per_copy = float(summarise(measured)) + for enzyme, copies in members.items(): + level = per_copy * copies + if enzyme in values and np.isfinite(values[enzyme]) and values[enzyme] > 0: + if not np.isclose(values[enzyme], level): + result.adjusted.append(enzyme) + values[enzyme] = level + elif fill_missing: + if enzyme not in values: + order.append(enzyme) + values[enzyme] = level + result.added.append(enzyme) + return order, np.array([values[p] for p in order], dtype=float), result + + +def level_complex_concs( + model: "EcModel", + *, + statistic: str = "median", + fill_missing: bool = True, +) -> LevelComplexResult: + """Give every complex in ``model.ec.concs`` one abundance. + + Run after ``fill_enz_concs`` and before ``constrain_enz_concs``. + Concentrations (mg/gDCW) are converted to molar amounts with + ``ec.mw`` for the levelling and back afterwards. Mutates + ``model.ec.concs``. + + Parameters + ---------- + model + EcModel with ``ec.concs`` filled. + statistic + ``"median"`` (default) or ``"mean"``. + fill_missing + Give unmeasured, non-shared subunits the complex's concentration. + """ + enzymes = list(model.ec.enzymes) + mw = np.asarray(model.ec.mw, dtype=float) + concs = np.asarray(model.ec.concs, dtype=float) + molar = concs / mw + measured = [i for i in range(len(enzymes)) if np.isfinite(molar[i])] + ids, levelled, result = level_complex_abundances( + model, [enzymes[i] for i in measured], molar[measured], + statistic=statistic, fill_missing=fill_missing, + ) + new = concs.copy() + index = {e: i for i, e in enumerate(enzymes)} + for enzyme, value in zip(ids, levelled): + i = index[enzyme] + new[i] = value * mw[i] + model.ec.concs = new + logger.info( + "level_complex_concs: %d subunits levelled, %d filled, %d shared left as measured", + len(result.adjusted), len(result.added), len(result.shared), + ) + return result diff --git a/tests/test_level_complex_concs.py b/tests/test_level_complex_concs.py new file mode 100644 index 0000000..c8bfe16 --- /dev/null +++ b/tests/test_level_complex_concs.py @@ -0,0 +1,171 @@ +"""Tests for complex_families, level_complex_concs and flexibilize_enz_concs(per_complex=True).""" +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter +from geckopy.ec_model.ec_data import EcData +from geckopy.limit_proteins import ( + complex_families, + flexibilize_enz_concs, + level_complex_abundances, + level_complex_concs, +) + +ENZYMES = ["A", "B", "C", "D", "F", "S", "E", "G", "H"] +MW = np.array([100.0, 200.0, 50.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]) + + +def _ec_only(concs, rows): + """An EcModel carrying only ec data: ``rows`` is {rxn: {enzyme: copies}}.""" + model = EcModel("toy") + matrix = np.zeros((len(rows), len(ENZYMES))) + for r, enzymes in enumerate(rows.values()): + for enzyme, copies in enzymes.items(): + matrix[r, ENZYMES.index(enzyme)] = copies + n = len(rows) + model.ec = EcData( + rxns=list(rows), kcat=np.ones(n), source=[""] * n, notes=[""] * n, eccodes=[""] * n, + genes=[f"g_{e}" for e in ENZYMES], enzymes=list(ENZYMES), mw=MW.copy(), + sequence=[""] * len(ENZYMES), concs=np.asarray(concs, dtype=float), + rxn_enz_mat=sparse.csr_matrix(matrix), + ) + return model + + +ROWS = { + "R1_EXP_1": {"A": 1, "B": 2}, + "R1_EXP_2": {"A": 1, "B": 2, "C": 1}, # isozyme variant: same complex + "R2": {"D": 1, "F": 1, "S": 1}, + "R3": {"E": 1, "G": 1, "S": 1}, # S is shared with R2's complex + "R4": {"H": 1}, # a single enzyme: not a complex +} + + +def test_isozyme_variants_are_one_complex_and_a_subunit_in_two_is_shared(): + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), ROWS)) + complexes = sorted(sorted(c) for c in families.complexes) + assert complexes == [["A", "B", "C"], ["D", "F", "S"], ["E", "G", "S"]] + assert families.shared == {"S"} + assert sorted(families.members("B")) == ["A", "B", "C"] + assert families.members("S") == ["S"] + assert families.members("H") == ["H"] + + +def test_an_enzyme_that_also_works_alone_is_shared(): + rows = dict(ROWS, R5={"A": 1}) + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert "A" in families.shared + + +def test_levelling_takes_the_median_per_copy_and_scales_by_stoichiometry(): + # molar per copy: A = 1, B = 4 / 2 = 2, C = 6; median 2. + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, result = level_complex_abundances(model, ["A", "B", "C"], [1.0, 4.0, 6.0]) + levelled = dict(zip(ids, values)) + assert levelled == pytest.approx({"A": 2.0, "B": 4.0, "C": 2.0}) + assert sorted(result.adjusted) == ["A", "C"] + + +def test_mean_is_available(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, _ = level_complex_abundances(model, ["A", "B", "C"], [1.0, 4.0, 6.0], statistic="mean") + assert dict(zip(ids, values))["A"] == pytest.approx(3.0) + + +def test_shared_subunits_stay_as_measured_and_missing_ones_are_filled(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, values, result = level_complex_abundances(model, ["D", "S"], [1.0, 9.0]) + levelled = dict(zip(ids, values)) + assert levelled["S"] == 9.0 # shared: untouched and not in the summary + assert levelled["D"] == levelled["F"] == 1.0 + assert "F" in result.added + + +def test_filling_can_be_switched_off(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + ids, _, _ = level_complex_abundances(model, ["D"], [1.0], fill_missing=False) + assert "F" not in ids + + +def test_level_complex_concs_works_in_mass_units(): + concs = np.full(len(ENZYMES), np.nan) + concs[ENZYMES.index("A")] = 100.0 # 1 mmol per copy at MW 100 + concs[ENZYMES.index("B")] = 1600.0 # 8 mmol, 4 per copy at MW 200 + concs[ENZYMES.index("C")] = 100.0 # 2 mmol per copy at MW 50 + model = _ec_only(concs, ROWS) + level_complex_concs(model) + new = dict(zip(ENZYMES, model.ec.concs)) + # median per copy is 2 mmol: A = 200 mg, B = 2 copies * 2 * 200 = 800 mg, C = 2 * 50 = 100 mg. + assert new["A"] == pytest.approx(200.0) + assert new["B"] == pytest.approx(800.0) + assert new["C"] == pytest.approx(100.0) + assert np.isnan(new["H"]) + + +def test_an_unknown_statistic_is_rejected(): + model = _ec_only(np.full(len(ENZYMES), np.nan), ROWS) + with pytest.raises(ValueError, match="statistic"): + level_complex_abundances(model, ["A"], [1.0], statistic="max") + + +# --------------------------------------------------------------------------- # +# flexibilize_enz_concs(per_complex=True) +# --------------------------------------------------------------------------- # + +def _adapter(tmp_path: Path, gr_exp: float) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text( + f'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = {gr_exp}\nbio_rxn = "biomass"\n' + ) + return ModelAdapter.from_folder(tmp_path) + + +def _complex_toy(adapter): + """A -> B catalysed by a two-subunit complex X1 + X2; growth equals the flux.""" + model = EcModel("toy", adapter=adapter) + a_e, a_c, b_c = (cobra.Metabolite(m, compartment="c") for m in ("A_e", "A_c", "B_c")) + pool = cobra.Metabolite("prot_pool", compartment="c") + x1, x2 = cobra.Metabolite("prot_X1", compartment="c"), cobra.Metabolite("prot_X2", compartment="c") + model.add_metabolites([a_e, a_c, b_c, pool, x1, x2]) + reactions = { + "EX_A": ({a_e: -1}, (-1000, 0)), + "TR_A": ({a_e: -1, a_c: 1}, (0, 1000)), + "R_AB": ({a_c: -1, b_c: 1, x1: -1, x2: -1}, (0, 1000)), + "biomass": ({b_c: -1}, (0, 1000)), + "prot_pool_exchange": ({pool: 1}, (0, 1000)), + "usage_prot_X1": ({pool: -1, x1: 1}, (0, 1.0)), + "usage_prot_X2": ({pool: -1, x2: 1}, (0, 3.0)), + } + for rid, (mets, bounds) in reactions.items(): + r = cobra.Reaction(rid) + r.add_metabolites(mets) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R_AB"], kcat=np.array([1.0]), source=[""], notes=[""], eccodes=[""], + genes=["g_X1", "g_X2"], enzymes=["X1", "X2"], mw=np.array([100.0, 100.0]), + sequence=["", ""], concs=np.array([1.0, 3.0]), + rxn_enz_mat=sparse.csr_matrix([[1.0, 1.0]]), + ) + return model + + +def test_without_per_complex_only_the_limiting_subunit_is_raised(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + result = flexibilize_enz_concs(model, verbose=False) + assert result.uniprot_ids == ["X1"] + assert model.reactions.usage_prot_X2.upper_bound == pytest.approx(3.0) + + +def test_per_complex_raises_every_subunit_by_the_same_factor(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + result = flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert sorted(result.uniprot_ids) == ["X1", "X2"] + factor_1 = model.reactions.usage_prot_X1.upper_bound / 1.0 + factor_2 = model.reactions.usage_prot_X2.upper_bound / 3.0 + assert factor_1 == pytest.approx(factor_2) + assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) From cfe7067d70ce1bbde3744cb6668980b7c46bcbd9 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 14 Oct 2026 23:46:14 +0200 Subject: [PATCH 2/5] fix: keep carriers out of complexes and relax per complex without lowering caps Two enzyme sets merge only when they share two or more subunits and half of the larger set, so electron carriers no longer chain complexes together; a subunit in two complexes belongs to the one with more reactions. per_complex ranks complexes by summed control coefficients, counts iterations per complex and never lowers a cap. --- .../limit_proteins/flexibilize_enz_concs.py | 52 +++++++---- .../limit_proteins/level_complex_concs.py | 89 +++++++++++++------ tests/test_level_complex_concs.py | 51 +++++++++++ 3 files changed, 150 insertions(+), 42 deletions(-) diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index 797a60b..cc76872 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -156,11 +156,15 @@ def flexibilize_enz_concs( verbose Whether per-iteration progress is logged at INFO. per_complex - When the most limiting enzyme is a non-shared subunit of a - complex (see ``complex_families``), raise every measured subunit - of that complex by the same factor, and keep one factor per - complex in the refinement pass, so that concentrations levelled - by ``level_complex_concs`` stay levelled. Off by default. + Treat every complex (see ``complex_families``) as one enzyme: + rank complexes by the summed control coefficients of their + subunits (the gain from raising them all together cannot exceed + that sum, and a levelled complex binds on all subunits at once, + so no single subunit shows a gain), raise every measured subunit + of the picked complex by the same factor, count iterations per + complex, and keep one factor per complex in the refinement pass, + so that concentrations levelled by ``level_complex_concs`` stay + levelled. Off by default. Returns ------- @@ -214,6 +218,20 @@ def flexibilize_enz_concs( touched = np.zeros(n_proteins, dtype=bool) position = {p: i for i, p in enumerate(proteins)} families: Optional[ComplexFamilies] = complex_families(model) if per_complex else None + # Indices into ``proteins`` of each unit that is relaxed together: a + # complex's measured, non-shared subunits, or a single enzyme. + units: list[list[int]] = [] + unit_of = np.zeros(n_proteins, dtype=int) + assigned: set[int] = set() + for i, protein in enumerate(proteins): + if i in assigned: + continue + members = [position[m] for m in families.members(protein) if m in position] if families else [i] + for j in members: + unit_of[j] = len(units) + assigned.update(members) + units.append(members) + unit_frequence = np.zeros(len(units), dtype=int) initial_sol = model.optimize() if ( @@ -264,9 +282,14 @@ def flexibilize_enz_concs( pred_growth = pool_pred_growth break - max_idx = int(np.argmax(control_coeffs)) + unit_coeffs = np.zeros(len(units)) + np.add.at(unit_coeffs, unit_of, control_coeffs) + unit = int(np.argmax(unit_coeffs)) + members = units[unit] + max_idx = max(members, key=lambda j: control_coeffs[j]) frequence[max_idx] += 1 - if frequence[max_idx] > iter_per_enzyme: + unit_frequence[unit] += 1 + if unit_frequence[unit] > iter_per_enzyme: logger.warning( "flexibilize_enz_concs: iter limit reached for protein %r; " "consider revising its kcat.", @@ -275,15 +298,14 @@ def flexibilize_enz_concs( flex_break = True break - increase = fold_change * frequence[max_idx] - members = families.members(proteins[max_idx]) if families else [proteins[max_idx]] - for member in members: - if member not in position: - continue - i = position[member] + increase = fold_change * unit_frequence[unit] + for i in members: touched[i] = True - usage_rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{member}") - usage_rxn.upper_bound = float(model.ec.concs[measured_idx[i]]) * (1.0 + increase) + usage_rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{proteins[i]}") + usage_rxn.upper_bound = max( + usage_rxn.upper_bound, + float(model.ec.concs[measured_idx[i]]) * (1.0 + increase), + ) new_sol = model.optimize() pred_growth = float(new_sol.objective_value or 0.0) diff --git a/src/geckopy/limit_proteins/level_complex_concs.py b/src/geckopy/limit_proteins/level_complex_concs.py index ae59f95..d060f04 100644 --- a/src/geckopy/limit_proteins/level_complex_concs.py +++ b/src/geckopy/limit_proteins/level_complex_concs.py @@ -9,11 +9,18 @@ stoichiometry. A complex is a set of enzymes that catalyses a reaction together, taken -from ``ec.rxn_enz_mat``. Enzyme sets that share at least half of the -smaller set's subunits are one complex (isozyme variants that swap a -subunit, the same complex on different reactions). A subunit that ends -up in more than one complex, or that also catalyses a reaction on its -own, is *shared*: it is left as measured and does not take part in the +from ``ec.rxn_enz_mat``. Two enzyme sets are one complex when they share +at least two subunits and at least half of the larger set (isozyme +variants that swap a subunit, the same complex on different reactions). +A single shared protein never merges two sets, so electron carriers and +shared subunits (cytochrome c, thioredoxin, the E3 subunit of the +2-oxo-acid dehydrogenases) do not chain unrelated complexes together. + +A subunit in exactly two complexes belongs to the one whose enzyme sets +catalyse more reactions (the core subunits COX1-3 of complex IV also +appear in one assembly reaction). A subunit in three or more complexes, +in two with equal weight, or that also catalyses a reaction on its own +is *shared*: it is left as measured and does not take part in the summary. """ from __future__ import annotations @@ -41,22 +48,26 @@ class ComplexFamilies: complexes One dict per complex, ``{enzyme: copies}`` over every subunit that appears in any of its enzyme sets. + owner + For every subunit that belongs to a complex, the index of that + complex in ``complexes``. shared - Enzymes that belong to more than one complex, or to a complex - and also catalyse a reaction alone. + Subunits of a complex that belong to none: carriers and subunits + shared between complexes. """ complexes: list[dict[str, float]] = field(default_factory=list) + owner: dict[str, int] = field(default_factory=dict) shared: set[str] = field(default_factory=set) + def owned(self, index: int) -> dict[str, float]: + """``{enzyme: copies}`` of the subunits complex ``index`` owns.""" + return {e: c for e, c in self.complexes[index].items() if self.owner.get(e) == index} + def members(self, enzyme: str) -> list[str]: - """The non-shared subunits of the complex ``enzyme`` belongs to, or ``[enzyme]``.""" - if enzyme in self.shared: - return [enzyme] - for complex_ in self.complexes: - if enzyme in complex_: - return [e for e in complex_ if e not in self.shared] - return [enzyme] + """The subunits of the complex ``enzyme`` belongs to, or ``[enzyme]``.""" + index = self.owner.get(enzyme) + return [enzyme] if index is None else list(self.owned(index)) @dataclass @@ -81,7 +92,9 @@ class LevelComplexResult: skipped: list[tuple[str, ...]] = field(default_factory=list) -def complex_families(model: "EcModel", min_overlap: float = 0.5) -> ComplexFamilies: +def complex_families( + model: "EcModel", min_overlap: float = 0.5, min_shared: int = 2, +) -> ComplexFamilies: """Group the enzymes of multi-subunit reactions into complexes. Parameters @@ -90,11 +103,14 @@ def complex_families(model: "EcModel", min_overlap: float = 0.5) -> ComplexFamil An EcModel with ``ec.rxn_enz_mat`` (copies per subunit). min_overlap Two enzyme sets are one complex when they share at least this - fraction of the smaller set's subunits. + fraction of the larger set's subunits... + min_shared + ...and at least this many subunits. """ matrix = model.ec.rxn_enz_mat.tocsr() enzymes = list(model.ec.enzymes) sets: dict[frozenset, dict[str, float]] = {} + reactions: dict[frozenset, int] = {} alone: set[str] = set() for row in range(matrix.shape[0]): start, end = matrix.indptr[row], matrix.indptr[row + 1] @@ -103,6 +119,7 @@ def complex_families(model: "EcModel", min_overlap: float = 0.5) -> ComplexFamil alone.update(present) elif present: sets.setdefault(frozenset(present), present) + reactions[frozenset(present)] = reactions.get(frozenset(present), 0) + 1 keys = list(sets) parent = list(range(len(keys))) @@ -121,22 +138,40 @@ def find(i: int) -> int: for x in range(len(indices)): for y in range(x + 1, len(indices)): i, j = indices[x], indices[y] - if len(keys[i] & keys[j]) >= min_overlap * min(len(keys[i]), len(keys[j])): + common = len(keys[i] & keys[j]) + if common >= min_shared and common >= min_overlap * max(len(keys[i]), len(keys[j])): parent[find(i)] = find(j) grouped: dict[int, dict[str, float]] = {} + weight: dict[int, dict[str, int]] = {} for i, key in enumerate(keys): + w = weight.setdefault(find(i), {}) + for enzyme in key: + w[enzyme] = w.get(enzyme, 0) + reactions[key] complex_ = grouped.setdefault(find(i), {}) for enzyme, copies in sets[key].items(): + if enzyme in complex_ and complex_[enzyme] != copies: + logger.warning( + "complex_families: %s has %g copies in one variant and %g in another; " + "the first is used", enzyme, complex_[enzyme], copies, + ) complex_.setdefault(enzyme, copies) - complexes = list(grouped.values()) - - count: dict[str, int] = {} - for complex_ in complexes: - for enzyme in complex_: - count[enzyme] = count.get(enzyme, 0) + 1 - shared = {e for e, n in count.items() if n > 1} | (alone & set(count)) - return ComplexFamilies(complexes=complexes, shared=shared) + roots = list(grouped) + complexes = [grouped[r] for r in roots] + + found: dict[str, list[tuple[int, int]]] = {} + for index, root in enumerate(roots): + for enzyme, w in weight[root].items(): + found.setdefault(enzyme, []).append((w, index)) + owner: dict[str, int] = {} + for enzyme, hits in found.items(): + if enzyme in alone or len(hits) > 2: + continue + hits.sort(reverse=True) + if len(hits) == 1 or hits[0][0] > hits[1][0]: + owner[enzyme] = hits[0][1] + shared = set(found) - set(owner) + return ComplexFamilies(complexes=complexes, owner=owner, shared=shared) def level_complex_abundances( @@ -178,8 +213,8 @@ def level_complex_abundances( order = list(uniprot_ids) result = LevelComplexResult(shared=sorted(families.shared & set(values))) - for complex_ in families.complexes: - members = {e: c for e, c in complex_.items() if e not in families.shared} + for index, complex_ in enumerate(families.complexes): + members = families.owned(index) measured = [values[e] / c for e, c in members.items() if np.isfinite(values.get(e, np.nan)) and values[e] > 0] if not measured: result.skipped.append(tuple(sorted(complex_))) diff --git a/tests/test_level_complex_concs.py b/tests/test_level_complex_concs.py index c8bfe16..416ce2d 100644 --- a/tests/test_level_complex_concs.py +++ b/tests/test_level_complex_concs.py @@ -169,3 +169,54 @@ def test_per_complex_raises_every_subunit_by_the_same_factor(tmp_path): factor_2 = model.reactions.usage_prot_X2.upper_bound / 3.0 assert factor_1 == pytest.approx(factor_2) assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) + + +def test_one_shared_protein_does_not_chain_complexes_together(): + rows = {"R1": {"A": 1, "B": 1}, "R2": {"A": 1, "C": 1}, "R3": {"A": 1, "D": 1, "E": 1, "F": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(sorted(c) for c in families.complexes) == [["A", "B"], ["A", "C"], ["A", "D", "E", "F"]] + assert families.shared == {"A"} + + +def test_an_electron_carrier_in_several_complexes_is_shared_not_a_bridge(): + # B, C, D form a complex that takes carrier S or carrier H (two variants); + # S also pairs with E, like cytochrome c with a dehydrogenase. + rows = {"R1_EXP_1": {"B": 1, "C": 1, "D": 1, "S": 1}, "R1_EXP_2": {"B": 1, "C": 1, "D": 1, "H": 1}, + "R2": {"E": 1, "S": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(sorted(c) for c in families.complexes) == [["B", "C", "D", "H", "S"], ["E", "S"]] + assert "S" in families.shared + assert "E" not in families.members("B") + + +def test_per_complex_relaxes_a_levelled_complex_whose_subunits_bind_together(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + model.reactions.usage_prot_X2.upper_bound = 1.0 + model.ec.concs = np.array([1.0, 1.0]) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound == pytest.approx(model.reactions.usage_prot_X2.upper_bound) + assert model.optimize().objective_value == pytest.approx(2.0, rel=5e-3) + + +def test_per_complex_never_lowers_a_cap_below_its_measurement(tmp_path): + model = _complex_toy(_adapter(tmp_path, gr_exp=2.0)) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound >= 1.0 + assert model.reactions.usage_prot_X2.upper_bound >= 3.0 + + +def test_a_subunit_in_two_complexes_belongs_to_the_one_with_more_reactions(): + # A, B, C, D is a complex on three reactions; A, B with E, F, G, H on one + # (an assembly step): A and B belong to the first, not shared. + rows = {"R1": {"A": 1, "B": 1, "C": 1, "D": 1}, "R2": {"A": 1, "B": 1, "C": 1, "D": 1}, + "R3": {"A": 1, "B": 1, "C": 1, "D": 1}, "R4": {"A": 1, "B": 1, "E": 1, "F": 1, "G": 1, "H": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert sorted(families.members("A")) == ["A", "B", "C", "D"] + assert sorted(families.members("E")) == ["E", "F", "G", "H"] + assert families.shared == set() + + +def test_a_subunit_in_two_equally_weighted_complexes_is_shared(): + rows = {"R1": {"A": 1, "B": 1, "C": 1, "D": 1}, "R2": {"A": 1, "B": 1, "E": 1, "F": 1, "G": 1, "H": 1}} + families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) + assert {"A", "B"} <= families.shared From 48e25d4ff2d4a782bfe0c5f39adf959debb256b0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 14 Oct 2026 23:56:15 +0200 Subject: [PATCH 3/5] fix: rank units by the cap duals, and probe with the pool lifted when no dual is positive --- .../limit_proteins/flexibilize_enz_concs.py | 108 ++++++++++++++---- tests/test_level_complex_concs.py | 40 +++++++ 2 files changed, 126 insertions(+), 22 deletions(-) diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index cc76872..4905278 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -157,14 +157,13 @@ def flexibilize_enz_concs( Whether per-iteration progress is logged at INFO. per_complex Treat every complex (see ``complex_families``) as one enzyme: - rank complexes by the summed control coefficients of their - subunits (the gain from raising them all together cannot exceed - that sum, and a levelled complex binds on all subunits at once, - so no single subunit shows a gain), raise every measured subunit - of the picked complex by the same factor, count iterations per - complex, and keep one factor per complex in the refinement pass, - so that concentrations levelled by ``level_complex_concs`` stay - levelled. Off by default. + rank enzymes and complexes by the growth gained per fold of + raising all their caps together (a levelled complex binds on all + subunits at once, so no single subunit shows a gain), raise every + measured subunit of the picked complex by the same factor, count + iterations per complex, and keep one factor per complex in the + refinement pass, so that concentrations levelled by + ``level_complex_concs`` stay levelled. Off by default. Returns ------- @@ -267,14 +266,17 @@ def flexibilize_enz_concs( ) flex_break = True break - _, control_coeffs = get_conc_control_coeffs( - model, - proteins=proteins, - fold_change=fold_change, - limit=_CONTROL_COEFF_LIMIT, - ) + if families is not None: + scores = _unit_scores(model, proteins, units, fold_change, _CONTROL_COEFF_LIMIT) + else: + _, scores = get_conc_control_coeffs( + model, + proteins=proteins, + fold_change=fold_change, + limit=_CONTROL_COEFF_LIMIT, + ) - if not (control_coeffs > 0).any(): + if not (scores > 0).any(): pool_old, pool_new, pool_pred_growth = _try_relax_pool( model, bio_rxn, exp_growth, pred_growth, verbose, ) @@ -282,18 +284,20 @@ def flexibilize_enz_concs( pred_growth = pool_pred_growth break - unit_coeffs = np.zeros(len(units)) - np.add.at(unit_coeffs, unit_of, control_coeffs) - unit = int(np.argmax(unit_coeffs)) + if families is not None: + unit = int(np.argmax(scores)) + else: + unit = int(unit_of[int(np.argmax(scores))]) members = units[unit] - max_idx = max(members, key=lambda j: control_coeffs[j]) + max_idx = members[0] frequence[max_idx] += 1 unit_frequence[unit] += 1 if unit_frequence[unit] > iter_per_enzyme: logger.warning( - "flexibilize_enz_concs: iter limit reached for protein %r; " + "flexibilize_enz_concs: iter limit reached for %s %r; " "consider revising its kcat.", - proteins[max_idx], + "complex" if len(members) > 1 else "protein", + [proteins[j] for j in members] if len(members) > 1 else proteins[max_idx], ) flex_break = True break @@ -455,6 +459,65 @@ def _build_result( return _finalize_and_sort(_maybe_add_pool(result, pool_old, pool_new)) +def _unit_scores( + model: "EcModel", + proteins: list[str], + units: list[list[int]], + fold_change: float, + limit: float, +) -> np.ndarray: + """Growth gained per fold of raising each unit's caps together. + + With duals, ``sum(reduced cost of usage_prot_i * its cap)`` over the + unit's subunits that run at more than ``limit`` of their cap. Unlike + the ``prot_`` shadow price, the usage reaction's reduced cost + leaves out the protein pool's price, which would otherwise be counted + once per subunit; and the weighted sum is the same for every optimal + dual when a levelled complex binds on all subunits at once. When no + dual is positive (no duals exposed, or the protein pool binds and the + solver puts the whole price on it), one solve per candidate unit with + all its caps raised by ``fold_change`` and the pool lifted. + """ + from .get_conc_control_coeffs import _solution_is_optimal, _solver_supports_duals + + scores = np.zeros(len(units)) + sol = model.optimize() + if not _solution_is_optimal(sol): + return scores + candidates = [] + for k, members in enumerate(units): + reactions = [model.reactions.get_by_id(f"{USAGE_PREFIX}{proteins[j]}") for j in members] + active = [r for r in reactions if r.upper_bound > 0 and float(sol.fluxes[r.id]) / r.upper_bound > limit] + if active: + candidates.append((k, reactions, active)) + if _solver_supports_duals(model): + reduced = sol.reduced_costs + for k, _, active in candidates: + scores[k] = sum(abs(float(reduced[r.id])) * r.upper_bound for r in active) + scores[scores < 1e-12] = 0.0 + if (scores > 0).any(): + return scores + # No dual singles a unit out (none exposed, or the protein pool binds as + # well and carries the whole price): probe each candidate with the pool + # lifted, so the probe finds the caps that limit once the pool is relaxed. + pool_rxn = model.reactions.get_by_id(POOL_EXCHANGE_ID) + with model: + pool_rxn.upper_bound = max(pool_rxn.upper_bound, 1000.0) + base = model.optimize() + if not _solution_is_optimal(base): + return scores + baseline = float(base.objective_value) + for k, reactions, _ in candidates: + with model: + for r in reactions: + r.upper_bound = r.upper_bound * fold_change + probe = model.optimize() + if _solution_is_optimal(probe): + gain = float(probe.objective_value) - baseline + scores[k] = gain / (fold_change - 1.0) if gain > _GROWTH_DELTA_THRESHOLD else 0.0 + return scores + + def _per_complex( proteins: list[str], old: np.ndarray, @@ -475,7 +538,8 @@ def _per_complex( ratios = [new[index[m]] / old[index[m]] for m in group if old[index[m]] > 0] factor = max(ratios) if ratios else 1.0 for m in group: - out[index[m]] = old[index[m]] * factor + if old[index[m]] > 0: + out[index[m]] = old[index[m]] * factor return out diff --git a/tests/test_level_complex_concs.py b/tests/test_level_complex_concs.py index 416ce2d..1a05254 100644 --- a/tests/test_level_complex_concs.py +++ b/tests/test_level_complex_concs.py @@ -220,3 +220,43 @@ def test_a_subunit_in_two_equally_weighted_complexes_is_shared(): rows = {"R1": {"A": 1, "B": 1, "C": 1, "D": 1}, "R2": {"A": 1, "B": 1, "E": 1, "F": 1, "G": 1, "H": 1}} families = complex_families(_ec_only(np.full(len(ENZYMES), np.nan), rows)) assert {"A", "B"} <= families.shared + + +def _series_toy(adapter): + """A -> B by complex (X1, X2, X3), B -> P by single enzyme E; the pool binds. + + At flux 1 every subunit uses 1 mg of a 1.25 mg cap (80%) and E uses + 1 mg of a 1 mg cap; the pool holds exactly 4 mg. Only E limits. + """ + model = EcModel("series", adapter=adapter) + mets = {n: cobra.Metabolite(n, compartment="c") for n in ("A", "B", "P", "pool", "prot_X1", "prot_X2", "prot_X3", "prot_E")} + model.add_metabolites(list(mets.values())) + spec = { + "EX_A": ({"A": 1}, (0, 1000)), + "R1": ({"A": -1, "B": 1, "prot_X1": -1, "prot_X2": -1, "prot_X3": -1}, (0, 1000)), + "R2": ({"B": -1, "P": 1, "prot_E": -1}, (0, 1000)), + "biomass": ({"P": -1}, (0, 1000)), + "prot_pool_exchange": ({"pool": 1}, (0, 4.0)), + **{f"usage_prot_{e}": ({"pool": -1, f"prot_{e}": 1}, (0, 1.25 if e != "E" else 1.0)) for e in ("X1", "X2", "X3", "E")}, + } + for rid, (st, bounds) in spec.items(): + r = cobra.Reaction(rid) + r.add_metabolites({mets[m]: c for m, c in st.items()}) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R1", "R2"], kcat=np.ones(2), source=[""] * 2, notes=[""] * 2, eccodes=[""] * 2, + genes=["g1", "g2", "g3", "gE"], enzymes=["X1", "X2", "X3", "E"], mw=np.full(4, 100.0), + sequence=[""] * 4, concs=np.array([1.25, 1.25, 1.25, 1.0]), + rxn_enz_mat=sparse.csr_matrix([[1.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]), + ) + return model + + +def test_per_complex_does_not_pick_a_large_complex_for_the_pool_price(tmp_path): + model = _series_toy(_adapter(tmp_path, gr_exp=1.1)) + flexibilize_enz_concs(model, verbose=False, per_complex=True) + for e in ("X1", "X2", "X3"): + assert model.reactions.get_by_id(f"usage_prot_{e}").upper_bound == pytest.approx(1.25) + assert model.optimize().objective_value >= 1.1 * 0.9975 - 1e-9 From 3990109e9a3d938374bf23ab06f4f4c2b67794be Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 15 Oct 2026 22:00:40 +0200 Subject: [PATCH 4/5] feat: relax_enz_concs, the least total relaxation of enzyme caps in one LP --- src/geckopy/__init__.py | 2 + src/geckopy/limit_proteins/__init__.py | 2 + src/geckopy/limit_proteins/relax_enz_concs.py | 196 ++++++++++++++++++ tests/test_relax_enz_concs.py | 109 ++++++++++ 4 files changed, 309 insertions(+) create mode 100644 src/geckopy/limit_proteins/relax_enz_concs.py create mode 100644 tests/test_relax_enz_concs.py diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index fd1b166..ad0a949 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -126,6 +126,7 @@ get_conc_control_coeffs, level_complex_abundances, level_complex_concs, + relax_enz_concs, relax_proteomics_greedy, ) @@ -225,6 +226,7 @@ "get_conc_control_coeffs", "level_complex_abundances", "level_complex_concs", + "relax_enz_concs", "relax_proteomics_greedy", # utilities "AddNewRxnsResult", diff --git a/src/geckopy/limit_proteins/__init__.py b/src/geckopy/limit_proteins/__init__.py index f7c5374..03b96dd 100644 --- a/src/geckopy/limit_proteins/__init__.py +++ b/src/geckopy/limit_proteins/__init__.py @@ -18,6 +18,7 @@ level_complex_abundances, level_complex_concs, ) +from .relax_enz_concs import relax_enz_concs from .relax_proteomics_greedy import ( GreedyRelaxResult, RelaxationStep, @@ -40,5 +41,6 @@ "get_conc_control_coeffs", "level_complex_abundances", "level_complex_concs", + "relax_enz_concs", "relax_proteomics_greedy", ] diff --git a/src/geckopy/limit_proteins/relax_enz_concs.py b/src/geckopy/limit_proteins/relax_enz_concs.py new file mode 100644 index 0000000..be3b8ba --- /dev/null +++ b/src/geckopy/limit_proteins/relax_enz_concs.py @@ -0,0 +1,196 @@ +"""Loosen proteomics constraints by the least total amount, in one LP. + +``flexibilize_enz_concs`` and ``relax_proteomics_greedy`` relax enzyme +caps one at a time: find the most limiting enzyme, raise it, re-solve, +repeat. Which enzymes end up raised then depends on how "most limiting" +is ranked and on the order in which they are found; on a large model +with noisy proteomics, small changes to that ranking can change the +result. ``relax_enz_concs`` asks the question once instead: what is the +smallest set of fold increases of the measured caps (and, if necessary, +of the protein pool) that lets the model reach the target growth rate? + +Every measured enzyme's cap becomes ``usage <= conc * (1 + x)`` and the +sum of the ``x`` is minimised under ``growth >= exp_growth``. The protein +pool may grow by a fold too, at ``pool_weight`` times the price of an +enzyme's fold, so that it grows only where raising enzymes does not +suffice. With ``per_complex=True`` the measured subunits of a complex +(see ``complex_families``) share one ``x``, so concentrations levelled +by ``level_complex_concs`` stay levelled. +""" +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING, Optional + +import numpy as np + +from ..ec_model.constants import POOL_EXCHANGE_ID, USAGE_PREFIX +from .flexibilize_enz_concs import ( + FlexEnzResult, + _finalize_and_sort, + _maybe_add_pool, +) +from .level_complex_concs import complex_families + +if TYPE_CHECKING: + from ..ec_model.ec_model import EcModel + +logger = logging.getLogger(__name__) + + +def relax_enz_concs( + model: "EcModel", + *, + exp_growth: Optional[float] = None, + bio_rxn: Optional[str] = None, + per_complex: bool = False, + weight: str = "relative", + pool_weight: float = 1000.0, + floor: float = 1e-9, +) -> FlexEnzResult: + """Relax enzyme caps by the least total amount that reaches ``exp_growth``. + + Parameters + ---------- + model + EcModel with ``model.ec.concs`` filled and the caps applied + (``fill_enz_concs`` and ``constrain_enz_concs``). The usage + upper bounds (and ``prot_pool_exchange``) are updated in place; + ``ec.concs`` itself is left unchanged, as in + ``flexibilize_enz_concs``. + exp_growth + Target growth rate. Defaults to ``params.gr_exp``. + bio_rxn + Biomass reaction id. Defaults to ``params.bio_rxn``. + per_complex + One fold per complex (its measured, owned subunits; see + ``complex_families``) instead of one per enzyme. + weight + ``"relative"`` minimises the summed fold increases, so an enzyme + measured at almost nothing is not relaxed freely because its + absolute cost is small; ``"absolute"`` minimises the mg/gDCW + added. + pool_weight + Price of one fold of the protein pool relative to one fold of an + enzyme (or, with ``weight="absolute"``, per mg). + floor + Concentration below which a measured value is treated as this + floor when scaling its fold, so a zero measurement can still be + relaxed (at a high relative price). + + Returns + ------- + FlexEnzResult + The relaxed enzymes, as from ``flexibilize_enz_concs`` + (``frequence`` is 1 for every entry), with a ``"prot_pool"`` + entry when the pool was raised. + + Raises + ------ + ValueError + If no enzyme is measured, the weighting is unknown, or no + relaxation reaches ``exp_growth`` (the target is then not limited + by protein). + """ + from optlang.symbolics import add + + from ..adapter import resolve_param + + if weight not in ("relative", "absolute"): + raise ValueError(f"weight must be 'relative' or 'absolute', not {weight!r}") + if exp_growth is None: + exp_growth = resolve_param(model, None, "gr_exp", purpose="relax_enz_concs needs `exp_growth`") + if exp_growth is None: + raise ValueError("exp_growth not provided and params.gr_exp is unset.") + exp_growth = float(exp_growth) + bio = model.reactions.get_by_id( + resolve_param(model, bio_rxn, "bio_rxn", purpose="relax_enz_concs needs the biomass reaction id") + ) + + concs = np.asarray(model.ec.concs, dtype=float) + measured = [i for i in range(len(concs)) if np.isfinite(concs[i])] + if not measured: + raise ValueError("No enzyme concentrations are measured. Run fill_enz_concs and constrain_enz_concs first.") + enzymes = list(model.ec.enzymes) + + old_caps = {i: model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound for i in measured} + # folds count from the current cap, so a cap raised before is never lowered + base = {i: min(old_caps[i], 1000.0) for i in measured} + + families = complex_families(model) if per_complex else None + position = {e: i for i, e in enumerate(enzymes)} + measured_set = set(measured) + units: list[list[int]] = [] + assigned: set[int] = set() + for i in measured: + if i in assigned: + continue + if families is None or not base[i] > 0: + unit = [i] + else: # a subunit measured at zero cannot share a fold with the others + unit = [position[m] for m in families.members(enzymes[i]) + if m in position and position[m] in measured_set + and position[m] not in assigned and base[position[m]] > 0] + assigned.update(unit) + units.append(unit) + + pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) + pool_cap = pool.upper_bound + old_bio_bounds = bio.bounds + previous_objective = model.objective.expression + previous_direction = model.objective.direction + + folds, added = [], [] + for k, unit in enumerate(units): + x = model.problem.Variable(f"relax_unit_{k}", lb=0.0, ub=None) + added.append(x) + for i in unit: + rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}") + rxn.upper_bound = 1000.0 + added.append(model.problem.Constraint( + rxn.flux_expression - max(base[i], floor) * x, ub=base[i], name=f"relax_cap_{enzymes[i]}", + )) + folds.append((x, unit)) + pool_x = model.problem.Variable("relax_pool", lb=0.0, ub=None) + pool.upper_bound = 1000.0 + added += [pool_x, model.problem.Constraint(pool.flux_expression - pool_cap * pool_x, ub=pool_cap, name="relax_cap_pool")] + bio.bounds = (exp_growth, max(bio.upper_bound, exp_growth)) + model.add_cons_vars(added) + + try: + terms = [(1.0 if weight == "relative" else sum(max(base[i], floor) for i in unit)) * x for x, unit in folds] + terms.append(pool_weight * (1.0 if weight == "relative" else pool_cap) * pool_x) + model.objective = model.problem.Objective(add(terms), direction="min") + solution = model.optimize() + feasible = solution.status == "optimal" + values = model.solver.primal_values if feasible else {} + finally: + model.remove_cons_vars(added) + model.objective = model.problem.Objective(previous_objective, direction=previous_direction) + bio.bounds = (old_bio_bounds[0], max(old_bio_bounds[1], exp_growth)) + + if not feasible: + for i, cap in old_caps.items(): + model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound = cap + pool.upper_bound = pool_cap + raise ValueError( + f"relax_enz_concs: no relaxation of the measured enzymes or the protein pool reaches " + f"growth {exp_growth:g} ({solution.status}); the target is not limited by protein." + ) + + ids, old, new = [], [], [] + for x, unit in folds: + fold = float(values.get(x.name, 0.0)) + for i in unit: + cap = base[i] + max(base[i], floor) * fold + model.reactions.get_by_id(f"{USAGE_PREFIX}{enzymes[i]}").upper_bound = cap if fold > 1e-9 else old_caps[i] + if fold > 1e-9: + ids.append(enzymes[i]); old.append(base[i]); new.append(cap) + pool_fold = float(values.get(pool_x.name, 0.0)) + pool_new = pool_cap * (1 + pool_fold) if pool_fold > 1e-9 else None + pool.upper_bound = pool_new if pool_new is not None else pool_cap + logger.info("relax_enz_concs: %d enzymes relaxed%s", len(ids), + f", protein pool {pool_cap:g} -> {pool_new:g}" if pool_new is not None else "") + result = FlexEnzResult(uniprot_ids=ids, old_concs=np.array(old, dtype=float), + flex_concs=np.array(new, dtype=float), frequence=np.ones(len(ids), dtype=int)) + return _finalize_and_sort(_maybe_add_pool(result, pool_cap if pool_new is not None else None, pool_new)) diff --git a/tests/test_relax_enz_concs.py b/tests/test_relax_enz_concs.py new file mode 100644 index 0000000..2b09e8f --- /dev/null +++ b/tests/test_relax_enz_concs.py @@ -0,0 +1,109 @@ +"""Tests for relax_enz_concs: the least total relaxation, in one LP.""" +from pathlib import Path + +import cobra +import numpy as np +import pytest +from scipy import sparse + +from geckopy import EcModel, ModelAdapter +from geckopy.ec_model.ec_data import EcData +from geckopy.limit_proteins import relax_enz_concs + + +def _adapter(tmp_path: Path, gr_exp: float) -> ModelAdapter: + (tmp_path / "model_adapter.toml").write_text( + f'conv_gem = "dummy.xml"\norg_name = "test"\ngr_exp = {gr_exp}\nbio_rxn = "biomass"\n' + ) + return ModelAdapter.from_folder(tmp_path) + + +def _two_routes(adapter, caps, pool=1000.0): + """A -> B by enzyme E1 (2 mg per unit flux) or E2 (1 mg per unit); growth = flux to B. + + ``caps`` gives the measured concentration of each enzyme (mg). + """ + model = EcModel("routes", adapter=adapter) + mets = {n: cobra.Metabolite(n, compartment="c") for n in ("A", "B", "pool", "prot_E1", "prot_E2")} + model.add_metabolites(list(mets.values())) + spec = { + "EX_A": ({"A": 1}, (0, 10)), + "R1": ({"A": -1, "B": 1, "prot_E1": -2}, (0, 1000)), + "R2": ({"A": -1, "B": 1, "prot_E2": -1}, (0, 1000)), + "biomass": ({"B": -1}, (0, 1000)), + "prot_pool_exchange": ({"pool": 1}, (0, pool)), + "usage_prot_E1": ({"pool": -1, "prot_E1": 1}, (0, caps["E1"])), + "usage_prot_E2": ({"pool": -1, "prot_E2": 1}, (0, caps["E2"])), + } + for rid, (st, bounds) in spec.items(): + r = cobra.Reaction(rid) + r.add_metabolites({mets[m]: c for m, c in st.items()}) + r.bounds = bounds + model.add_reactions([r]) + model.objective = "biomass" + model.ec = EcData( + rxns=["R1", "R2"], kcat=np.ones(2), source=[""] * 2, notes=[""] * 2, eccodes=[""] * 2, + genes=["g1", "g2"], enzymes=["E1", "E2"], mw=np.full(2, 100.0), sequence=["", ""], + concs=np.array([caps["E1"], caps["E2"]], dtype=float), rxn_enz_mat=sparse.csr_matrix(np.eye(2)), + ) + return model + + +def test_the_least_summed_fold_increase_reaches_the_target(tmp_path): + # 8 units: x on E1 (2 mg each) and 8 - x on E2 (1 mg each), both capped at 2 mg. + # Folds (x - 1)+ and (3 - x/2)+ sum to 2 + x/2 on [1, 6], 3 - x/2 below: least is 2.5 at x = 1, + # E1 kept at 2 mg and E2 raised to 7 mg. + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + result = relax_enz_concs(model) + assert model.optimize().objective_value == pytest.approx(8.0) + folds = model.reactions.usage_prot_E1.upper_bound / 2.0 - 1 + model.reactions.usage_prot_E2.upper_bound / 2.0 - 1 + assert folds == pytest.approx(2.5) + assert model.reactions.usage_prot_E2.upper_bound == pytest.approx(7.0) + assert "prot_pool" not in result.uniprot_ids + + +def test_enough_capacity_relaxes_nothing(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=3.0), {"E1": 2.0, "E2": 2.0}) + result = relax_enz_concs(model) + assert result.uniprot_ids == [] + assert model.reactions.usage_prot_E1.upper_bound == 2.0 + + +def test_the_pool_grows_only_when_enzymes_cannot_help(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 100.0, "E2": 100.0}, pool=4.0) + result = relax_enz_concs(model) + assert result.uniprot_ids == ["prot_pool"] + assert model.reactions.prot_pool_exchange.upper_bound == pytest.approx(8.0) # all on E2 + + +def test_a_cap_raised_before_is_not_lowered(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + model.reactions.usage_prot_E2.upper_bound = 6.0 # raised earlier + relax_enz_concs(model) + assert model.reactions.usage_prot_E2.upper_bound >= 6.0 + assert model.optimize().objective_value == pytest.approx(8.0) + + +def test_an_unreachable_target_is_an_error_and_leaves_the_caps(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=20.0), {"E1": 2.0, "E2": 2.0}) # supply caps growth at 10 + with pytest.raises(ValueError, match="not limited by protein"): + relax_enz_concs(model) + assert model.reactions.usage_prot_E1.upper_bound == 2.0 + assert model.reactions.prot_pool_exchange.upper_bound == 1000.0 + + +def test_per_complex_raises_the_subunits_of_a_complex_together(tmp_path): + from test_level_complex_concs import _adapter as adapter2, _complex_toy + + model = _complex_toy(adapter2(tmp_path, gr_exp=2.0)) + model.ec.concs = np.array([1.0, 1.0]) + model.reactions.usage_prot_X2.upper_bound = 1.0 + relax_enz_concs(model, per_complex=True) + assert model.reactions.usage_prot_X1.upper_bound == pytest.approx(2.0) + assert model.reactions.usage_prot_X2.upper_bound == pytest.approx(2.0) + + +def test_an_unknown_weighting_is_rejected(tmp_path): + model = _two_routes(_adapter(tmp_path, gr_exp=8.0), {"E1": 2.0, "E2": 2.0}) + with pytest.raises(ValueError, match="weight"): + relax_enz_concs(model, weight="max") From 62b64f8323ddfdceb938144205704ee2f0e93437 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 15 Oct 2026 22:01:10 +0200 Subject: [PATCH 5/5] docs: describe relax_enz_concs and the complex options --- CHANGELOG.md | 12 +++++++ docs/relaxation_methods.md | 69 +++++++++++++++++++++++++++++++++++++- 2 files changed, 80 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e57433..065508a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,18 @@ ## Unreleased +* New: `level_complex_concs` and `level_complex_abundances` give every enzyme + complex one abundance, the median per copy over its measured subunits, so a + single badly quantified subunit no longer sets a complex's capacity. + `complex_families` groups isozyme variants into complexes and leaves shared + subunits and electron carriers as measured. See `docs/relaxation_methods.md`. +* New: `relax_enz_concs` relaxes enzyme caps (and, if need be, the protein + pool) by the least total fold increase that reaches the target growth, in one + LP, so the result does not depend on the order in which limiting enzymes are + found. +* `flexibilize_enz_concs` gained `per_complex=True`: complexes are ranked and + raised as one unit, keeping levelled subunits levelled. Off by default. + * Fix: `write_dlkcat_input` and the OpenKineticsPredictor input builder raised `TypeError: cannot use 'list' as a set element` on any model that had been saved and reloaded. cobra returns an annotation value as a list diff --git a/docs/relaxation_methods.md b/docs/relaxation_methods.md index 29eb2d8..3b823d5 100644 --- a/docs/relaxation_methods.md +++ b/docs/relaxation_methods.md @@ -1,4 +1,4 @@ -# `flexibilize_enz_concs` vs `relax_proteomics_greedy` +# Relaxing proteomics constraints: `flexibilize_enz_concs`, `relax_proteomics_greedy` and `relax_enz_concs` Both functions answer the same question โ€” *the measured proteomics constraints make the ecModel grow too slowly (or not at all); which enzyme bounds do we @@ -6,6 +6,8 @@ loosen, and by how much, to restore feasible growth?* โ€” but they are **not** duplicates. They use different selection signals, different relaxation step sizes, and produce results of very different stringency. This note records the differences and the recommendation (keep both; do not merge). +`relax_enz_concs`, a third method that solves the same question as one LP, and +the per-complex options for all of them are described at the end. ## Side by side @@ -50,3 +52,68 @@ follow-ups would make them more interchangeable: the over-relaxation point in the modeling-intent review, ยง3.3). 2. Align the failure convention โ€” currently one returns a `converged=False` result and the other raises; pick one so callers don't special-case each. + +## `relax_enz_concs`: the least total relaxation, in one LP + +Both methods above relax one enzyme at a time: rank, raise the top one, +re-solve, repeat. Which enzymes end up raised then depends on the ranking +signal and on the order in which enzymes become limiting. On a large model +with noisy proteomics, a small change to the ranking (for example, control +coefficient per mg against gain per fold) can change which enzymes are +raised and, downstream, predictions that depend on the remaining caps. + +`relax_enz_concs` asks the question once. Every measured cap becomes +`usage <= conc * (1 + x)`, the protein pool may grow by a fold too, and one LP +minimises the sum of the `x` subject to `growth >= exp_growth`: + +- `weight="relative"` (default) minimises summed fold increases, so an enzyme + measured at almost nothing is not relaxed freely because its absolute cost is + small; `weight="absolute"` minimises the mg added. +- The pool's fold costs `pool_weight` (default 1000) times an enzyme's, so the + pool grows only where raising enzymes cannot reach the target, the role of + the pool fallback in `flexibilize_enz_concs`. +- Folds count from the current caps, so a cap raised earlier is never lowered. +- If no relaxation reaches the target, it raises `ValueError`: growth is then + not limited by protein. + +It returns a `FlexEnzResult` like `flexibilize_enz_concs` (`frequence` is 1 for +every entry). Use it when the outcome should not depend on the order of +relaxation, and to compare with the iterative methods: if they disagree on +which enzymes need raising, conclusions that hinge on those enzymes are +fragile. + +## Enzyme complexes + +Absolute proteomics recovers the subunits of one complex unevenly: subunits +present 1:1 can come out several-fold apart while their replicates agree within +a few percent. Capped one by one, a complex can then carry no more flux than +its least abundant measurement allows, and relaxation methods that raise one +subunit at a time leave its subunits at unrelated levels, or raise a different +subunit in each condition. + +- `level_complex_concs` (after `fill_enz_concs`, before `constrain_enz_concs`) + gives every complex one abundance: the median per copy over its measured + subunits, scaled back to each subunit by its copies in `ec.rxn_enz_mat`; + unmeasured subunits are filled. `level_complex_abundances` does the same on + molar data before it is written into the model. +- `complex_families` defines the complexes. Two enzyme sets are one complex when + they share at least two subunits and half the larger set (isozyme variants, + the same complex on several reactions). A protein in three or more complexes, + in two with equal weight, or that also catalyses a reaction alone is *shared* + and keeps its own measurement: electron carriers (cytochrome c), thioredoxin, + or the E3 subunit of the 2-oxo-acid dehydrogenases, which would otherwise + chain unrelated complexes together. A subunit in two complexes belongs to the + one with more reactions (core subunits that also appear in an assembly step). +- `flexibilize_enz_concs(per_complex=True)` and `relax_enz_concs(per_complex=True)` + raise the measured subunits of a complex by one factor, so levelled + concentrations stay levelled. In `flexibilize_enz_concs` complexes are ranked + by the growth gained per fold of raising all their caps (the usage reactions' + reduced costs times the caps), which leaves out the protein pool's price and + does not depend on which of several simultaneously binding subunits the solver + puts the dual on. + +Levelling relies on the stoichiometry in `ec.rxn_enz_mat`; where +`apply_complex_data` found no exact ComplexPortal match, every subunit counts +as one copy. Per-complex relaxation without levelling raises every subunit by +the fold its worst-measured subunit needs, which inflates the well-measured +ones; use the two together.