Repository navigation
Phase 3: condition-specific proteome-constrained ecModels - #5
Merged
Merged
Conversation
filter_prot_data ports GECKO 2's filter_ProtData: a protein is kept when it is detected in two thirds of the replicates and its relative standard deviation is below 1, and is represented by mean + 1.96 SD rather than by its mean. A protein with a missing replicate is dropped, because its spread cannot be evaluated. Abundances are converted from mmol/gDW to mg/gDW with UniProt molecular masses, which covers every measured protein and removes the mean-mass fallback the MATLAB implementation used for proteins outside the model. P61830 is listed once per gene of the duplicated histone; the reader collapses it and fails if the rows ever disagree.
edkerk
force-pushed
the
feat/py-phase3-condition-models
branch
from
September 19, 2026 23:08
7858de3 to
4134be3
Compare
This was referenced Sep 19, 2026
Merged
Phase 2: proteomics reading, filtering and unit conversion
Base automatically changed from
feat/py-phase2-data
to
feat/py-phase1-skeleton
September 25, 2026 10:59
The model and its conventional GEM now come from SysBioChalmers/ecModels (branch gecko4, yeast-GEM 9.1.1) instead of the GECKO tutorial model. Provenance records the commit, blob and checksum.
Protein content and growth-associated maintenance become explicit policies. Maintenance is written by shifting the biomass proton coefficient rather than setting it: yeast-GEM's biomass releases three fewer protons than it hydrolyses ATP, and setting the proton to the maintenance value instead adds three protons to every gram of biomass. That is invisible in the stoichiometry and costs about one percent of growth, which is the difference between CN4 reaching its dilution rate and missing it.
Undetected byproducts are blocked, detected ones keep ten percent of headroom, and the chemostat is set by minimising uptake first and protein second, so that protein cannot be bought with extra carbon.
Proteomics recovers membrane-bound subunits unevenly, leaving a complex whose subunits disagree about how much of it there is. Each subunit is set to its stoichiometric share of one abundance per complex, in molar units, because the stoichiometry it is divided by is molar.
A minimum-usage pass raises abundances the model cannot run the condition on, then limiting caps are released one at a time until the model reaches the dilution rate and pulled back to the usage a minimum-protein solution gives them. The limiting cap is found from shadow prices where they are informative and by probing each binding cap where they are not. Shadow prices go silent on a degenerate optimum, which reads as though nothing were limiting: the search then grows the protein pool instead, and on CN4 it grew it 457 times without ever gaining growth. The tests run against a hand-built ecModel small enough that every quantity is worked out by hand.
The residual is relative over glucose, CO2 and oxygen together, so that no rate dominates the fit by being large.
One model per chemostat condition, written with the tables describing how it was built. A condition the model cannot reach fails with what it reached and what was already tried, rather than surfacing three steps later as an infeasible chemostat.
The pipeline scripts take --solver. The test suite defaults to GLPK: it is small enough not to need a commercial solver, and a licence that is busy or absent should not turn into a wall of failures that look like code.
A model file records objective coefficients but not the direction, so minimising protein came back as maximising it: every model this pipeline wrote reported the most protein the budget allowed rather than the least the condition needed, 334.7 against 61.9 mg/gDW for CN4. Written now as maximising the negative, which survives the round trip.
The pipeline constrains uptake, growth and the measured byproducts, but leaves CO2 and oxygen free and every other secretion open, and minimises protein. Respiring carbon is the most protein-expensive way to dispose of it, so the model vents the surplus through whichever exit is cheapest instead: at CN38 a third of the carbon taken up leaves as acetaldehyde, and both gas rates fall short of their measurements. --fit-rates holds CO2, oxygen and the byproducts within a tolerance of their measurements, then raises the measured abundances by the least that makes that feasible. Each cap becomes 'usage <= measured + slack' and the slacks are minimised in fold-change, so the answer is the smallest departure from the proteomics that accounts for the rates rather than one enzyme released outright. An enzyme measured at almost nothing is not relaxed freely just because its absolute cost is small.
--uptake-flex sets the cap on glucose uptake; the 9.1.1 model needs 6.5% more than CN4 measured. --objective picks the smallest protein pool, or the smallest total flux over all or only the metabolic reactions.
edkerk
force-pushed
the
feat/py-phase3-condition-models
branch
from
September 25, 2026 11:33
e007e51 to
7577392
Compare
Merged
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stacked on #4. Builds one ecModel per chemostat condition: biomass and
maintenance set for the condition, measured rates applied, enzyme abundances
integrated and reconciled, maintenance fitted to the gas rates.
python -m overflow.buildwritesmodels/ecModel_P_<cond>.ymlplus the tablesdescribing how each was built. All five conditions build, with
Model
The model and its conventional GEM now come from
SysBioChalmers/ecModels, branch
gecko4at8802887: yeast-GEM 9.1.1, 1144 enzymes over 4850 enzyme-constrainedreactions, kcats tuned against measured growth and exchange fluxes. This replaces
the GECKO tutorial model that #3 vendored.
models/PROVENANCE.mdrecords thecommit, blob hash and yeast-GEM checksum. Everything below was built and measured
on it, and the tests that named the old model's size are updated.
The tutorial model had cytochrome c oxidase at 6000 /s (a kcat-matcher bug, see
#7); here it is 1376 /s.
Options
--fit-ratesholds CO2, oxygen and the byproducts within--rate-toleranceof their measurements and raises the measured abundances bythe least that makes that feasible: each cap becomes
usage <= measured + slackand the slacks are minimised in fold-change, so the answer is the smallest
departure from the proteomics that accounts for the rates. Minimising
milligrams instead would relax whichever enzyme is cheapest in absolute terms,
usually the one measured at almost nothing.
--uptake-flexcaps glucose uptake at that multiple of the measured rate(default 1.05). CN4 needs 1.08, see below.
--objectivechooses what the finished model minimises: the protein pool(
protein, default), or total flux at the dilution rate over every reaction(
flux) or over the metabolic reactions only (metabolic-flux).--scale-proteinrescales biomass to the measured protein content. Off bydefault, see the last section.
Fit to the measurements
Measured / GECKO 2 / this pipeline (built as above), mmol/gDW/h, and the relative
residual over glucose, CO2 and oxygen together. GECKO 2 values are read from
legacy_matlab/results/modelSimulation/.Mean residual 0.487 → 0.123. Every condition grows at exactly its dilution
rate.
The build without
--fit-rateswas not rerun on this model. On the earliertutorial model it left CO2 and oxygen far short at CN38 and CN75: with those
rates free, the protein-minimising objective disposed of surplus carbon through
unmeasured exits (36% of the carbon at CN38) instead of respiring it, and closing
the exits made the conditions infeasible. That is why
--fit-ratesexists. Thosenumbers belong to the previous model and are not repeated here.
CN4 needs more glucose than was measured
With protein unlimited and the measured byproducts held, the model reaches
0.0939 /h at CN4's measured glucose; growing at 0.1 /h needs 1.188 mmol/gDW/h
against 1.116 measured, 6.5% more. The 5% default is not enough, so CN4 fails
with a message naming the carbon budget.
--uptake-flex 1.08is used for allconditions; building the other four at 1.05 gives essentially the same results.
Maintenance follows from this: with no carbon to spare, the fitted NGAM sits at
its lower bound of 0 in every condition, and in CN22 and CN75 the scan finds no
feasible point at all (error
inf), so it carries no information here.The objective does not decide the pentose phosphate pathway
The enzyme-constrained solution puts no flux through the pentose phosphate
pathway (PP shunt) at CN38 and CN75, where the published analysis had about 20%
median capacity usage. To check whether that comes from minimising protein, the
same builds were run with the two flux objectives:
proteinfluxmetabolic-fluxAll three agree, so it is not the objective. NADPH demand is the same as in the
previous model (1.9–2.0 mmol/gDW/h at CN38 and CN75) and is met by cytosolic
isocitrate dehydrogenase (IDP2) and aldehyde dehydrogenase instead of the shunt.
The kcats of ZWF1, GND1 and IDP2 are identical in both models. What blocks the shunt
is the measured cap of RKI1 (about 0.01 mg/gDW at CN38, all of it used for nucleotide
precursors), with SOL3 next: lifting them gives the shunt 0.21 at CN38 under every
objective, and at CN75 it is also limited by the CO2 band. This is accepted. The model
is built from the measured abundances, and it need not carry flux through the shunt
or match the published shunt usage.
Two bugs worth reading about
The biomass proton.
raven_toolbox.biomass.set_gamsets every listedcofactor's coefficient to ±GAM, and yeast-GEM's biomass carries H+ at 52.3
against ATP at 55.3 — three protons fewer, from the NADPH it also consumes.
Passing
H+in the cofactor list therefore rewrote 52.3 to 55.3 even when theGAM was unchanged, silently adding three protons per gram of biomass. It cost
CN4 about 1.3% of growth and made it infeasible. Maintenance is now written by
shifting the proton, and a test pins the offset. Worth a look on the
raven-toolbox side — the sharp edge is that the no-op case is not a no-op.
Shadow prices go silent on a degenerate optimum. Finding the limiting
enzyme from
get_conc_control_coeffsalone fails when several caps bind atonce and none carries a non-zero dual: the search concludes nothing is limiting
and grows the protein pool instead. On CN4 it grew it 457 times, to 92× its
size, gaining nothing. The original used a finite-difference probe
(
findLimitingUBs, flexFactor 100) and was never exposed to this. Shadowprices are now the fast path and the probe the fallback, plus a stall detector
so a pool that is not the constraint cannot be inflated indefinitely.
The one judgement call
generateProtModels.mrescales biomass protein to the measured content. Doingthat here costs carbon the model does not have: CN4 already needs 6.5% more
glucose than measured before any rescaling. Rescaling is therefore off by
default and available as
--scale-protein; GECKO 4's own yeast tutorialdisables the same step (
YeastGEMAdapter.changeProteinBiomass). The measuredprotein content still sets the protein pool either way.
Tests
114 in this branch, including the integration build of CN4. The flexibilization
ones run against a hand-built five-reaction ecModel where a unit of flux costs
exactly 10 mg of one enzyme or 5 mg of the other, so the expected minimum usage,
the released bound and the pool draw are all worked out by hand rather than
asserted from a previous run. The subset-parsimony test uses a two-route network
where minimising every reaction and minimising only the metabolic ones must pick
different routes.