Phase 6: random sampling - #8
Merged
Merged
Conversation
Two passes per condition, one with formate left out of the measured set so the sampler can say which byproducts the model makes unprompted, and the ATP and redox budget over the sampled means. yeast-GEM has merged glucokinase (r_4235) into hexokinase, so the glycolytic ATP expenditure is now read from r_0534 alone; the accounting takes the reaction sets as arguments, which is what lets it be checked against the published summary on that summary's own reaction set.
Screening every reaction for loop involvement is the dominant cost of a sampling run and depends only on the network, so it is kept across runs. Opening the arbitrary 1000 bounds to infinity, as the published analysis did, is now off by default: the screening that picks the random objectives runs on the finite bounds, so once they are opened an objective can turn out unbounded, and this sampler raises on that where RAVEN's returned no solution and moved on.
The screening is a flux variability analysis, which reads its process count from cobra's configuration rather than from the sampler argument. Left unset it ran on one core and took most of the run.
The published second pass minimised total flux within each draw to suppress loops. Here one numerically awkward draw in several thousand raises out of the parsimonious solve and aborts the whole run, where RAVEN's solveLP returned no solution and the loop simply retried, so it is off unless asked for.
The malate dehydrogenases cycled against each other at several hundred mmol/gDW/h, against 0.1 and 0.9 in the published result: ordinary variability analysis lets a reaction in a thermodynamically infeasible cycle reach the model's arbitrary 1000 bound, and sampling then draws from that room. Loop-free variability analysis puts the same pair at -2.1 to 2.9 and -4.6 to 4.5. Every flux distribution free of closed cycles already lies inside those bounds, so nothing sensible is excluded. It costs one loop-free variability analysis per condition, about a minute on six cores. Also: have multiprocessing fork rather than use a forkserver. Python 3.14 made forkserver the default and its helper process does not survive here, which killed runs partway through as a broken pipe.
The reactions worth using as random objectives are the ones that can carry flux without a closed cycle, which is what the loop-free variability analysis already reports. Deriving both from it costs nothing extra and, more to the point, makes the objective set belong to the condition being sampled: reusing one condition's set for another leaves reactions in it that are pinned here, and each time one is picked the sampler retries a hundred times before giving up on that draw. hGR spent thirty minutes on what takes the other conditions two.
The sampled means drift into high-flux routes without the second step the published analysis had: the pentose phosphate pathway and a reversed NADP isocitrate dehydrogenase carry several times the published flux. The step now tolerates a draw whose optimum cannot be pinned exactly. cobra picks the cyclic reactions with unseeded random LP weights, and one seed in a while runs for hours. The search now has a time limit and retries with a new seed.
edkerk
force-pushed
the
feat/py-phase5-enzyme-usage
branch
from
September 25, 2026 11:33
309c3a3 to
d6f163e
Compare
edkerk
force-pushed
the
feat/py-phase6-sampling
branch
from
September 25, 2026 11:33
3cd51e1 to
4c90c01
Compare
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
1 participant
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 #7. Samples the conventional yeast-GEM under the measured rates,
with and without formate in the measured set, and summarises the ATP and redox
budget over the sampled means.
python -m overflow.run_sampling --procs 12— about 4.5 minutes per condition at5000 samples per pass.
The ATP budget reproduces the published one
Published / this pipeline, mmol ATP/gDW/h:
Glucose uptake is identical to two to three figures, respiratory and
growth-associated ATP within a few percent. Maintenance comes out lower
throughout, which is where the remaining difference in the total sits; the model
needs more carbon for growth than was measured (#5), so less is left for it.
The redox shuttle now runs the published way:
Cytosolic malate dehydrogenase sits near 0 here where the published values are
−0.1 to −1.9.
Measured byproducts reproduce; the predicted one does not
Both passes hold the measured rates to within 5%, and the sampled means sit inside
those bands in all five conditions. The first pass leaves formate out of the
measured set and asks what the model secretes unprompted. The published answer was
formate, at 0.45 to 3.67 mmol/gDW/h. Here formate is much weaker (0.20 at CN22,
0.35 at CN75, 0.95 at hGR, nothing above the detection threshold at CN4 and CN38),
and glyoxylate (0.26 at CN22, 0.78 at CN75) and L-glycine (up to 1.0 at hGR)
appear. So the model still wants an overflow route, but a different one.
Exact parity where it can be had
The ATP accounting is arithmetic on a flux vector, so it is checked against the
published summary rather than described: fed the committed sampling means, it
reproduces
selectedFluxes.txtfor all five conditions, to that file's foursignificant figures.
Three things the port had to do differently
Minimising total flux in each draw. The published second pass minimised total
flux within each draw. Without it the sampled means drift into high-flux routes:
at CN38 the pentose phosphate pathway carries 2.47 (published 0.24) and NADP
isocitrate dehydrogenase runs backwards at −3.5 (published −0.08), turning NADPH
into NADH. With it they are 0.73 and −0.07 (at 300 samples). The step is on by
default (
--no-min-fluxskips it) and tolerates the occasional draw whoseoptimum cannot be pinned exactly, keeping the plain solution for it.
Loops. Ordinary variability analysis lets a reaction in a closed cycle reach
the model's arbitrary 1000 bound, so the malate dehydrogenases cycled against
each other at several hundred mmol/gDW/h. Reactions are now held to their
loop-free variability range before sampling, which puts that pair at −2.1 to 2.9
and −4.6 to 4.5. Every cycle-free flux distribution already lies inside those
bounds, so nothing sensible is excluded. The same analysis supplies the random
objectives, which has to be done per condition: a reaction pinned in one condition
wastes a whole draw when reused in another.
cobra picks the cyclic reactions with unseeded random LP weights, and how long
that takes depends on them: most seeds finish in a minute or two, an occasional
one ran for hours. The search now has a 15-minute limit and retries with a new
seed, up to four times.
Vertices, not interior.
random_objectivedraws polytope vertices, as thepublished analysis did, so these means are vertex means rather than an average
over the flux space.
Tests
182 tests in this branch, including the loop-free bounding, the retry of an
overrunning search and the parity check above.