Repository navigation
Conversation
This reverts commit c8ef390.
…e. In intracellular.pp added more comments that refers to unit measures. In Model, swtiched from primal to dual solver, and reverted to default settings. Added warning when a reaction is not found while setting the lower and upper bounds. Finally, switched intracellulr solver in the cell container to be called before the BioFVM solver.
…el to reduce model degeneracy and using primal if problem infeasible to make sure cell is dead
Every object rule named only its own .cpp as a prerequisite:
custom.o: ./custom_modules/custom.cpp
$(COMPILE_COMMAND) $(INC) -c ./custom_modules/custom.cpp
so editing a header rebuilt nothing. That is not merely a stale-build
annoyance here: addons/dFBA/src/dfba_intracellular.h defines the layout of
dFBAIntracellular, and sizeof changes when members are added to it or to the
dFBAModel it holds by value. An incremental `make` after such an edit produced
a binary whose translation units disagreed about that layout -- one object
allocating the cell's intracellular model at one size while another read
substrate_exchanges at an offset from a different version of the header.
The result was a read past the end of the object, a std::map iterated through a
pointer from unrelated memory, and a segfault three seconds into a run. It
disappeared on any clean rebuild and in any other build directory, which made
it look convincingly like latent undefined behaviour in the source. It was not;
the committed code is fine. Valgrind on the bad binary:
Invalid read of size 8 at post_update_intracellular(...)
Address 0x51bf3478 is 24 bytes after a block of size 1,632
Invalid read of size 8 at post_update_intracellular(...)
Address 0x708 is not stack'd, malloc'd or (recently) free'd
1,632 matches no layout in the tree: sizeof(dFBAIntracellular) is 2520 on
development and 4432 with the pFBA members. Three layouts in one binary.
The fix is -MMD -MP into COMPILE_COMMAND, `-include $(wildcard *.d)` to pull
the generated prerequisites in, and *.d removed by `clean` and gitignored.
Vendored libSBML and Clp move from -I to -isystem at the same time. They are
never edited, so -MMD should not record them: they were 304 of the 328 entries
in custom.d, and dropping them takes each dependency file from 40K to 12K.
-isystem also silences warnings originating inside them. addons/dFBA/src stays
on -I, since those are exactly the headers that change.
Verified by A/B on custom.o:
with fix without fix
touch dfba_intracellular.h rebuilt NOT rebuilt
touch coin/ClpSimplex.hpp not rebuilt not rebuilt
and the dependency file still lists all five addons/dFBA/src/dfba_*.h.
Clean builds and the gate pass: cancer_core_metabolism builds with no errors
and runs the slice with no memory errors; the addon test binaries build; and
dfba_unit_test builds and runs Unit_test_mass_conversion.xml to completion
through 25 save points.
Two pre-existing staging warts this made visible but does NOT touch, because
they live in the root Makefile, which is outside the agent modification
boundary: the staging targets are not .PHONY, so `make cancer-core-metabolism`
is skipped when a binary of that name exists, and the target ends in an
unguarded `mkdir ./scripts/` that fails whenever the directory is already
there.
Co-Authored-By: Claude Opus 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01PLyUEkoZC8GSqSdZgzWB7o
Header dependency tracking in the project Makefiles, so editing an addon header rebuilds the objects that include it. Fixes the mixed-class-layout binaries that produced a segfault in post_update_intracellular which looked like a source bug but was a stale-object build artifact. Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01PLyUEkoZC8GSqSdZgzWB7o
Brings in the header dependency tracking hotfix. The only conflict was .gitignore, where both sides had appended a block; both are kept -- the dFBA test binaries and artifacts from this branch, and *.d from the hotfix.
Callers no longer orchestrate the two stages. optimize() performs them and
returns one consolidated solution:
status stage one's -- a cell lives or dies on FBA feasibility
objective_value stage one's, the growth rate; never a total flux
fluxes stage two's when it solved, else stage one's
secondary_status "none" / "optimal" / "infeasible" / "unknown"
total_absolute_flux total |flux| of the distribution actually returned
optimizeParsimonious() is gone as public API, replaced by a private
solveParsimoniousStage() that writes in place. There is now only one solution
object, so a value that is a total flux can no longer appear where a caller
expects growth -- it is never put there. dFBAIntracellular::update() is one
optimize() call and run_secondary_objective() is removed.
This also fixes a real gap in the previous split: optimize_for_objective() and
custom_optimization never called run_secondary_objective(), so pFBA was
silently skipped on those paths. It cannot be now.
Two defects found while verifying, both by sweeping the bound space rather than
by reading the code:
secondary_status and total_absolute_flux were only set inside the optimal
branch, so an infeasible or unknown solve left whatever the previous call had
put there. A cell that failed stage one reported secondary_status "optimal" and
a total_absolute_flux from an earlier timestep. They are now reset for every
outcome at the top of optimize().
The second stage now retries once from a clean basis before giving up.
Warm-starting it across timesteps can drive Clp into status 4, "stopped due to
errors", on a problem that is perfectly well posed: neither primal nor dual
infeasible, tens of iterations rather than the 10000 limit, plausible objective
value. Measured on MCF7 core, 3 of 81 bound combinations failed with one model
reused across all of them and 0 of 81 with a fresh model per combination, at
identical bounds. Sweeping fraction_of_optimum from 1.0 down to 0.99 changed
nothing, which ruled out the objective floor sitting on the optimal face -- the
obvious suspect and the wrong one. allSlackBasis(true) plus a re-solve fixes
every case, and runs only on failure, so the warm start is untouched normally.
It mattered more than a status code: these failures cluster where the flux
distribution is worst conditioned, which is where degeneracy is worst, so the
silent fallback skipped the second stage for exactly the cells that needed it.
The affected distributions had total |flux| 3011.8 against the 11.8 the second
stage produces -- 256x less parsimonious.
Verification:
- test_pfba 19/19 on E. coli core and on MCF7 core, including a new
81-combination sweep asserting the second stage solves for every feasible
bound set on a reused model, which is the regression test for the above
- cobrapy cross-check on E. coli core still exact: total flux agrees to
4.4e-16, per-reaction max difference 5.3e-14 over 95 reactions, none above
1e-6. On MCF7 core the totals agree to 6.4e-06 with the addon marginally
LOWER, and 61 of 374 reactions differ above 1e-6 -- expected on a
degenerate problem where distinct distributions tie on total flux
- the unified reversible flux is asserted, not assumed: Sum(p+m) from the LP
equals Sum|v| read back off the reactions, which holds only if each
original column really is p - m with at most one of the pair nonzero
- stage two costs 0.115 ms on E. coli core and 0.471 ms on MCF7 core,
1000/1000 solved; the warm start is unchanged at 0.005 iterations/solve
- on MCF7 core the second stage cuts total |flux| from 3026.8 to 21.5, a
99.3% reduction, so the arbitrary vertex plain FBA returned carried 140x
more flux than necessary at the same growth rate
Manuscript-sensitive: this changes the flux distribution reported by any
configuration with <secondary_objective> enabled, and the stale-status fix
changes what a failed stage reports. cancer_core_metabolism is the validation
case to rerun.
The FVA analysis that produced the alternative-optima evidence was tracked only
on the unmerged `pfba` branch: `development` and `feature/pfba` carried no
fva_* files at all, while experiments_fva/README.md claimed the canonical copy
lived in cancer_core_metabolism/scripts/. The copies inside each experiment
folder are snapshots in a gitignored directory, so the code behind the
reviewer evidence was effectively untracked. fva_degeneracy.py, fva_figures.py
and fva_degeneracy.ipynb now live here.
reviewer_response_figures.ipynb builds the figure for the alternative-optima
comment. Four panels, each chosen for the job its data does rather than by
preference:
A degeneracy is system-dependent -- one row per use case, one dot per
exchange. Deliberately not a boxplot: two of the three systems are at
exactly zero, and a box drawn at zero reads as missing data rather than
as a result.
B where it exists it is spatially structured -- the five widest exchanges
as median and interquartile ribbon against depth, the other eighteen as
one grey band rather than eighteen more colours.
C growth is preserved -- the Dirichlet-clamped band nearest the boundary is
the control, since feedback cannot reach it.
D what pFBA changes -- diverging log2 ratio, because direction matters and
zero is meaningful.
S1, supplementary: exchange by time as a heatmap on a single sequential hue.
A handful of byproducts carry essentially all the degeneracy, and from early
on rather than only once the core turns hypoxic.
Numbers the figure reports, tumour slice at 4320 min, pFBA against standard
FBA at identical omp_num_threads:
mean growth 0.00421475 -> 0.00435489 (+3.33%)
total |exchange flux| 0.455998 -> 0.415824 (-8.81%)
glucose, O2 ratios 0.994, 1.005 -- the primary uptakes barely move
glycine, threonine 0.408, 0.408 -- the freedom is in the byproducts
Three things worth knowing about how it is built:
Panel D plots an exchange only when BOTH runs carry flux above 1e-3
mmol/gDW/h. Requiring only one of them keeps cases like glutamate, which moves
from 1.7e-4 to 1.4e-3 -- a ratio of 8.4 that would be the largest bar in the
panel, off a denominator that is effectively zero against a glucose uptake of
2.5e-2. The seven excluded exchanges are printed below the figure rather than
silently dropped.
Panels B and S1 report the ABSOLUTE interval width while panel A reports the
relative one. The first version used the relative measure throughout, which
sent panel B above 20 for a metabolite panel A showed at 2.7, because dividing
per-row by a small |flux| inflates without bound. Two different denominators
for one word is a figure that argues with itself; the panels now carry two
units and say so.
Panels C and D compare depth-binned distributions, not matched cells: over
4320 min the two runs diverge through BioFVM feedback and division
stochasticity, ending at 308 and 309 cells with different IDs.
The palette is Okabe-Ito reordered so the weakest colour-vision-deficient
adjacent pair is not adjacent, checked with a validator rather than by eye;
every series that uses the low-contrast orange is directly labelled.
Figures are written to experiments_fva/figures_reviewer_response/ as PNG and
PDF. Outputs are stripped from the committed notebook.
The dFBA uptake cap is meant to stop cells taking more of a substrate than their voxel holds. Two defects make it ineffective, and this test fails on both of them at the current development head: * update_dfba_inputs converts the voxel content to mmol by MULTIPLYING by liter_micron_cubes_conversion (1e15). One um^3 is 1e-15 L, so this needs a division; the cap is off by 1e30 and never engages. Visible with a single cell (n=1: glucose -1.8e-2 mM). * The cap is applied to each cell independently, so n cells sharing a voxel can draw n times its content (n=8: glucose -0.62 mM, exactly 0.165 - 8 * 0.0981). The two configs isolate glucose per voxel (D = 0) and put n cells at one position, so the check exposes the cap arithmetic alone. Not yet fixed; nothing in addons/dFBA/src is changed here. Co-Authored-By: Claude Sonnet 5.5 <[email protected]>
total_substrate converted mM * um^3 to mmol with a multiplication by liter_micron_cubes_conversion (1e15); one um^3 is 1e-15 L so it must be a division. The cap was ~1e30 too large and never engaged, even for a single cell. Verified: * voxel-sharing test, n=1 control: -1.8e-2 mM -> -5e-6 mM (PASS). * n=8 still FAILS (-0.62 mM): every cell's own request fits in the voxel, so a per-cell cap cannot engage. That is the separate no-sharing defect and is not addressed here. * deterministic sweep (D=0, 1 thread, S0 0.0005..1 mM, n=1/3/8): the fixed build is never worse than the unfixed one. Where the cap engages the overdraw is exactly -(n-1) x voxel content, against -4x to -42x unfixed. * mass-conversion gate with dense snapshots through the depletion window (n=6 each): worst excursion -0.066..-0.145 mM fixed vs -0.072..-0.190 mM unfixed. Peaks are smaller, but the fixed runs spend more snapshots below -1e-4 mM (145-202/351 vs 52-69/351). A coarse (10 min) snapshot comparison suggested the fix made the minimum worse; that was an artefact of sampling the two builds at different phases of the depletion and is not a property of the change. Multi-cell overdraw (-(n-1) x content) remains until a per-voxel share is applied. That is the next change. Co-Authored-By: Claude Sonnet 5.5 <[email protected]>
The uptake cap compared each cell's request with the whole voxel content, so n cells in one voxel could draw n times what it holds and BioFVM, which applies the sum of the net export rates explicitly, drove the concentration negative. The cap is now proportional: each live dFBA cell's bound is scaled by min(1, content / sum of the requests of all live dFBA cells in its voxel). The request is the cell's own Michaelis-Menten uptake over one dFBA window, so a cell that cannot take the substrate (Vmax = 0) asks for nothing and does not dilute the others. With a single cell it reduces to the previous cap. The request is computed in one helper shared by the per-cell bound and the table. The per-voxel table is shared state: it is built once per simulated time by the first cell to need it, under a named omp critical section with an atomic double-checked time stamp, and is read-only afterwards. It is not per-cell state, so the four ownership paths are unaffected. Assumes intracellular_dt equals the diffusion dt so every cell updates every step. With a larger intracellular_dt, cells that are not updating still draw on their persisted rates, which the sum does not include (colony, Ecoli_large_scale and cancer_metabolism configs are in that situation). Results (unit-test build, deterministic where stated): * voxel-sharing test, min glucose: n=1 -5e-6, n=4 -8e-6, n=8 -4e-6 mM (before: -1.8e-2, -0.23, -0.62); mixed voxel, 4 takers + 4 non-takers, draws exactly what 4 takers alone draw. * sweep S0 = 5e-4..1 mM, n = 1/3/8 (D=0, 1 thread): residual <= 1.0e-5 mM everywhere, against -(n-1) x content with the unit fix alone. * 4 OpenMP threads, n=8, x5: all pass. * mass-conversion gate, dense snapshots through the depletion window, n=6: worst excursion -3.8e-4 mM (two runs never negative) against -0.07..-0.19 mM unfixed. Cost: no measurable difference at 3000 cells. Changes simulated results wherever a substrate is limiting and cells share a voxel: manuscript-sensitive. Cases to rerun: bacterial_colony, cancer core metabolism, crossfeeding, and any acetate-switch run reaching depletion. docs/agent/ is not on this branch; architecture.md section 4 documents the old formula and must be updated when this is merged into feature/pfba. Co-Authored-By: Claude Sonnet 5.5 <[email protected]>
update setup_tissue to use cell.csv
Faster per-solve bookkeeping, header dependency tracking, custom_data corruption fix
Fix the substrate-uptake cap: a unit error that disabled it, and cells in one voxel overdrawing it
…om this branch fva_degeneracy.py, fva_figures.py, fva_degeneracy.ipynb and reviewer_response_figures.ipynb are analysis code for the reviewer response, not part of the pFBA implementation, and keeping them here made this branch's diff against development larger than the feature it carries. They are not lost: commit c4a1913 still adds them, and the first three are identical on feature/fva-degeneracy-diagnostic. reviewer_response_figures.ipynb is only in c4a1913, which stays reachable from the provenance tag. c4a1913 is cited as provenance by several result READMEs, so history is left as it is and the files are removed in a new commit instead.
dFBAModel relies on LP column i being reaction i in three places: initProblem() writes bounds and objective at the indexer's position but appends the matrix columns in vector order, and optimize() reads fluxes as columnPrimal[i] for reactions[i] (the per-solve bookkeeping change, b780351, made that last dependency explicit). addReaction() keeps the vector and the indexer in step, so nothing breaks today, but nothing checked it either. If anything ever reordered or removed reactions without rebuilding the indexer, fluxes would be assigned to the wrong reactions without any error. initProblem() now verifies that reactions[i] is at indexer position i and exits with the position and reaction id if not. The check is O(reactions) with a map lookup each, once per model construction, which is small against building the LP. This replaces the "reactionsIndexer is insertion order" line printed by the solve-path benchmark that was drafted as 328a6190 (never merged), which only printed yes/NO and could not fail a build. Verified: test_fba, test_sbml and test_clp give the same output as before; dfba_unit_test builds and the voxel-sharing suite passes 5/5; a one-cell mass-conversion run starts and completes; and a model whose reactions are swapped behind the indexer's back stops at initProblem with exit code 1 and the message above.
Integrate 92 upstream commits while retaining dFBA sources and sample projects. Adopt upstream intracellular scheduling, preserving secretion before intracellular updates. Keep the existing dFBA compiler flags and header dependency tracking in the two conflicting sample Makefiles.
…der-check Pr/development lp column order check
Pr/development pfba
…i-20261014 Fix Make CI and sync PhysiCell development
…ty sample The four glucose source patches of setup_secretion_nodes() belong to the multiple-sources example (300 x 300 um domain, centres at +-90 um), but main.cpp registered them for every config. The motility configs have a domain of y in [-50, 50], so BioFVM's nearest_voxel_index() truncated the centres onto the edge row and created 0.79 mM glucose hot spots above the 0.6 mM switch threshold. With the Km 0.8 mM config the cell therefore left the motility mode at 415 min and divided at 550 min, instead of staying motile for the whole run. The default PhysiCell_settings.xml, which has no glucose_source_* parameters, stopped at start-up with "Unknown parameter glucose_source_radius". - setup_secretion_nodes() is gated by the bool user parameter glucose_sources_enabled; a missing parameter counts as false. - When enabled, a warning names every source centre that lies outside the domain and the voxel it is truncated to (negative out-of-range indices are wrapped to the opposite edge by BioFVM, so both bottom centres end up on the top row). - PhysiCell_settings_multiple_sources.xml sets the parameter to true; the motility, motility_proliferation and no_motility configs set it to false. - PhysiCell_settings.xml is now a copy of PhysiCell_settings_motility_proliferation.xml, so the default behaviour of the sample is to move and then proliferate. Tested in a stage copy outside the repository: Km 0.8 stays motile and undivided for 600 min (x = 129.4 um, glucose maximum 0.55 mM); the default config moves, switches by about 300 min and reaches 16 cells; multiple_sources with the flag on is bit-identical to the previous build over 60 min (7 snapshots); with the flag off or absent no patch is registered and glucose stays at 0; the warning fires for the +-90 centres in the +-50 domain. Manuscript-sensitive: results of the Km 0.8 mM motility case change (no longer divides) and the Km 0.02 mM case switches about 25 min earlier. Re-run the motility figures (Fig. 6).
…ry 0.65 mM, Km 0.5 mM With the glucose diffusion coefficient of the shipped configs (46800 um^2/min) the time step of 0.01 min gives D dt / dx^2 = 18.7. The BioFVM solver then flattens the left end of the glucose gradient (0.13 mM at x = -140 um instead of 0.03 mM, although the left Dirichlet value is 0) and the cell trajectories change with it. A Km of 0.8 mM only gives a motile cell because of that distortion: with the field resolved (dt 0.001, D/dt as in the colony use case) a 0.8 mM cell cannot meet the maintenance ATP and never moves. - dt_diffusion and intracellular_dt: 0.01 -> 0.001 min in the motility and the motility_proliferation configs (mechanics 0.1 and phenotype 6 unchanged). - Motility-only config: glucose boundary at xmax 0.55 -> 0.65 mM (the same as the growth-recovery config, as the Fig. 6 caption says), glucose Km 0.8 -> 0.5 mM. - PhysiCell_settings.xml is again a copy of PhysiCell_settings_motility_proliferation.xml. - PhysiCell_settings_no_motility.xml is not changed. Results (stage copy, seed 60): the motility-only cell stays motile for the 600 min and never divides (about +122 um at 10 h, speed reaching 1.0); the growth-recovery cell switches at 290 min and divides at 355 min; FBA and pFBA give the same trajectories. Manuscript-sensitive: Km, boundary and dt of the Fig. 6 use case change (text, Supplementary S5.6, Tables S7/S8 and Fig. 6 have to be updated); the shipped results at dt 0.01 are not reproduced.
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
10 participants
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.
I've updated the parameters from the metabolic-drive mobility example and fiexed dirichlet nodes left from the multi sources examples.