From 7d5dab30cb29faaa0fe5734eb0380b5ce623860e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 01:10:28 +0200 Subject: [PATCH 1/4] feat: put the ribosome in the way of protein synthesis An ecModel charges protein cost to the enzymes carrying metabolic flux but not to the machinery that makes them, which understates the cost of growth. The protein pseudoreaction now makes a pool of amino acids, and a translation reaction turns that pool into protein while drawing on the ribosomal subunits, so every gram of protein pays for the ribosomes that made it. The subunits are selected on their average abundance across all conditions: 48 of the 253 candidates, which is the same 48 UniProt entries ribosome.m added to each of the five published condition models. Masses in ribosome.txt are written with thousands separators, so 43,758 is one number and not two. --- src/overflow/ribosome.py | 294 ++++++++++++++++++++++++++++++ tests/test_ribosome.py | 375 +++++++++++++++++++++++++++++++++++++++ tests/tiny.py | 38 ++++ 3 files changed, 707 insertions(+) create mode 100644 src/overflow/ribosome.py create mode 100644 tests/test_ribosome.py diff --git a/src/overflow/ribosome.py b/src/overflow/ribosome.py new file mode 100644 index 0000000..7e5e5dd --- /dev/null +++ b/src/overflow/ribosome.py @@ -0,0 +1,294 @@ +"""Putting the ribosome in the way of protein synthesis. + +An ecModel charges protein cost to the enzymes that carry metabolic +flux, but not to the machinery that makes those enzymes. Translation is +the largest single protein expense a growing cell has, so leaving it out +understates the cost of growth. + +The protein pseudoreaction is split in two: it now produces a pool of +amino acids, and a new ``translation`` reaction turns that pool into +protein while drawing on the ribosomal subunits. Every gram of protein +the cell makes therefore pays for the ribosomes that made it. +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path +from typing import TYPE_CHECKING, Iterable, Optional, Sequence + +import cobra +import numpy as np +import pandas as pd +from geckopy import set_kcat_for_reactions +from geckopy.utilities import NewEnzyme, add_new_rxns_to_ec + +from overflow.config import PROTEIN_RXN, RIBOSOME_DATA + +if TYPE_CHECKING: # pragma: no cover - typing only + from geckopy import EcModel + +#: Elongation rate of the yeast ribosome [amino acids/s], doi:10.1042/bj1680409. +AA_PER_SECOND = 10.5 + +#: Subunits averaging less than this across all conditions [mmol/gDW] are +#: treated as accessory rather than part of the core ribosome. +MIN_MEAN_ABUNDANCE = 1e-5 + +TRANSLATION_RXN = "translation" +AMINO_ACID_MET = "s_aminoAcids" +AMINO_ACID_NAME = "amino acids for protein" +SUBSYSTEM = "sce03010 Ribosome" + + +@dataclass +class RibosomeSubunits: + """The subunits taken to make up the core ribosome.""" + + uniprot_ids: list[str] + genes: dict[str, str] + masses: dict[str, float] + mean_abundance: dict[str, float] + n_candidates: int = 0 + + def __len__(self) -> int: + return len(self.uniprot_ids) + + +@dataclass +class RibosomeResult: + """What adding the ribosome changed.""" + + subunits: list[str] = field(default_factory=list) + enzymes_added: list[str] = field(default_factory=list) + reactions_added: list[str] = field(default_factory=list) + amino_acid_demand: float = 0.0 + kcat: float = 0.0 + + +def read_ribosome(path: Path | str = RIBOSOME_DATA) -> pd.DataFrame: + """Read the ribosomal subunit table. + + Masses are written with thousands separators, so ``43,758`` is one + number rather than two; they are daltons, matching ``ec.mw``. + """ + table = pd.read_csv(path, sep="\t") + table = table.rename( + columns={ + "Entry": "uniprot", + "Gene names (primary )": "gene_name", + "Gene names (ordered locus )": "gene", + "Mass": "mass", + "Sequence": "sequence", + } + ) + table["mass"] = pd.to_numeric( + table["mass"].astype(str).str.replace(",", "", regex=False), errors="coerce" + ) + return table + + +def core_subunits( + ribosome_table: pd.DataFrame, + proteomics_table: pd.DataFrame, + threshold: float = MIN_MEAN_ABUNDANCE, +) -> RibosomeSubunits: + """Select the subunits abundant enough to count as the core ribosome. + + Abundance is averaged over every replicate of every condition: a + subunit belongs to the core or it does not, and that should not be + decided one condition at a time. + """ + value_columns = [c for c in proteomics_table.columns if c.endswith("_abs")] + measured = proteomics_table.set_index("Protein.IDs")[value_columns] + + ids, genes, masses, means = [], {}, {}, {} + for _, row in ribosome_table.iterrows(): + uniprot = str(row["uniprot"]) + if uniprot not in measured.index: + continue + values = pd.to_numeric(measured.loc[uniprot], errors="coerce").to_numpy(float) + if np.isnan(values).all(): + continue + mean = float(np.nanmean(values)) + if not np.isfinite(mean) or mean < threshold: + continue + ids.append(uniprot) + genes[uniprot] = str(row["gene"]) + masses[uniprot] = float(row["mass"]) + means[uniprot] = mean + + return RibosomeSubunits( + uniprot_ids=ids, + genes=genes, + masses=masses, + mean_abundance=means, + n_candidates=len(ribosome_table), + ) + + +def amino_acid_demand(model: cobra.Model, protein_rxn: str = PROTEIN_RXN) -> float: + """Amino acids consumed per unit of protein [mmol/mmol]. + + Read from the protein pseudoreaction's substrates, so it follows any + rescaling of the biomass composition. + """ + reaction = model.reactions.get_by_id(protein_rxn) + return float(-sum(c for c in reaction.metabolites.values() if c < 0)) + + +def translation_kcat( + model: cobra.Model, + aa_per_second: float = AA_PER_SECOND, + protein_rxn: str = PROTEIN_RXN, +) -> float: + """Turnover of the ribosome in units of protein [1/s]. + + A ribosome adds ``aa_per_second`` amino acids per second, and one + unit of protein takes ``amino_acid_demand`` of them, so in units of + the protein pseudoreaction it turns over that many times slower. + """ + return aa_per_second / amino_acid_demand(model, protein_rxn) + + +def add_ribosome( + model: "EcModel", + subunits: RibosomeSubunits, + aa_per_second: float = AA_PER_SECOND, + protein_rxn: str = PROTEIN_RXN, +) -> RibosomeResult: + """Route protein synthesis through a ribosome-catalysed reaction. + + The protein pseudoreaction is left making everything it made before + except protein itself, which now comes out of ``translation``. + """ + if AMINO_ACID_MET in {m.id for m in model.metabolites} or any( + r.id == TRANSLATION_RXN or r.id.startswith(TRANSLATION_RXN + "_") + for r in model.reactions + ): + raise ValueError("the ribosome is already in this model") + if not subunits.uniprot_ids: + raise ValueError("no ribosomal subunits to add") + + demand = amino_acid_demand(model, protein_rxn) + kcat = aa_per_second / demand + + reaction = model.reactions.get_by_id(protein_rxn) + protein = _protein_metabolite(reaction) + produced = reaction.metabolites[protein] + + pool = cobra.Metabolite( + AMINO_ACID_MET, name=AMINO_ACID_NAME, compartment=protein.compartment + ) + model.add_metabolites([pool]) + reaction.add_metabolites({protein: -produced, pool: produced}) + + translation = cobra.Reaction( + TRANSLATION_RXN, name=TRANSLATION_RXN, lower_bound=0.0, upper_bound=1000.0 + ) + translation.add_metabolites({pool: -produced, protein: produced}) + translation.subsystem = SUBSYSTEM + translation.gene_reaction_rule = " and ".join( + subunits.genes[u] for u in subunits.uniprot_ids + ) + + added = add_new_rxns_to_ec( + model, + [translation], + [ + NewEnzyme(enzyme=u, gene=subunits.genes[u], mw=subunits.masses[u]) + for u in subunits.uniprot_ids + ], + ) + set_kcat_for_reactions(model, [TRANSLATION_RXN], kcat, apply=True) + + return RibosomeResult( + subunits=list(subunits.uniprot_ids), + enzymes_added=list(added.enz_added), + reactions_added=list(added.rxns_added), + amino_acid_demand=demand, + kcat=kcat, + ) + + +def _protein_metabolite(reaction: cobra.Reaction) -> cobra.Metabolite: + """The protein the pseudoreaction produces.""" + produced = [m for m, c in reaction.metabolites.items() if c > 0 and m.name == "protein"] + if len(produced) != 1: + raise ValueError( + f"{reaction.id} does not produce exactly one metabolite named " + f"'protein' (found {len(produced)})" + ) + return produced[0] + + +def subunit_coefficient( + mass: float, demand: float, aa_per_second: float = AA_PER_SECOND +) -> float: + """Ribosomal protein needed per unit of protein made [mg/mmol].""" + return mass * demand / (aa_per_second * 3600.0) + + +@dataclass +class SubunitConstraints: + """How the measured subunit abundances were applied.""" + + applied: dict[str, float] = field(default_factory=dict) + adjusted: list[str] = field(default_factory=list) + measured: list[float] = field(default_factory=list) + required: list[float] = field(default_factory=list) + + def table(self) -> pd.DataFrame: + return pd.DataFrame( + { + "protein": self.adjusted, + "measured_mg_gDW": self.measured, + "adjusted_mg_gDW": self.required, + "fold_change": [ + r / m if m else np.inf + for m, r in zip(self.measured, self.required) + ], + } + ) + + +def constrain_subunits( + model: "EcModel", + abundances: dict[str, float], + subunits: Optional[Iterable[str]] = None, + tolerance: float = 1.01, +) -> SubunitConstraints: + """Cap the ribosomal subunits at their measured abundances [mg/gDW]. + + What the model needs is read first, with the subunits unconstrained. + A subunit measured at less than that is raised to what the model + needs rather than being allowed to stop translation, on the grounds + that the cell was demonstrably growing. + + A subunit with no measurement keeps drawing on the protein pool. + """ + if subunits is None: + subunits = list(abundances) + result = SubunitConstraints() + + solution = model.optimize() + for unit in subunits: + reaction = model.reactions.get_by_id(f"usage_prot_{unit}") + measured = abundances.get(unit) + if measured is None: + continue + required = float(solution.fluxes[reaction.id]) if solution.status == "optimal" else 0.0 + if required > measured: + bound = tolerance * required + result.adjusted.append(unit) + result.measured.append(measured) + result.required.append(bound) + else: + bound = measured + reaction.upper_bound = bound + result.applied[unit] = bound + + concentrations = np.asarray(model.ec.concs, dtype=float).copy() + for unit, bound in result.applied.items(): + concentrations[model.ec.enzymes.index(unit)] = bound + model.ec.concs = concentrations + return result diff --git a/tests/test_ribosome.py b/tests/test_ribosome.py new file mode 100644 index 0000000..59a698c --- /dev/null +++ b/tests/test_ribosome.py @@ -0,0 +1,375 @@ +"""Putting the ribosome in the way of protein synthesis. + +Every number here is worked out by hand from the tiny model: one unit of +protein costs two amino acids, so a ribosome elongating at 10.5 aa/s +turns over 5.25 times a second in units of protein, and a 30 kDa subunit +costs 30000 x 2 / (10.5 x 3600) = 1.5873 mg per mmol of protein made. +""" +import numpy as np +import pandas as pd +import pytest + +from overflow.ribosome import ( + AA_PER_SECOND, + AMINO_ACID_MET, + MIN_MEAN_ABUNDANCE, + RibosomeSubunits, + add_ribosome, + amino_acid_demand, + core_subunits, + read_ribosome, + subunit_coefficient, + translation_kcat, +) +from tiny import TINY_AA_PER_PROTEIN, tiny_translation_model + +SUBUNIT_MASS = 30_000.0 +EXPECTED_COEFFICIENT = SUBUNIT_MASS * TINY_AA_PER_PROTEIN / (AA_PER_SECOND * 3600.0) + + +@pytest.fixture +def model(): + return tiny_translation_model(pool=200.0, supply=10.0) + + +@pytest.fixture +def subunits(): + ids = ["R1", "R2"] + return RibosomeSubunits( + uniprot_ids=ids, + genes={u: f"gene_{u}" for u in ids}, + masses={u: SUBUNIT_MASS for u in ids}, + mean_abundance={u: 1e-4 for u in ids}, + n_candidates=2, + ) + + +def _translation(model): + return next( + r for r in model.reactions + if r.id == "translation" or r.id.startswith("translation_") + ) + + +# --- the arithmetic --------------------------------------------------- + +def test_amino_acid_demand_is_read_from_the_pseudoreaction(model): + assert amino_acid_demand(model, "r_protein") == pytest.approx(TINY_AA_PER_PROTEIN) + + +def test_the_ribosome_turns_over_once_per_protein_it_finishes(model): + """10.5 amino acids a second, two per protein, so 5.25 proteins a second.""" + assert translation_kcat(model, protein_rxn="r_protein") == pytest.approx(5.25) + + +def test_a_protein_that_takes_more_amino_acids_costs_more_ribosome(): + cheap = subunit_coefficient(SUBUNIT_MASS, demand=2.0) + dear = subunit_coefficient(SUBUNIT_MASS, demand=4.0) + assert dear == pytest.approx(2 * cheap) + + +# --- the rewiring ----------------------------------------------------- + +def test_protein_now_comes_out_of_translation(model, subunits): + add_ribosome(model, subunits, protein_rxn="r_protein") + + pseudoreaction = model.reactions.get_by_id("r_protein") + protein = model.metabolites.get_by_id("protein") + assert protein not in pseudoreaction.metabolites + assert model.metabolites.get_by_id(AMINO_ACID_MET) in pseudoreaction.metabolites + + translation = _translation(model) + assert translation.metabolites[protein] == pytest.approx(1.0) + assert translation.metabolites[model.metabolites.get_by_id(AMINO_ACID_MET)] == pytest.approx(-1.0) + + +def test_each_subunit_becomes_an_enzyme_drawing_on_the_pool(model, subunits): + add_ribosome(model, subunits, protein_rxn="r_protein") + + for unit in subunits.uniprot_ids: + assert unit in model.ec.enzymes + index = model.ec.enzymes.index(unit) + assert model.ec.mw[index] == pytest.approx(SUBUNIT_MASS) + usage = model.reactions.get_by_id(f"usage_prot_{unit}") + assert any(m.id == "prot_pool" for m in usage.metabolites) + + +def test_the_subunit_cost_matches_the_elongation_rate(model, subunits): + """The round trip that matters: a kcat written through geckopy has to + come back out as the protein cost the elongation rate implies. This + is where a confusion between daltons and kDa would show up, as a + thousandfold error in the cost of translation.""" + add_ribosome(model, subunits, protein_rxn="r_protein") + + translation = _translation(model) + for unit in subunits.uniprot_ids: + coefficient = translation.metabolites[model.metabolites.get_by_id(f"prot_{unit}")] + assert coefficient == pytest.approx(-EXPECTED_COEFFICIENT) + assert EXPECTED_COEFFICIENT == pytest.approx(1.5873, abs=1e-4) + + +# --- what it does to the model --------------------------------------- + +def test_growth_is_unchanged_when_ribosomes_are_affordable(model, subunits): + before = model.optimize().objective_value + add_ribosome(model, subunits, protein_rxn="r_protein") + after = model.optimize().objective_value + assert after == pytest.approx(before) + + +def test_the_ribosome_is_paid_for_out_of_the_protein_pool(model, subunits): + before = model.optimize().fluxes["prot_pool_exchange"] + add_ribosome(model, subunits, protein_rxn="r_protein") + after = model.optimize().fluxes["prot_pool_exchange"] + + growth = model.optimize().objective_value + expected = before + len(subunits.uniprot_ids) * growth * EXPECTED_COEFFICIENT + assert after == pytest.approx(expected) + + +def test_without_ribosomes_there_is_no_protein_and_no_growth(model, subunits): + """The point of the rewiring: protein cannot bypass translation.""" + add_ribosome(model, subunits, protein_rxn="r_protein") + for unit in subunits.uniprot_ids: + model.reactions.get_by_id(f"usage_prot_{unit}").upper_bound = 0.0 + + solution = model.optimize() + assert (solution.objective_value or 0.0) == pytest.approx(0.0, abs=1e-9) + + +def test_a_capped_subunit_limits_growth_in_proportion(model, subunits): + """Half the ribosome, half the protein: growth 5 becomes 2.5.""" + add_ribosome(model, subunits, protein_rxn="r_protein") + full = model.optimize().objective_value + for unit in subunits.uniprot_ids: + model.reactions.get_by_id(f"usage_prot_{unit}").upper_bound = ( + 0.5 * full * EXPECTED_COEFFICIENT + ) + + assert model.optimize().objective_value == pytest.approx(0.5 * full) + + +def test_the_result_reports_what_it_added(model, subunits): + result = add_ribosome(model, subunits, protein_rxn="r_protein") + assert set(result.enzymes_added) == {"R1", "R2"} + assert result.amino_acid_demand == pytest.approx(TINY_AA_PER_PROTEIN) + assert result.kcat == pytest.approx(5.25) + assert result.reactions_added + + +# --- refusals --------------------------------------------------------- + +def test_adding_the_ribosome_twice_is_refused(model, subunits): + add_ribosome(model, subunits, protein_rxn="r_protein") + with pytest.raises(ValueError, match="already"): + add_ribosome(model, subunits, protein_rxn="r_protein") + + +def test_an_empty_subunit_set_is_refused(model): + empty = RibosomeSubunits(uniprot_ids=[], genes={}, masses={}, mean_abundance={}) + with pytest.raises(ValueError, match="no ribosomal subunits"): + add_ribosome(model, empty, protein_rxn="r_protein") + + +def test_a_pseudoreaction_without_protein_is_refused(model, subunits): + model.reactions.get_by_id("r_protein").subtract_metabolites( + {model.metabolites.get_by_id("protein"): 1.0} + ) + with pytest.raises(ValueError, match="protein"): + add_ribosome(model, subunits, protein_rxn="r_protein") + + +# --- selecting the core ---------------------------------------------- + +def _ribosome_table(ids): + return pd.DataFrame( + { + "uniprot": ids, + "gene_name": [f"N{i}" for i in ids], + "gene": [f"Y{i}" for i in ids], + "mass": [30000.0] * len(ids), + "sequence": ["M"] * len(ids), + } + ) + + +def _proteomics(values): + return pd.DataFrame( + { + "Protein.IDs": list(values), + "Gene": [f"Y{i}" for i in values], + "CN4_1_abs": [v[0] for v in values.values()], + "CN4_2_abs": [v[1] for v in values.values()], + } + ) + + +def test_a_subunit_below_the_threshold_is_left_out(): + table = _ribosome_table(["A", "B"]) + data = _proteomics({"A": (1e-4, 1e-4), "B": (1e-6, 1e-6)}) + selected = core_subunits(table, data, threshold=MIN_MEAN_ABUNDANCE) + assert selected.uniprot_ids == ["A"] + + +def test_a_subunit_that_was_never_measured_is_left_out(): + table = _ribosome_table(["A", "B"]) + data = _proteomics({"A": (1e-4, 1e-4)}) + assert core_subunits(table, data).uniprot_ids == ["A"] + + +def test_a_subunit_measured_only_as_missing_is_left_out(): + table = _ribosome_table(["A"]) + data = _proteomics({"A": (np.nan, np.nan)}) + assert core_subunits(table, data).uniprot_ids == [] + + +def test_selection_averages_over_every_replicate_of_every_condition(): + """A subunit seen strongly in one condition and not at all in another + still belongs to the ribosome; the average decides, not one column.""" + table = _ribosome_table(["A"]) + data = _proteomics({"A": (2e-5, 2e-6)}) + assert core_subunits(table, data).uniprot_ids == ["A"] + assert core_subunits(table, data).mean_abundance["A"] == pytest.approx(1.1e-5) + + +# --- against the real data ------------------------------------------- + +@pytest.mark.slow +def test_the_real_ribosome_table_reads(): + table = read_ribosome() + assert {"uniprot", "gene", "mass", "sequence"} <= set(table.columns) + assert len(table) == 253 + assert table["mass"].min() > 1000 + + +@pytest.mark.slow +def test_the_core_ribosome_is_selected_from_the_real_data(): + from overflow.proteomics import read_proteomics + + selected = core_subunits(read_ribosome(), read_proteomics()) + assert len(selected) == 48 + assert all(selected.mean_abundance[u] >= MIN_MEAN_ABUNDANCE for u in selected.uniprot_ids) + + +@pytest.mark.slow +@pytest.mark.parametrize("condition", ["CN4", "CN22", "CN38", "CN75", "hGR"]) +def test_the_core_is_exactly_what_the_matlab_run_added(condition): + """The GECKO 2 ecModel carries 816 enzymes and each finished + condition model carries 864. Those 48 are the subunits ribosome.m + selected, and this selection reproduces them one for one.""" + import scipy.io as sio + + from overflow.config import ROOT + from overflow.proteomics import read_proteomics + + models = ROOT / "legacy_matlab" / "models" + base = sio.loadmat( + models / "ecYeastGEM.mat", struct_as_record=False, squeeze_me=True + )["ecModel"].enzymes + built = sio.loadmat( + models / f"ecModel_P_{condition}.mat", struct_as_record=False, squeeze_me=True + )[f"ecModelP_{condition}"].enzymes + added = set(map(str, built)) - set(map(str, base)) + + selected = core_subunits(read_ribosome(), read_proteomics()) + assert set(selected.uniprot_ids) == added + + +@pytest.mark.slow +def test_the_real_protein_pseudoreaction_sets_the_elongation_cost(ec_model): + demand = amino_acid_demand(ec_model) + assert demand == pytest.approx(4.191011, rel=1e-5) + assert translation_kcat(ec_model) == pytest.approx(10.5 / demand) + + +# --- applying the measurements --------------------------------------- + +def test_a_measured_subunit_is_capped_at_its_measurement(model, subunits): + from overflow.ribosome import constrain_subunits + + add_ribosome(model, subunits, protein_rxn="r_protein") + growth = model.optimize().objective_value + generous = 2 * growth * EXPECTED_COEFFICIENT + + result = constrain_subunits(model, {"R1": generous, "R2": generous}) + assert result.adjusted == [] + assert model.reactions.get_by_id("usage_prot_R1").upper_bound == generous + assert model.ec.concs[model.ec.enzymes.index("R1")] == generous + + +def test_a_subunit_measured_below_what_translation_needs_is_raised(model, subunits): + """The cell was growing, so a measurement that would stop it is + taken as an underestimate rather than as a limit.""" + from overflow.ribosome import constrain_subunits + + add_ribosome(model, subunits, protein_rxn="r_protein") + growth = model.optimize().objective_value + needed = growth * EXPECTED_COEFFICIENT + + result = constrain_subunits(model, {"R1": 0.1 * needed, "R2": 10 * needed}) + assert result.adjusted == ["R1"] + assert result.required[0] == pytest.approx(1.01 * needed) + assert model.optimize().objective_value == pytest.approx(growth, rel=1e-6) + + +def test_an_unmeasured_subunit_keeps_drawing_on_the_pool(model, subunits): + from overflow.ribosome import constrain_subunits + + add_ribosome(model, subunits, protein_rxn="r_protein") + constrain_subunits(model, {"R1": 1.0}) + assert model.reactions.get_by_id("usage_prot_R2").upper_bound == 1000.0 + + +# --- reading the condition's subunit measurements -------------------- + +def test_subunit_abundances_are_filtered_like_any_other_protein(): + """A subunit belongs to the core on its average across conditions, + but its cap in one condition still has to survive that condition's + own quality filter.""" + from overflow.build_ribosome import subunit_abundances + from overflow.config import load_conditions + + table = pd.DataFrame( + { + "Protein.IDs": ["A", "B"], + "Gene": ["YA", "YB"], + "CN4_1_abs": [1e-5, 1e-7], + "CN4_2_abs": [1.1e-5, 9e-5], + "CN4_3_abs": [1.05e-5, 1e-7], + } + ) + selected = RibosomeSubunits( + uniprot_ids=["A", "B"], + genes={"A": "YA", "B": "YB"}, + masses={"A": 30000.0, "B": 30000.0}, + mean_abundance={"A": 1e-5, "B": 3e-5}, + ) + result = subunit_abundances( + load_conditions()["CN4"], selected, table=table, + masses={"A": 30000.0, "B": 30000.0}, + ) + assert set(result) == {"A"}, "B varies more than its median and should drop out" + assert result["A"] > 0 + + +def test_subunit_abundances_are_milligrams(): + from overflow.build_ribosome import subunit_abundances + from overflow.config import load_conditions + + table = pd.DataFrame( + { + "Protein.IDs": ["A"], + "Gene": ["YA"], + "CN4_1_abs": [1e-5], + "CN4_2_abs": [1e-5], + "CN4_3_abs": [1e-5], + } + ) + selected = RibosomeSubunits( + uniprot_ids=["A"], genes={"A": "YA"}, masses={"A": 30000.0}, + mean_abundance={"A": 1e-5}, + ) + result = subunit_abundances( + load_conditions()["CN4"], selected, table=table, masses={"A": 30000.0} + ) + assert result["A"] == pytest.approx(1e-5 * 30000.0) diff --git a/tests/tiny.py b/tests/tiny.py index a89eaf1..ad0c668 100644 --- a/tests/tiny.py +++ b/tests/tiny.py @@ -104,3 +104,41 @@ def tiny_ec_model(pool: float = 100.0, supply: float = 10.0) -> EcModel: ec_model.ec.validate() return ec_model +#: Amino acids consumed per unit of protein in the tiny model's protein +#: pseudoreaction. Two, so that the translation kcat works out to a +#: number that is easy to check by hand. +TINY_AA_PER_PROTEIN = 2.0 + + +def tiny_translation_model(pool: float = 200.0, supply: float = 10.0) -> EcModel: + r"""The tiny model with a protein pseudoreaction in the way of growth. + + supply -> S -+- R1 (E1) -+-> P -2P-> protein -> biomass + \- R2 (E2) -/ + + Growth now needs protein, and protein needs the protein + pseudoreaction, which is where the ribosome gets inserted. + """ + model = tiny_ec_model(pool=pool, supply=supply) + inner = cobra.Model("tiny_translation") + inner.add_metabolites( + [cobra.Metabolite("protein", name="protein", compartment="c")] + ) + protein = inner.metabolites.get_by_id("protein") + model.add_metabolites([protein]) + + product = model.metabolites.get_by_id("P") + biomass = model.metabolites.get_by_id("biomass") + + protein_rxn = cobra.Reaction( + "r_protein", name="protein pseudoreaction", lower_bound=0.0, upper_bound=1000.0 + ) + protein_rxn.add_metabolites({product: -TINY_AA_PER_PROTEIN, protein: 1.0}) + + growth = model.reactions.get_by_id("BIO") + growth.subtract_metabolites({product: growth.metabolites[product]}) + growth.add_metabolites({protein: -1.0}) + + model.add_reactions([protein_rxn]) + model.objective = "BIO" + return model From c8a75e5bf83e6dbdc3ce939d9538f43d5dc274f0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 01:10:28 +0200 Subject: [PATCH 2/4] feat: add the ribosome to each condition model Measured subunit abundances cap the subunits, except where the model needs more of one than was measured: the cell was demonstrably growing, so such a measurement is treated as an underestimate rather than as a limit on translation. --- src/overflow/build_ribosome.py | 201 +++++++++++++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 src/overflow/build_ribosome.py diff --git a/src/overflow/build_ribosome.py b/src/overflow/build_ribosome.py new file mode 100644 index 0000000..58675eb --- /dev/null +++ b/src/overflow/build_ribosome.py @@ -0,0 +1,201 @@ +"""Adding the ribosome to each condition model. + +Run after ``overflow.build``, on the models it wrote: +``python -m overflow.build_ribosome``. +""" +from __future__ import annotations + +import argparse +from dataclasses import dataclass +from pathlib import Path +from typing import Optional + +import numpy as np +import pandas as pd +from geckopy import load_ec_model, save_ec_model + +from overflow.adapter import build_adapter +from overflow.config import ( + BIO_RXN, + CONDITION_ORDER, + MODELS_DIR, + POOL_RXN, + RESULTS_DIR, + Condition, + load_conditions, +) +from overflow.proteomics import ( + filter_prot_data, + molecular_masses, + read_proteomics, + replicate_matrix, + to_mass, +) +from overflow.ribosome import ( + RibosomeResult, + RibosomeSubunits, + SubunitConstraints, + add_ribosome, + constrain_subunits, + core_subunits, + read_ribosome, +) + + +@dataclass +class RibosomeBuild: + condition: str + model: object = None + ribosome: Optional[RibosomeResult] = None + constraints: Optional[SubunitConstraints] = None + growth_before: float = 0.0 + growth_after: float = 0.0 + pool_before: float = 0.0 + pool_after: float = 0.0 + ribosome_mass: float = 0.0 + + +def subunit_abundances( + condition: Condition, + subunits: RibosomeSubunits, + table: Optional[pd.DataFrame] = None, + masses: Optional[dict[str, float]] = None, +) -> dict[str, float]: + """Measured subunit abundances for one condition [mg/gDW]. + + Filtered on the same terms as every other protein, so a subunit that + was measured unreliably in this condition carries no cap here even + though its average across conditions put it in the core. + """ + if table is None: + table = read_proteomics() + if masses is None: + masses = molecular_masses() + + ids, replicates = replicate_matrix(table, condition) + filtered = filter_prot_data(ids, replicates) + wanted = set(subunits.uniprot_ids) + keep = [(p, v) for p, v in zip(filtered.uniprot_ids, filtered.abundances) if p in wanted] + if not keep: + return {} + proteins = [p for p, _ in keep] + values = np.array([v for _, v in keep], dtype=float) + return dict(zip(proteins, to_mass(proteins, values, masses))) + + +def add_ribosome_to_condition( + condition: Condition, + models_dir: Path = MODELS_DIR, + subunits: Optional[RibosomeSubunits] = None, + table: Optional[pd.DataFrame] = None, + masses: Optional[dict[str, float]] = None, + verbose: bool = True, +) -> RibosomeBuild: + """Load a condition model, put the ribosome in it and constrain it.""" + if subunits is None: + subunits = core_subunits(read_ribosome(), table if table is not None else read_proteomics()) + + path = models_dir / f"ecModel_P_{condition.name}.yml" + model = load_ec_model(str(path), adapter=build_adapter(condition)) + result = RibosomeBuild(condition=condition.name, model=model) + + solution = model.optimize() + result.growth_before = float(solution.fluxes[BIO_RXN]) + result.pool_before = float(solution.fluxes[POOL_RXN]) + + result.ribosome = add_ribosome(model, subunits) + result.constraints = constrain_subunits( + model, subunit_abundances(condition, subunits, table, masses), subunits.uniprot_ids + ) + + solution = model.optimize() + result.growth_after = float(solution.fluxes[BIO_RXN]) + result.pool_after = float(solution.fluxes[POOL_RXN]) + result.ribosome_mass = sum( + abs(float(solution.fluxes[f"usage_prot_{u}"])) for u in subunits.uniprot_ids + ) + + if verbose: + print( + f"[{condition.name}] {len(subunits)} subunits, kcat " + f"{result.ribosome.kcat:.4f} /s ({result.ribosome.amino_acid_demand:.4f} " + f"amino acids per protein)\n" + f" {len(result.constraints.adjusted)} subunits raised to what " + f"translation needs\n" + f" growth {result.growth_before:.5f} -> {result.growth_after:.5f}; " + f"protein pool {result.pool_before:.1f} -> {result.pool_after:.1f} mg/gDW, " + f"of which ribosome {result.ribosome_mass:.1f}", + flush=True, + ) + return result + + +def summary_table(results: list[RibosomeBuild]) -> pd.DataFrame: + return pd.DataFrame( + [ + { + "condition": r.condition, + "subunits": len(r.ribosome.subunits), + "translation_kcat": round(r.ribosome.kcat, 4), + "aa_per_protein": round(r.ribosome.amino_acid_demand, 4), + "subunits_raised": len(r.constraints.adjusted), + "growth_before": round(r.growth_before, 5), + "growth_after": round(r.growth_after, 5), + "pool_before": round(r.pool_before, 1), + "pool_after": round(r.pool_after, 1), + "ribosome_mg_gDW": round(r.ribosome_mass, 1), + "ribosome_fraction": round(r.ribosome_mass / r.pool_after, 4) + if r.pool_after + else float("nan"), + } + for r in results + ] + ) + + +def main(argv: Optional[list[str]] = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("conditions", nargs="*", default=list(CONDITION_ORDER)) + parser.add_argument("--models-dir", type=Path, default=MODELS_DIR) + parser.add_argument("--results-dir", type=Path, default=RESULTS_DIR) + parser.add_argument( + "--solver", + help="LP solver to use, e.g. gurobi or glpk (default: cobrapy's)", + ) + args = parser.parse_args(argv) + + if args.solver: + import cobra + + cobra.Configuration().solver = args.solver + + conditions = load_conditions() + table = read_proteomics() + masses = molecular_masses() + subunits = core_subunits(read_ribosome(), table) + print(f"core ribosome: {len(subunits)} of {subunits.n_candidates} subunits", flush=True) + + generation = args.results_dir / "modelGeneration" + generation.mkdir(parents=True, exist_ok=True) + + results = [] + for name in args.conditions or CONDITION_ORDER: + result = add_ribosome_to_condition( + conditions[name], models_dir=args.models_dir, + subunits=subunits, table=table, masses=masses, + ) + save_ec_model(result.model, args.models_dir / f"ecModel_P_{name}_ribosome.yml") + result.constraints.table().to_csv( + generation / f"modifiedRibosomeSubunits_{name}.tsv", sep="\t", index=False + ) + results.append(result) + + summary = summary_table(results) + summary.to_csv(generation / "ribosomeSummary.tsv", sep="\t", index=False) + print() + print(summary.to_string(index=False)) + return 0 + + +if __name__ == "__main__": # pragma: no cover + raise SystemExit(main()) From c2918925a1d4adfb750be6d744bda76219cfccad Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 01:10:57 +0200 Subject: [PATCH 3/4] docs: how to add the ribosome --- README.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/README.md b/README.md index 8d37f1b..1b933ae 100644 --- a/README.md +++ b/README.md @@ -88,6 +88,18 @@ Options change what the models are asked to account for: default: the model has no carbon to spare at the measured glucose uptake, and CN4 then falls short of its dilution rate. +## Adding the ribosome + +```bash +python -m overflow.build_ribosome +``` + +Writes `models/ecModel_P__ribosome.yml`. The protein pseudoreaction is +split so that protein comes out of a `translation` reaction catalysed by the +ribosomal subunits, at 10.5 amino acids per second per ribosome. The core is +the 48 subunits whose average abundance across all conditions reaches +1e-5 mmol/gDW. + ## Tests ```bash From 67eaec36c173114108646a3bf5b205b8834a1c6c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 20 Sep 2026 01:16:38 +0200 Subject: [PATCH 4/4] test: pin that the ribosome is charged per unit of protein made --- tests/test_ribosome.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/test_ribosome.py b/tests/test_ribosome.py index 59a698c..2da96c6 100644 --- a/tests/test_ribosome.py +++ b/tests/test_ribosome.py @@ -373,3 +373,18 @@ def test_subunit_abundances_are_milligrams(): load_conditions()["CN4"], selected, table=table, masses={"A": 30000.0} ) assert result["A"] == pytest.approx(1e-5 * 30000.0) + + +def test_ribosome_cost_is_proportional_to_how_fast_protein_is_made(model, subunits): + """Translation is charged per unit of protein, so halving the supply + halves the ribosome. Across the real conditions this shows up as + 11.5 mg/gDW at D=0.1 and 33.3 at D=0.29, a ratio of 2.90 against a + dilution-rate ratio of 2.9.""" + add_ribosome(model, subunits, protein_rxn="r_protein") + full = model.optimize() + full_usage = abs(full.fluxes["usage_prot_R1"]) + + model.reactions.get_by_id("EX_S").upper_bound /= 2 + half = model.optimize() + assert half.objective_value == pytest.approx(full.objective_value / 2) + assert abs(half.fluxes["usage_prot_R1"]) == pytest.approx(full_usage / 2)