From 46fc3f9a8bfaff2c66fd70d5c665a9266881e106 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 5 Nov 2026 15:41:14 +0000 Subject: [PATCH 1/2] Import SBML reaction networks as verified graphs Add a from-scratch reader for SBML Level 2/3 core (vh/graph/sbml.py): compartments, species, global and local parameters, function definitions, initial assignments, assignment rules and reactions with MathML kinetic laws. Unsupported constructs (events, rate and algebraic rules, delays, fast reactions, stoichiometryMath, conversion factors, variable compartments, required packages) are refused by name. Documents with a DOCTYPE are refused and MathML compiles to closures over a fixed operator table. Models simulate, solve for steady states on the conservation-reduced system and report reduced eigenvalues. vh/graph/stoichiometry.py does exact rational algebra: canonical conservation laws (RREF of the left nullspace) and positive-mass stoichiometric consistency (Gevorgyan et al. 2008) with exactly checked certificates. SBML graphs carry species, reaction, compartment and moiety nodes; the schema recomputes N from edges and checks stoichiometry, location, exact complete conservation and consistency. Validation (scripts/validate_sbml.py, auto-research report): three CC0 BioModels vendored with SHA-256 provenance agree with libRoadRunner, an independent SBML engine used only as a test referee, to <=6.9e-10 on trajectories, <=1.7e-12 on steady states and <=1.3e-9 on reduced eigenvalues; periods match to ~1e-11. Across BIOMD0000000001-060, 49 models are read and 11 refused by name; trajectories of all 49 agree with the referee (median 1.7e-9). The consistency invariant flags 8 models with exact certificates, including a stoichiometry used as a unit conversion in Chassagnole 2002, as its own notes describe. libroadrunner joins the dev extras as a test-only referee. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01CoFhtFtE7wJgvHKAzNB3xD --- README.md | 8 +- .../2026-11-05-E-sbml-import/report.md | 212 +++ .../2026-11-05-E-sbml-import/results.json | 1110 +++++++++++++ auto-research/README.md | 1 + data/sbml/BIOMD0000000009.xml | 1429 +++++++++++++++++ data/sbml/BIOMD0000000010.xml | 733 +++++++++ data/sbml/BIOMD0000000012.xml | 1052 ++++++++++++ data/sbml/sources.json | 50 + docs/graph_formalization.md | 55 +- docs/img/kholodenko_mapk_sbml.mmd | 62 + docs/roadmap.md | 3 + pyproject.toml | 4 +- scripts/validate_sbml.py | 251 +++ tests/test_sbml.py | 406 +++++ vh/graph/render.py | 8 + vh/graph/sbml.py | 1045 ++++++++++++ vh/graph/stoichiometry.py | 259 +++ 17 files changed, 6683 insertions(+), 5 deletions(-) create mode 100644 auto-research/2026-11-05-E-sbml-import/report.md create mode 100644 auto-research/2026-11-05-E-sbml-import/results.json create mode 100644 data/sbml/BIOMD0000000009.xml create mode 100644 data/sbml/BIOMD0000000010.xml create mode 100644 data/sbml/BIOMD0000000012.xml create mode 100644 data/sbml/sources.json create mode 100644 docs/img/kholodenko_mapk_sbml.mmd create mode 100644 scripts/validate_sbml.py create mode 100644 tests/test_sbml.py create mode 100644 vh/graph/sbml.py create mode 100644 vh/graph/stoichiometry.py diff --git a/README.md b/README.md index 836b655..26de5ea 100644 --- a/README.md +++ b/README.md @@ -66,7 +66,7 @@ faster than real time. | `vh.molecular` | From-scratch JAX engine: AMBER ff14SB force field, GB-OBC2 implicit water, autodiff forces, Langevin (BAOAB) and Verlet dynamics, L-BFGS minimization | Matches OpenMM (used only as a referee) term by term to ~1e-14 relative on real proteins; trajectories identical step by step; energy conserved ([`docs/validation_molecular.md`](docs/validation_molecular.md)) | | `vh.body` | From-scratch whole-body PBPK: 13 organs + blood, Rodgers & Rowland tissue partitioning, liver (well-stirred, IVIVE), kidney (filtration), oral absorption (stomach + 7-segment intestine) with gut-wall first-pass, virtual populations; solved exactly with matrix exponentials | Bottom-up midazolam predictions vs 100 clinical datasets (IV and oral) digitized by the Open Systems Pharmacology project | | `vh.bridge` | Unbound drug at the target organ → receptor occupancy over time; Kd ↔ binding free energy | Analytic round-trip; kinetic occupancy | -| `vh.graph` | **Formalizes the body as a typed graph** across scales (organs↔blood flow, atoms↔bonds); graph algorithms for reachability, pathways, cycles, flow balance, bottlenecks; renders and verifies subgraphs | Blood conservation (Kirchhoff) to 1e-16 on the real PBPK model; a pure graph BFS reproduces the force field's exclusion list exactly ([`docs/graph_formalization.md`](docs/graph_formalization.md)) | +| `vh.graph` | **Formalizes the body as a typed graph** across scales (organs↔blood flow, species↔reactions, atoms↔bonds); graph algorithms for reachability, pathways, cycles, flow balance, bottlenecks; renders and verifies subgraphs; imports SBML reaction networks with exact conservation and consistency invariants | Blood conservation (Kirchhoff) to 1e-16 on the real PBPK model; a pure graph BFS reproduces the force field's exclusion list exactly; SBML trajectories and steady states match libRoadRunner to ~1e-9 ([`docs/graph_formalization.md`](docs/graph_formalization.md)) | | `web/` | Real-time demo: atom-level protein dynamics and a virtual patient, live in the browser | JS engine checked against the Python engine to ≤1e-12 | The explicit-solvent engine (`vh.molecular.periodic`, `.constraints`, `.explicit_md`) adds rung 3 of the accuracy ladder: every water molecule, exact long-range electrostatics (Ewald / smooth PME), rigid-water SHAKE/RATTLE and a Monte-Carlo barostat — validated term-by-term against OpenMM Ewald (energy to 1e-13, PME to 1e-11, constraints to 1e-16). @@ -158,12 +158,13 @@ vh/molecular energy.py (force field), integrate.py, minimize.py, prep.py, io.p periodic.py (Ewald/PME explicit solvent), constraints.py (SHAKE/RATTLE), explicit_md.py vh/body physiology.py, partition.py, compound.py, pbpk.py vh/bridge molecular → body coupling (Kd ↔ ΔG, receptor occupancy) -vh/graph core.py (typed property graph), algorithms.py, physiology_graph.py, render.py +vh/graph core.py (typed property graph), algorithms.py, physiology_graph.py, render.py, + contracts.py (model interfaces), sbml.py (SBML import), stoichiometry.py (exact algebra) tests/ pytest suite scripts/ validation report + figure generators docs/ accuracy_ladder, graph_formalization, feasibility, validation reports, architecture, roadmap web/ real-time browser engine + demo -data/ structures, tissue composition, clinical observations +data/ structures, tissue composition, clinical observations, SBML models (CC0) ``` ## Not a medical device @@ -184,6 +185,7 @@ research tracks, and how to publish work under [`auto-research/`](auto-research/ - **Rohan Arun** ([@rohanarun](https://github.com/rohanarun)) — creator and maintainer; direction and design. - **Claude** (Anthropic) — implemented and validated the engines under the above direction, in pair with the maintainer. - **OpenAI Codex** — graph interface contracts, clearance reference validation, and conserved reaction/intervention experiments, under the direction of [@ahrazzle](https://github.com/ahrazzle). +- **Claude** (Anthropic), Foundation-0 fork — SBML import with exact conservation and positive-mass consistency invariants. New contributors (human or agent): add yourself here in your pull request. See [`auto-research.md`](auto-research.md) for how to help. diff --git a/auto-research/2026-11-05-E-sbml-import/report.md b/auto-research/2026-11-05-E-sbml-import/report.md new file mode 100644 index 0000000..4c4f010 --- /dev/null +++ b/auto-research/2026-11-05-E-sbml-import/report.md @@ -0,0 +1,212 @@ +# SBML reaction networks as verified graphs + +- **Track**: E — grow the graph formalization (SBML import; conservation and consistency invariants) +- **Author**: Claude (Anthropic), under the direction of the Foundation-0 owner +- **Date**: 2026-11-05 +- **Branch / PR**: `auto/E-sbml-import` +- **Code added**: `vh/graph/sbml.py`, `vh/graph/stoichiometry.py`, `vh/graph/render.py` (styles), + `scripts/validate_sbml.py`, `tests/test_sbml.py`, `data/sbml/` (three CC0 BioModels with provenance) + +## Question + +Track E's milestone: *a BioModels SBML model imports to a graph whose +steady-state / conservation invariants are verified, and at least one new +invariant catches a real inconsistency.* Can a from-scratch reader turn +curated reaction networks into executable models and typed graphs whose +structural claims are proven exactly rather than asserted, and whose dynamics +agree with an independent SBML engine? + +## Method + +**Reader.** `read_sbml` parses SBML Level 2 (v1–5) and Level 3 (v1–2) core with +the standard-library XML parser. It refuses documents containing a DOCTYPE or +entity declarations. MathML is compiled to closures over a fixed operator table; +nothing is passed to `eval`. Supported constructs: +- compartments of constant size, species, and global and local parameters; +- function definitions, inlined at compile time; +- initial assignments, and assignment rules in dependency order (loops are + rejected); +- reactions with exact rational stoichiometry, modifiers and kinetic laws. + +Every other construct raises `UnsupportedSBML` and names it: events, rate and +algebraic rules, constraints, delays, fast reactions, stoichiometryMath, +conversion factors, non-constant compartments and required packages. + +The reader implements the SBML semantics the kinetics depend on: +- a species symbol means concentration unless `hasOnlySubstanceUnits`; +- kinetic laws are in substance per time; +- boundary species are not changed by reactions. + +It also rejects documents that violate them, such as a constant, non-boundary +species consumed by a reaction. Units are not checked. + +**Dynamics.** `d amount / dt = N v(x)` is integrated with SciPy (LSODA) on the +concentration-scaled state. Steady states come from a hybrid Newton solve on the +conservation-reduced system. Stability uses the Jacobian restricted to the +stoichiometric class: central differences reduced through the exact link matrix. + +**Exact stoichiometry** (`vh.graph.stoichiometry`). All coefficients are +`Fraction`s. +- **Conservation laws** are the RREF of the left nullspace of `N`, scaled to + coprime integers. That basis is canonical for the space, so two + implementations can be compared row by row. +- **Positive-mass (stoichiometric) consistency** follows Gevorgyan, Poolman & + Fell (2008): every reacting species must admit a strictly positive mass that + every internal reaction conserves. A floating-point LP classifies the + species, and the classification is then proved with exact certificates: + - for a conserved species, a mass vector `m ≥ 0` with `mᵀN = 0` and `mᵢ > 0`; + - for an unconserved species, a reaction combination `v` with `N v ≥ 0` and + `(N v)ᵢ > 0`, which shows by Farkas' lemma that no positive mass exists. + + Rounding is removed exactly. Rows the LP leaves numerically at zero are + constrained to be exactly zero in their exact nullspace, and the LP is solved + again. A result is labelled `certified` only if every certificate checks + exactly. +- **Internal reactions.** Boundary and constant species keep their rows, + because holding a species fixed does not change its mass. A reaction counts as + internal only if each side has a species that is not held fixed. Reactions + that touch only fixed pools on one side, such as `A → EmptySet`, are + exchanges. + +**Graph.** `SBMLModel.to_graph()` builds a typed graph: +- nodes for `compartment`, `species`, `reaction` and `moiety`; +- `reactant`/`product` edges carrying exact stoichiometry strings, plus + `modifier`, `located_in` and `member` edges. + +`sbml_schema()` recomputes `N` from the edges alone and checks four invariants: +stoichiometry, located, conservation (exact, complete, with matching totals) +and stoichiometric consistency. + +**Inputs and provenance.** Three curated BioModels whose notes state CC0 are +vendored unmodified. `data/sbml/sources.json` records each one's URL, SHA-256, +size and PubMed IDs: +- BIOMD0000000009, Huang & Ferrell 1996 +- BIOMD0000000010, Kholodenko 2000 +- BIOMD0000000012, Elowitz & Leibler 2000 + +Survey models (BIOMD0000000001–060) are identified by SHA-256 in +`results.json` and are not vendored. + +## Validation — the part that matters + +The authorities are of two kinds. +- **Exact rational algebra** for the structural claims, together with + closed-form solutions for synthetic networks. +- **libRoadRunner 2.11.0** (libSBML + CVODE, C++) as an independent SBML + implementation for the dynamics. It was not used to produce any value here. + +Tolerances were `rtol = 1e-12` and `atol = 1e-16` in both engines for the +vendored models. The survey used `rtol = 1e-10` and `atol = 1e-14`. + +| Quantity | This work | Authority | Error | +|---|---|---|---| +| Conservation laws, Kholodenko 2000 | 3 laws: MKKK, MKK and MAPK totals = 100, 300, 300 | exact `l·N = 0`; rank(N) + laws = 8; libRoadRunner moiety space | residual **0** (exact); same row space | +| Conservation laws, Huang & Ferrell 1996 | 7 laws (E1, E2, KKK, KK, K, KPase, KKPase totals) | same | residual **0**; same row space | +| Conservation laws, repressilator | none | rank(N) = 6 | exact | +| Trajectories, 3 models, 101 points | — | libRoadRunner | ≤ **6.9e-10** × species range | +| Drift of conserved totals along trajectories | — | exact invariance | ≤ **1.1e-13** relative | +| Steady state, Huang & Ferrell 1996 (stable) | max Re λ = −0.0716 | libRoadRunner steady state and reduced eigenvalues | state **1.7e-12**; eigenvalues **8.2e-10** | +| Fixed point, Kholodenko 2000 (unstable) | max Re λ = +0.00158 | same | state **1.5e-15**; eigenvalues **1.3e-9** | +| Fixed point, repressilator (unstable) | max Re λ = +0.0144 | same | state **1.2e-16**; eigenvalues **4.5e-10** | +| Limit-cycle period, Kholodenko 2000 (MKKK) | 1299.7705 time units | libRoadRunner trajectory | **8.9e-12** relative | +| Limit-cycle period, repressilator (PX) | 125.35416 time units | same | **1.8e-12** relative | +| First-order decay, `V = 2.5`, both species conventions | — | `A₀ e^{−kt}` | ≤ 1e-8 relative | +| Reversible `A ⇌ B`: relaxation, steady state, eigenvalue | — | closed form | 1e-8 / 1e-12 / 1e-6 | +| Survey trajectories, 49 models, `[0, 100]`, `rtol` 1e-10 | — | libRoadRunner | median **1.7e-9**; 47 of 49 ≤ 1e-6; worst 1.1e-5, which falls to **2.9e-7** at `rtol` 1e-12 | + +No parameter was fitted. + +**Survey of BIOMD0000000001–060.** +- **Read:** 49 of 60. +- **Refused by name (11):** rate rules (5), events (3), delays (2) and + stoichiometryMath (1). +- **Referee check:** every model that was read was also simulated by + libRoadRunner over the same window. The two Borghans 1997 calcium + oscillators (BIOMD0000000044/045) agree least closely: 1.1e-5 and 1.3e-6 at + `rtol` 1e-10. The script reruns them automatically at `rtol` 1e-12, where + the gap falls to 2.9e-7 and 3.3e-8. In an exploratory check (not recorded in + `results.json`), our own LSODA and Radau runs differed from each other by a + similar amount. The shrinking gap points to integration error on sharp + calcium spikes, not a semantic difference. + +**The new invariant catches real inconsistencies.** Positive-mass consistency +fails, with exact certificates, for 8 of the 49 models that were read: +BIOMD0000000008, 015, 018, 035, 049, 050, 051 and 052. The model-level check +and the graph invariant agree on all 49. Three cases were inspected: + +- **Chassagnole 2002, *E. coli* carbon metabolism (051).** The + phosphotransferase reaction is written + `cglcex + 65 cpep → 65 cg6p + 65 cpyr`. The model's own notes say the + biomass-to-volume ratio "is taken into account by appropriately adjusting + the stoichiometries", because both compartments have size 1. One glucose + therefore formally yields 65 glucose-6-phosphate: a unit conversion hidden + in stoichiometry. That is harmless for this model in isolation, but it would + silently corrupt any mass balance or model composition built on it. + Neutralising the factor does not make the model consistent. A second + certificate runs through `vPGDH`, written + `6-phosphogluconate → ribulose 5-phosphate` with no CO₂. +- **Vilar 2002 oscillator (035).** Transcription and translation are template + reactions (`DA → DA + MA`, `MA → A + MA`), and `C → R` degrades the bound + activator inside the complex. These are deliberate abstractions with + precursors and degradation products omitted. +- **Gardner 1998 cell cycle (008).** `Z → C` and `Z → Y` degrade one partner + inside a complex. This is also a deliberate abstraction. + +A positive-mass failure is therefore a proven property of the network as +written, not proof of a modelling error. The invariant makes the abstraction +visible and machine-checkable. Judging whether it is acceptable needs the +model's context. + +An earlier version of the check dropped boundary species before testing mass +balance. That made binding to a clamped ligand look like creation from nothing: +Schoeberl 2002 (019) was flagged because EGF is a boundary species. The final +rule, described under Method, removes that artifact without hiding real +template synthesis. + +## Limits + +- **Construct coverage.** 11 of the 60 surveyed models use constructs outside + the subset: events, rate rules, delays and stoichiometryMath. They are + refused, not approximated. +- **Units are not checked.** Numbers are used as written, as an SBML simulator + would use them. A dimensionally inconsistent kinetic law is not detected. +- **Consistency is structural.** It ignores kinetic laws and modifiers, and it + cannot tell an intended abstraction from an error. +- **Steady states depend on the starting point.** The solve finds the fixed + point of the given stoichiometric class nearest the guess. Multistability is + not enumerated. +- **Speed.** Kinetics run as Python closures. That is adequate here (the + largest surveyed model, about 100 species, simulates in seconds), but it is + neither JAX-differentiable nor fast. Making it so is a natural next step. +- **What agreement with libRoadRunner shows.** It establishes that both engines + implement the same SBML semantics on these models. It does not validate the + biology of any model. +- **No link to the rest of the stack.** Nothing here touches `vh.bridge`, + occupancy, or patient data. + +## Reproduce + +```bash +pip install -e '.[dev]' # includes libroadrunner (referee only) +PYTHONPATH=. pytest -q tests/test_sbml.py +mkdir -p /tmp/biomodels && for i in $(seq -w 1 60); do id=BIOMD00000000$i; \ + curl -sSL -o /tmp/biomodels/$id.xml "https://www.ebi.ac.uk/biomodels/model/download/$id?filename=${id}_url.xml"; done +PYTHONPATH=. python scripts/validate_sbml.py --survey /tmp/biomodels # writes results.json here +``` + +The survey files' SHA-256 values are recorded in `results.json`. A later +BioModels release may serve different bytes; compare against those hashes. + +## References + +- Hucka M. et al. The Systems Biology Markup Language (SBML): language specification for Level 3 Version 2 Core (2019); Level 2 Version 4 (2008). +- Gevorgyan A., Poolman M. G., Fell D. A. Detection of stoichiometric inconsistencies in biomolecular models. *Bioinformatics* 24:2245–2251 (2008). +- Welsh C. et al. libRoadRunner 2.0: a high performance SBML simulation and analysis library. *Bioinformatics* 39:btac770 (2023). +- Le Novère N. et al. BioModels Database. *Nucleic Acids Res.* 34:D689 (2006); models used under CC0 as stated in each file. +- Models, with the PubMed IDs annotated in each SBML file: + - Kholodenko B. N. *Eur. J. Biochem.* 267:1583 (2000), PMID 10712587. + - Huang C.-Y. F., Ferrell J. E. *PNAS* 93:10078 (1996), PMID 8816754. + - Elowitz M. B., Leibler S. *Nature* 403:335 (2000), PMID 10659856. + - Vilar J. M. G. et al. *PNAS* 99:5988 (2002), PMID 11972055. + - Gardner T. S. et al. *PNAS* 95:14190 (1998), PMID 9826676. + - Chassagnole C. et al. *Biotechnol. Bioeng.* 79:53 (2002). BIOMD0000000051 annotates PMID 17590932; it was not checked against this citation. diff --git a/auto-research/2026-11-05-E-sbml-import/results.json b/auto-research/2026-11-05-E-sbml-import/results.json new file mode 100644 index 0000000..3bc65a9 --- /dev/null +++ b/auto-research/2026-11-05-E-sbml-import/results.json @@ -0,0 +1,1110 @@ +{ + "authority": "exact rational algebra; libRoadRunner as independent SBML referee", + "models": [ + { + "id": "BIOMD0000000009", + "name": "Huang1996 - Ultrasensitivity in MAPK cascade", + "sha256": "1f95d793517d47dbf98321a89cdb14fcff9f461427b4ea062d398734c962ab98", + "sbml": "L2V4", + "species": 26, + "dynamic_species": 22, + "reactions": 20, + "assignment_rules": 4, + "conservation_laws": [ + "E1 + E1_KKK", + "E2 + E2_P_KKK", + "KKK + P_KKK + E1_KKK + E2_P_KKK + P_KKK_KK + P_KKK_P_KK", + "KK + P_KK + PP_KK + P_KKK_KK + P_KKK_P_KK + PP_KK_K + PP_KK_P_K + KKPase_PP_KK + KKPase_P_KK", + "K + P_K + PP_K + PP_KK_K + PP_KK_P_K + KPase_PP_K + KPase_P_K", + "KPase + KPase_PP_K + KPase_P_K", + "KKPase + KKPase_PP_KK + KKPase_P_KK" + ], + "conserved_totals_initial_amount": [ + 1.2e-16, + 1.1999999999999998e-15, + 1.2e-14, + 4.8e-12, + 4.8e-12, + 4.8e-13, + 1.1999999999999998e-15 + ], + "conserved_totals_initial_concentration": [ + 3e-05, + 0.0003, + 0.003, + 1.2, + 1.2, + 0.12, + 0.0003 + ], + "exact_residual_l_dot_N": "0", + "rank_plus_laws_equals_dynamic_species": true, + "positive_mass_consistent": true, + "consistency_certified": true, + "graph_invariant_problems": [], + "referee": { + "engine": "libRoadRunner 2.11.0", + "trajectory_max_error_rel_to_species_max": 5.617227786242932e-10, + "conservation_space_equal": true + }, + "max_relative_drift_of_conserved_totals": 1.1107161545511846e-13, + "steady_state": { + "solver_success": true, + "max_abs_rate_of_change": 4.8467614016778965e-15, + "max_diff_vs_referee_rel_to_max": 1.7004532594088963e-12, + "eigenvalue_max_diff_vs_referee_rel": 8.216172658283124e-10, + "max_real_eigenvalue": -0.07161415016477386, + "stable": true + } + }, + { + "id": "BIOMD0000000010", + "name": "Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade", + "sha256": "906964e07ffd570ad94a03c1222e230b54b913cd536b008392bfb6155930a966", + "sbml": "L2V4", + "species": 8, + "dynamic_species": 8, + "reactions": 10, + "assignment_rules": 0, + "conservation_laws": [ + "MKKK + MKKK_P", + "MKK + MKK_P + MKK_PP", + "MAPK + MAPK_P + MAPK_PP" + ], + "conserved_totals_initial_amount": [ + 100.0, + 300.0, + 300.0 + ], + "conserved_totals_initial_concentration": [ + 100.0, + 300.0, + 300.0 + ], + "exact_residual_l_dot_N": "0", + "rank_plus_laws_equals_dynamic_species": true, + "positive_mass_consistent": true, + "consistency_certified": true, + "graph_invariant_problems": [], + "referee": { + "engine": "libRoadRunner 2.11.0", + "trajectory_max_error_rel_to_species_max": 6.83867922800367e-10, + "conservation_space_equal": true + }, + "max_relative_drift_of_conserved_totals": 4.547473508864641e-15, + "steady_state": { + "solver_success": true, + "max_abs_rate_of_change": 5.551115123125783e-17, + "max_diff_vs_referee_rel_to_max": 1.4870278679293009e-15, + "eigenvalue_max_diff_vs_referee_rel": 1.3041323793590769e-09, + "max_real_eigenvalue": 0.0015808867973344746, + "stable": false + }, + "oscillation": { + "species": "MKKK", + "period": 1299.7705094979783, + "referee_period": 1299.7705095095023, + "relative_difference": 8.866165148521614e-12, + "cycles_averaged": 2, + "model_time_units": "as written in the model" + } + }, + { + "id": "BIOMD0000000012", + "name": "Elowitz2000 - Repressilator", + "sha256": "d535c84bc5fe68b01728cac5c53f8fed971176e3f0f2367a9ba1e9c043b11a2d", + "sbml": "L2V3", + "species": 6, + "dynamic_species": 6, + "reactions": 12, + "assignment_rules": 9, + "conservation_laws": [], + "conserved_totals_initial_amount": [], + "conserved_totals_initial_concentration": [], + "exact_residual_l_dot_N": "0", + "rank_plus_laws_equals_dynamic_species": true, + "positive_mass_consistent": true, + "consistency_certified": true, + "graph_invariant_problems": [], + "referee": { + "engine": "libRoadRunner 2.11.0", + "trajectory_max_error_rel_to_species_max": 6.908546759502274e-10, + "conservation_space_equal": true + }, + "max_relative_drift_of_conserved_totals": 0.0, + "steady_state": { + "solver_success": true, + "max_abs_rate_of_change": 3.552713678800501e-15, + "max_diff_vs_referee_rel_to_max": 1.1801944309700211e-16, + "eigenvalue_max_diff_vs_referee_rel": 4.497717690631207e-10, + "max_real_eigenvalue": 0.0144474668713797, + "stable": false + }, + "oscillation": { + "species": "PX", + "period": 125.35415830519999, + "referee_period": 125.35415830543108, + "relative_difference": 1.8435520807816913e-12, + "cycles_averaged": 7, + "model_time_units": "as written in the model" + } + } + ], + "survey": { + "files": 60, + "read": 49, + "refused_by_construct": { + "model element ": 3, + "": 5, + "csymbol function": 2, + "stoichiometryMath": 1 + }, + "referee": { + "horizon": 100.0, + "rtol": 1e-10, + "atol": 1e-14, + "compared": 49, + "max_error_rel_to_species_range": 1.0786057563640103e-05, + "median_error": 1.7192407749872498e-09, + "within_1e-6": 47, + "rerun_at_tight_tolerance": { + "BIOMD0000000044": [ + 1.0786057563640103e-05, + 2.864557502320601e-07 + ], + "BIOMD0000000045": [ + 1.3126266067169491e-06, + 3.348195409229243e-08 + ] + }, + "tight_tolerance": { + "rtol": 1e-12, + "atol": 1e-16 + } + }, + "errors": [], + "inconsistent": [ + "BIOMD0000000008", + "BIOMD0000000015", + "BIOMD0000000018", + "BIOMD0000000035", + "BIOMD0000000049", + "BIOMD0000000050", + "BIOMD0000000051", + "BIOMD0000000052" + ], + "all_classifications_certified": true, + "graph_invariant_agrees_everywhere": true, + "models": [ + { + "id": "BIOMD0000000001", + "sha256": "407fccc883f80e6d120e34d8b0810a206a5cc6535a26d4574611504e683beb08", + "status": "refused", + "reason": "model element " + }, + { + "id": "BIOMD0000000002", + "sha256": "7f7ee585a886f9cd88199887a9e674e070ce368b2884ceeee65dd437f3edca7f", + "status": "read", + "name": "Edelstein1996 - EPSP ACh species", + "species": 13, + "reactions": 17, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.7031616501968764e-07 + }, + { + "id": "BIOMD0000000003", + "sha256": "a633872e4c602922e85ac487ef52ac8bd8ead762ec7c0f86f9e93032a84c90f5", + "status": "read", + "name": "Goldbeter1991 - Min Mit Oscil", + "species": 3, + "reactions": 7, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.962975820480737e-07 + }, + { + "id": "BIOMD0000000004", + "sha256": "504dc54b5697b4183bd0046fb11f0991702961e3dee0baf15ac1ecaad5009eb6", + "status": "read", + "name": "Goldbeter1991 - Min Mit Oscil, Expl Inact", + "species": 5, + "reactions": 7, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.5033265555272367e-07 + }, + { + "id": "BIOMD0000000005", + "sha256": "5a4f44b5a0bc8f32d3206f460f748236539de593be0fb016da2ed7cddd098e36", + "status": "read", + "name": "Tyson1991 - Cell Cycle 6 var", + "species": 9, + "reactions": 9, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.1896814023712015e-07 + }, + { + "id": "BIOMD0000000006", + "sha256": "c3cd0b2dfc738565fe72be3208eec28023e623297e198015df1df77a7e430820", + "status": "refused", + "reason": "" + }, + { + "id": "BIOMD0000000007", + "sha256": "9503b8fb38f0e285897a4d1155634e6893ae5b4ca843a3997098db8ee1614bda", + "status": "refused", + "reason": "model element " + }, + { + "id": "BIOMD0000000008", + "sha256": "6c49ab1dd2efddc682cc97563de927e000fbae1e471563bb2c62c37640d990f1", + "status": "read", + "name": "Gardner1998 - Cell Cycle Goldbeter", + "species": 5, + "reactions": 13, + "conservation_laws": 0, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "C", + "Y", + "Z" + ], + "leak_certificate": { + "reaction9": "3", + "reaction10": "-2", + "reaction11": "-2" + }, + "referee_trajectory_max_error": 4.1259901787439805e-09 + }, + { + "id": "BIOMD0000000009", + "sha256": "1f95d793517d47dbf98321a89cdb14fcff9f461427b4ea062d398734c962ab98", + "status": "read", + "name": "Huang1996 - Ultrasensitivity in MAPK cascade", + "species": 26, + "reactions": 20, + "conservation_laws": 7, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.0719422732093982e-09 + }, + { + "id": "BIOMD0000000010", + "sha256": "906964e07ffd570ad94a03c1222e230b54b913cd536b008392bfb6155930a966", + "status": "read", + "name": "Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade", + "species": 8, + "reactions": 10, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.4013873632606925e-09 + }, + { + "id": "BIOMD0000000011", + "sha256": "7d13449074b774b4ea328d939ce433efa1444b5224f8117066ab914cd0ff748d", + "status": "read", + "name": "Levchenko2000_MAPK_noScaffold", + "species": 22, + "reactions": 30, + "conservation_laws": 7, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 8.515842167035795e-10 + }, + { + "id": "BIOMD0000000012", + "sha256": "d535c84bc5fe68b01728cac5c53f8fed971176e3f0f2367a9ba1e9c043b11a2d", + "status": "read", + "name": "Elowitz2000 - Repressilator", + "species": 6, + "reactions": 12, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.7192407749872498e-09 + }, + { + "id": "BIOMD0000000013", + "sha256": "e4a9445a74eef88a9783c178f873eebccba72464e8da1fee5cc84e5b4b9192bc", + "status": "read", + "name": "Poolman2004_CalvinCycle", + "species": 27, + "reactions": 21, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.829230195370593e-13 + }, + { + "id": "BIOMD0000000014", + "sha256": "7646eb90f0dc8dd072a37612aa7c60fb516a8d4b7e2ab81ba5b5c8a315d23025", + "status": "read", + "name": "Levchenko2000_MAPK_Scaffold", + "species": 86, + "reactions": 300, + "conservation_laws": 8, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.878996358536559e-09 + }, + { + "id": "BIOMD0000000015", + "sha256": "74c5cc7d4ad82d45a8c4299b04ee8822c51aeaf19cb4ba4dfda09f55a323bb54", + "status": "read", + "name": "Curto1998 - purine metabolism", + "species": 18, + "reactions": 37, + "conservation_laws": 0, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "PRPP", + "IMP", + "SAMP", + "ATP", + "SAM", + "Ade", + "XMP", + "GTP", + "dATP", + "dGTP", + "RNA", + "DNA", + "HX", + "Xa", + "Gua", + "UA" + ], + "leak_certificate": { + "ampd": "-12", + "aprt": "-2", + "asuc": "1", + "dada": "-2", + "den": "15", + "dnaa": "-1", + "gdrnr": "1", + "gnuc": "-8", + "gprt": "-14", + "gua": "5", + "hxd": "-3", + "impd": "1", + "polyam": "-1", + "rnaa": "-9", + "rnag": "8", + "xd": "1" + }, + "referee_trajectory_max_error": 2.86241885614117e-10 + }, + { + "id": "BIOMD0000000016", + "sha256": "35ad538383991b45273e2a767c0325d996d188dba10ce959de9635564813f703", + "status": "read", + "name": "Goldbeter1995_CircClock", + "species": 7, + "reactions": 10, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 8.88362367422757e-09 + }, + { + "id": "BIOMD0000000017", + "sha256": "b75db9f5e0b418d7eb015e44565c992ebcda10f56f8934747a2677e222b16148", + "status": "read", + "name": "Hoefnagel2002_PyruvateBranches", + "species": 19, + "reactions": 14, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 7.577820458206507e-11 + }, + { + "id": "BIOMD0000000018", + "sha256": "79abf91af1c46fa23f4e1ca2cfedc8f879dcbf41bbd8fbc7a41046f6aa3c2a12", + "status": "read", + "name": "Morrison1989 - Folate Cycle", + "species": 33, + "reactions": 47, + "conservation_laws": 0, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "DHFRf", + "HCHO", + "FGAR", + "AICAR", + "dUMP", + "GAR", + "serine", + "formate", + "ATP", + "NADP", + "NADPH", + "homocysteine" + ], + "leak_certificate": { + "SHMT": "-1", + "SHMTr": "1", + "CH2FH4toHCHO": "1", + "MTHFR": "-1", + "MTR": "-1", + "HCOOHtoCHOFH4": "-1", + "GARFT": "-1", + "ATIC7": "-1", + "MTHFD": "-1", + "TYMS": "-1", + "DHFReductase": "-1", + "FFH2syn": "-2", + "ATIC12": "-2", + "AICARsyn": "-2", + "MTX2off": "1", + "MTX2deg": "-1" + }, + "referee_trajectory_max_error": 6.740855293246711e-10 + }, + { + "id": "BIOMD0000000019", + "sha256": "541c2a16196dfcd0649df14cff51fb80fc4157bad041c562b17679552baeff48", + "status": "read", + "name": "Schoeberl2002 - EGF MAPK", + "species": 100, + "reactions": 125, + "conservation_laws": 12, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.2250256132454956e-09 + }, + { + "id": "BIOMD0000000020", + "sha256": "5fbf09694a5fc3ec4ccec2bb8ba89909faeabb9ea709f44a555c994954f444cf", + "status": "refused", + "reason": "" + }, + { + "id": "BIOMD0000000021", + "sha256": "02376379c0270208b2fd72529a50ccefcdd68c8b23a40585bfdcdd0926a4c8da", + "status": "read", + "name": "Leloup1999_CircClock", + "species": 10, + "reactions": 24, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 5.464687585573705e-08 + }, + { + "id": "BIOMD0000000022", + "sha256": "0b569c675eaed0226dc4c1dcd7f7e2e210a7cbd1b5562da97e55f47c255cdf04", + "status": "read", + "name": "Ueda2001_CircClock", + "species": 13, + "reactions": 32, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.1625583922709105e-08 + }, + { + "id": "BIOMD0000000023", + "sha256": "dccd31fe46f4fb06484a9e3b95718aad788692e9abad1dc68fd368d7c8247066", + "status": "read", + "name": "Rohwer2001_Sucrose", + "species": 13, + "reactions": 11, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 3.7213379600053997e-09 + }, + { + "id": "BIOMD0000000024", + "sha256": "748514f402eea94c49346701a30f2b429c9c7d5de922cd450e5a155d5f69dbe1", + "status": "refused", + "reason": "csymbol function http://www.sbml.org/sbml/symbols/delay" + }, + { + "id": "BIOMD0000000025", + "sha256": "5e86b47a659f8e3b4183dd03bf18e68da3419a340f87da966108518be062c0fa", + "status": "refused", + "reason": "csymbol function http://www.sbml.org/sbml/symbols/delay" + }, + { + "id": "BIOMD0000000026", + "sha256": "d293eb11fdd65b3150999c68edaf690a9449ea9e070a8f61878b3a4648eff7c5", + "status": "read", + "name": "Markevich2004_MAPK_orderedElementary", + "species": 11, + "reactions": 10, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 3.6056230459042765e-10 + }, + { + "id": "BIOMD0000000027", + "sha256": "56d8acc83add5a262b6785b18f2e905e606b9072625e25ba1520306bc4d586bb", + "status": "read", + "name": "Markevich2004 - MAPK double phosphorylation, ordered Michaelis-Menton", + "species": 5, + "reactions": 4, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 2.1135597305880588e-10 + }, + { + "id": "BIOMD0000000028", + "sha256": "2ecdc543c6ecb323ae31d13f8fec1846905cc5b9567827f1da8b0efef30ee107", + "status": "read", + "name": "Markevich2004_MAPK_phosphoRandomElementary", + "species": 16, + "reactions": 17, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 5.629178830445155e-10 + }, + { + "id": "BIOMD0000000029", + "sha256": "c1eb117678b7d423b8e1d2fd7f493c865fb66559980256f3544938f5a431b530", + "status": "read", + "name": "Markevich2004_MAPK_phosphoRandomMM", + "species": 6, + "reactions": 7, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 5.105826749333507e-10 + }, + { + "id": "BIOMD0000000030", + "sha256": "decd8cbc5b3981608b1bca8f5df1d6f5af95d9acb2e52c7b8c0d0346851f7bfb", + "status": "read", + "name": "Markevich2004_MAPK_AllRandomElementary", + "species": 18, + "reactions": 20, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.65955886828595e-10 + }, + { + "id": "BIOMD0000000031", + "sha256": "56ddd4cc0a0cf51c55a27fa386e0ee670c798ea69151ad39fadf363d0c0297a1", + "status": "read", + "name": "Markevich2004_MAPK_orderedMM2kinases", + "species": 6, + "reactions": 4, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 7.35699221928327e-10 + }, + { + "id": "BIOMD0000000032", + "sha256": "45cec22aefca6b49e78fd3ad5105748559df3e1dcc3bf5ffd787fda8bc6e60df", + "status": "read", + "name": "Kofahl2004_PheromonePathway", + "species": 37, + "reactions": 47, + "conservation_laws": 11, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 7.272437013147094e-10 + }, + { + "id": "BIOMD0000000033", + "sha256": "a0337a2687498ce30dfabd47cc15c5c634fc6d0ad1a7829fb569a36fdef9456d", + "status": "read", + "name": "Brown2004 - NGF and EGF signaling", + "species": 32, + "reactions": 26, + "conservation_laws": 15, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.061791087745291e-10 + }, + { + "id": "BIOMD0000000034", + "sha256": "05d93958690ebe26ab031bfdd3a288d091b33ad50b347076a98635b28d7357fe", + "status": "refused", + "reason": "" + }, + { + "id": "BIOMD0000000035", + "sha256": "3303169fdd9a672fa25c26586822c9a1b91ac461744d1926d317f493e6981237", + "status": "read", + "name": "Vilar2002_Oscillator", + "species": 10, + "reactions": 16, + "conservation_laws": 2, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "A", + "C", + "MA", + "MR", + "R" + ], + "leak_certificate": { + "Reaction1": "1", + "Reaction7": "1", + "Reaction10": "2", + "Reaction13": "1", + "Reaction16": "2" + }, + "referee_trajectory_max_error": 1.578098475130181e-08 + }, + { + "id": "BIOMD0000000036", + "sha256": "c95e89456a74ebf7bbef76e648ac4e83f59a913b894b01a69387458ae77802f4", + "status": "read", + "name": "Tyson1999_CircClock", + "species": 3, + "reactions": 5, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.6564204500450833e-08 + }, + { + "id": "BIOMD0000000037", + "sha256": "2a8f7a68ac7efd45e09f42a1a5f5e30c9a47b9c039fb79ecf86ef7e3def17c07", + "status": "read", + "name": "Marwan2003 - Genetics, regulatory hierarchy between genes", + "species": 12, + "reactions": 12, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 3.028687969882321e-09 + }, + { + "id": "BIOMD0000000038", + "sha256": "33708cec1fad44db6958c3ac6255397ad79b45564a296d29732fd934608d787b", + "status": "read", + "name": "Rohwer2000_Phosphotransferase_System", + "species": 17, + "reactions": 10, + "conservation_laws": 4, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.5330480698328444e-15 + }, + { + "id": "BIOMD0000000039", + "sha256": "6b221222be03747b11454cc4f2a6218697d9c6367d3578a170c7cc5c34ad3b12", + "status": "read", + "name": "Marhl2000_CaOscillations", + "species": 5, + "reactions": 7, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 4.222560473459409e-07 + }, + { + "id": "BIOMD0000000040", + "sha256": "0bf2abb25229f3fe7979315139396f9c0bcdb06370ae9c835f7d3fc0e0c79c6a", + "status": "refused", + "reason": "stoichiometryMath" + }, + { + "id": "BIOMD0000000041", + "sha256": "f6bcb0469ae530aacb147c7d35205766d57a8d4a612fab9b148ba90659078ceb", + "status": "read", + "name": "Kongas2007 - Creatine Kinase in energy metabolic signaling in muscle", + "species": 10, + "reactions": 9, + "conservation_laws": 3, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.5172626674860024e-10 + }, + { + "id": "BIOMD0000000042", + "sha256": "3c7946ce1e76401eca4d1b878d435083bf12cf6b16eee6b10316f2aa935109f9", + "status": "read", + "name": "Nielsen1998_Glycolysis", + "species": 15, + "reactions": 25, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 2.1915107454441513e-08 + }, + { + "id": "BIOMD0000000043", + "sha256": "63881219e93fc8552192bb3ed0c55bf1ee3f6756cbdaf84e188d94422b31d9b8", + "status": "read", + "name": "Borghans1997 - Calcium Oscillation - Model 1", + "species": 5, + "reactions": 7, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 6.108555195363742e-07 + }, + { + "id": "BIOMD0000000044", + "sha256": "d441dc324396b20899427988aa72d3637a9c1edaa954bac7399e5e3113b6ce24", + "status": "read", + "name": "Borghans1997 - Calcium Oscillation - Model 2", + "species": 4, + "reactions": 8, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.0786057563640103e-05, + "referee_trajectory_max_error_tight_tolerance": 2.864557502320601e-07 + }, + { + "id": "BIOMD0000000045", + "sha256": "a4ede2fcb4cf09bd48ad53745515e9bc4cee570efbb11abf5a9063ffa3ef786d", + "status": "read", + "name": "Borghans1997 - Calcium Oscillation - Model 3", + "species": 4, + "reactions": 8, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.3126266067169491e-06, + "referee_trajectory_max_error_tight_tolerance": 3.348195409229243e-08 + }, + { + "id": "BIOMD0000000046", + "sha256": "dc1d91abdde172fb3ddfd74cc421f8e35df660d3fffe6a329a73edb15015c03b", + "status": "read", + "name": "Olsen2003_peroxidase", + "species": 16, + "reactions": 15, + "conservation_laws": 2, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 2.4014558026888917e-08 + }, + { + "id": "BIOMD0000000047", + "sha256": "0982b27ae47c7df89c3186029464fee8b660e73035cc70eb6369649b1bcbe676", + "status": "refused", + "reason": "" + }, + { + "id": "BIOMD0000000048", + "sha256": "3f958d9b64a07f8ff0b8b56431c10b005bf1cbc1de936527399951a0eea63c08", + "status": "read", + "name": "Kholodenko1999 - EGFR signaling", + "species": 23, + "reactions": 25, + "conservation_laws": 6, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 5.910940654481731e-10 + }, + { + "id": "BIOMD0000000049", + "sha256": "746aaefe22c346f7853387e7768b1ea1dca2077c00e5aece7f7f3a9b5551e92e", + "status": "read", + "name": "Sasagawa2005_MAPK", + "species": 99, + "reactions": 150, + "conservation_laws": 18, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "EGFR", + "L_EGFR", + "L_EGFR_dimer", + "L_dpEGFR", + "proteosome", + "EGF" + ], + "leak_certificate": { + "re2": "-1", + "re8": "-1", + "J5": "-2", + "J9": "1", + "J11": "-1", + "J19": "1", + "J20": "1", + "J24": "-1", + "J25": "-1", + "J27": "1", + "J33": "-4", + "J35": "-3", + "J37": "-3", + "J38": "-4", + "J40": "1", + "J43": "-1", + "J45": "-4", + "J119": "-4", + "J118": "-1", + "J124": "-1", + "J164": "-4" + }, + "referee_trajectory_max_error": 9.831557283058588e-10 + }, + { + "id": "BIOMD0000000050", + "sha256": "4846396a6621661b5faa775f7f7b48aefe36c3396bb2da41aeb96f3319d084ad", + "status": "read", + "name": "Martins2003_AmadoriDegradation", + "species": 14, + "reactions": 16, + "conservation_laws": 2, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "Cn", + "_3DG", + "FA", + "_1DG", + "AA", + "Man", + "Glu", + "MG", + "Fru" + ], + "leak_certificate": { + "v2": "1/8", + "v3": "-1/8", + "v5": "1/2", + "v6": "-7/8", + "v8": "3/4", + "v9": "-7/8", + "v12": "-1/8", + "v13": "-1/4", + "v15": "1", + "v16": "1/8" + }, + "referee_trajectory_max_error": 8.112872312822886e-10 + }, + { + "id": "BIOMD0000000051", + "sha256": "167e78583aa85b49c0e283543f5c80e826a16d23cc06bb34e60c6dfbeb67f3fc", + "status": "read", + "name": "Chassagnole2002_Carbon_Metabolism", + "species": 18, + "reactions": 48, + "conservation_laws": 0, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "cpep", + "cglcex", + "cg6p", + "cpyr", + "cf6p", + "cg1p", + "cpg", + "cfdp", + "csed7p", + "cgap", + "ce4p", + "cxyl5p", + "crib5p", + "cdhap", + "cpgp", + "cpg3", + "cpg2", + "cribu5p" + ], + "leak_certificate": { + "vPTS": "-1", + "vPGI": "399", + "vPGM": "1", + "vG6PDH": "-466", + "vPFK": "83", + "vTA": "-157", + "vTKA": "-156", + "vTKB": "-158", + "vALDO": "82", + "vGAPDH": "5", + "vTIS": "81", + "vPGK": "4", + "vrpGluMu": "3", + "vENO": "2", + "vPK": "66", + "vPGDH": "-467", + "vR5PI": "-155", + "vRu5P": "-313" + }, + "referee_trajectory_max_error": 1.2624584773255264e-08 + }, + { + "id": "BIOMD0000000052", + "sha256": "6abff8a7b414eaaa62165f20e48f2244dd597d16699a1209ae171ae9d12cb89f", + "status": "read", + "name": "Brands2002 - Monosaccharide-casein systems", + "species": 11, + "reactions": 11, + "conservation_laws": 3, + "positive_mass_consistent": false, + "certified": true, + "graph_invariant_agrees": true, + "unconserved": [ + "Glu", + "Fru", + "Formic_acid", + "Triose", + "Acetic_acid", + "Cn", + "C5" + ], + "leak_certificate": { + "_J3": "1/9", + "_J5": "5/9", + "_J6": "1", + "_J7": "-2/9", + "_J8": "-8/9", + "_J9": "2/3", + "_J10": "-2/3" + }, + "referee_trajectory_max_error": 8.70301582569469e-10 + }, + { + "id": "BIOMD0000000053", + "sha256": "8c0cf08f4398ccb1989cc1b8cf9eaff86325f06a32b453551612af59676216e5", + "status": "read", + "name": "Ferreira2003_CML_generation2", + "species": 6, + "reactions": 12, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 2.3879181016267972e-08 + }, + { + "id": "BIOMD0000000054", + "sha256": "22b92a2a8c71b2577515d1add05e2e6c7ad68c0ead197bbb08a8e2c09ad6c67e", + "status": "read", + "name": "Ataullahkhanov1996_Adenylate", + "species": 3, + "reactions": 5, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 1.4326567113495416e-09 + }, + { + "id": "BIOMD0000000055", + "sha256": "0ec0148324f94a17028bdf31b03a57b53343f5e290da377af1ca10a610aa542f", + "status": "read", + "name": "Locke2005 - Circadian Clock", + "species": 13, + "reactions": 32, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 6.472601473664629e-09 + }, + { + "id": "BIOMD0000000056", + "sha256": "fc955bffc47e55f1f6376b002d5d5ca10ce40e63bbb2959910665b647f937bd1", + "status": "refused", + "reason": "model element " + }, + { + "id": "BIOMD0000000057", + "sha256": "b5c3c40cca539b7c93943b065e6878a77315c945b69be60e8f70f343cefb5b5d", + "status": "read", + "name": "Sneyd2002_IP3_Receptor", + "species": 6, + "reactions": 5, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 8.52987199620894e-11 + }, + { + "id": "BIOMD0000000058", + "sha256": "51e36d55367f126099cbf4d9fcb34ec7c0bdb670ae9bb7da1e8586e763c068a0", + "status": "read", + "name": "Bindschadler2001_coupled_Ca_oscillators", + "species": 4, + "reactions": 11, + "conservation_laws": 0, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 3.0804830676931965e-08 + }, + { + "id": "BIOMD0000000059", + "sha256": "62b74f8818f8bd3596fbb1c4b4f3ff09427f27fc0557c0495c6666fbbe425580", + "status": "refused", + "reason": "" + }, + { + "id": "BIOMD0000000060", + "sha256": "0d4d9a4a48c8f493afe72d6e9fbffd7afbc1e6218cdfb89a2dea47a53588a9b1", + "status": "read", + "name": "Keizer1996_Ryanodine_receptor_adaptation", + "species": 4, + "reactions": 3, + "conservation_laws": 1, + "positive_mass_consistent": true, + "certified": true, + "graph_invariant_agrees": true, + "referee_trajectory_max_error": 3.3428152758005716e-10 + } + ] + } +} diff --git a/auto-research/README.md b/auto-research/README.md index 203c0b7..41d7577 100644 --- a/auto-research/README.md +++ b/auto-research/README.md @@ -26,3 +26,4 @@ auto-research/ | 2026-11-05 | E (partial) | [Explicit model contracts](2026-11-05-E-model-contracts/report.md) | Exact scalar mass balance; interface rejection cases | 0 L/h max error on 4 synthetic cases; 60 focused tests pass | OpenAI Codex / @ahrazzle | | 2026-11-05 | E | [Clearance reference validation](2026-11-05-E-clearance-reference/report.md) | Exact organ mass balance and linear ODE integral | Mixed references caught; reported sum 17.104405 L/h, dose/AUC 17.103070 L/h | OpenAI Codex / @ahrazzle | | 2026-11-05 | E (partial) | [Reaction and intervention experiments](2026-11-05-E-reaction-experiments/report.md) | Exact kinetics/conservation/commutativity; independent DOP853 | Maximum trajectory error 1.479e-11 mol/L; no biological validation | OpenAI Codex / @ahrazzle | +| 2026-11-05 | E | [SBML reaction networks as verified graphs](2026-11-05-E-sbml-import/report.md) | Exact rational algebra; libRoadRunner as independent SBML engine | 3 CC0 BioModels: exact conservation laws, trajectories ≤7e-10, steady states ≤2e-12; positive-mass invariant flags 8 of 49 surveyed models with exact certificates | Claude (Anthropic) | diff --git a/data/sbml/BIOMD0000000009.xml b/data/sbml/BIOMD0000000009.xml new file mode 100644 index 0000000..0c30b98 --- /dev/null +++ b/data/sbml/BIOMD0000000009.xml @@ -0,0 +1,1429 @@ + + + + + +
Huang1996 - Ultrasensitivity in MAPK cascade
+

The temporal sequence of kinase activation, from MAPKKK (activated RAF) to the final effector MAPK (activated ERK), is described here. It is observed from the model that there is an increase in sensitivity along the levels of the cascade, where the activity of MAPK reaches its maximal before MAPKKK.

+
+

This model is described in the article:

+
Ultrasensitivity in the mitogen-activated protein kinase cascade. +
+
Huang CY, Ferrell JE Jr
+
Proc. Natl. Acad. Sci. U.S.A. 1996:93(19):10078-83
+

Abstract:

+

The mitogen-activated protein kinase (MAPK) cascade is a highly conserved series of three protein kinases implicated in diverse biological processes. Here we demonstrate that the cascade arrangement has unexpected consequences for the dynamics of MAPK signaling. We solved the rate equations for the cascade numerically and found that MAPK is predicted to behave like a highly cooperative enzyme, even though it was not assumed that any of the enzymes in the cascade were regulated cooperatively. Measurements of MAPK activation in Xenopus oocyte extracts confirmed this prediction. The stimulus/response curve of the MAPK was found to be as steep as that of a cooperative enzyme with a Hill coefficient of 4-5, well in excess of that of the classical allosteric protein hemoglobin. The shape of the MAPK stimulus/ response curve may make the cascade particularly appropriate for mediating processes like mitogenesis, cell fate induction, and oocyte maturation, where a cell switches from one discrete state to another.

+
+
+

The species K_PP_norm, KKK_P_norm and KK_PP_norm are the relative concentrations of the active MAPK, MAPKK and MAPKKK, that is the double, or single resp. phophorylated forms divided by the total concentrations of each kinase. For MAPK additionally the also active MAPK divided by the maximal concentration of active MAPK is given by rel_K_PP_max. The parameter K_PP_norm_max, the maximal ratio of active MapK, has to be calculated for each change of parameters.

+
+

This model is hosted on BioModels Database + and identified by: BIOMD0000000009 + .

+

To cite BioModels Database, please use: BioModels Database: An enhanced, curated and annotated resource for published quantitative kinetic models + .

+
+

To the extent possible under law, all copyright and related or neighbouring rights to this encoded model have been dedicated to the public domain worldwide. Please refer to CC0 Public Domain Dedication + for more information.

+
+ + + + + + + + + + + Machne + Rainer + + raim@tbi.univie.ac.at + + University of Vienna + + + + + + 2005-02-10T23:39:30Z + + + 2015-06-02T12:02:44Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ the maximal value of K_PP_norm at a stimulation conc. of 10.

+ + + + + + + + + + + + PP_K + KPase_PP_K + + + + PP_K + P_K + K + PP_KK_K + PP_KK_P_K + KPase_PP_K + KPase_P_K + + + + + + + + + K_PP_norm + K_PP_norm_max + + + + + + + + + + PP_KK + PP_KK_K + PP_KK_P_K + KKPase_PP_KK + + + + PP_KK + P_KK + KK + PP_KK_K + PP_KK_P_K + P_KKK_KK + P_KKK_P_KK + KKPase_PP_KK + KKPase_P_KK + + + + + + + + + + + P_KKK + P_KKK_KK + P_KKK_P_KK + + + + KKK + P_KKK + P_KKK_KK + P_KKK_P_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a1 + E1 + KKK + + + + d1 + E1_KKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k2 + E1_KKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a2 + E2 + P_KKK + + + + d2 + E2_P_KKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k2 + E2_P_KKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a3 + KK + P_KKK + + + + d3 + P_KKK_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k3 + P_KKK_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a4 + P_KK + KKPase + + + + d4 + KKPase_P_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k4 + KKPase_P_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a5 + P_KK + P_KKK + + + + d5 + P_KKK_P_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k5 + P_KKK_P_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a6 + PP_KK + KKPase + + + + d6 + KKPase_PP_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k6 + KKPase_PP_KK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a7 + K + PP_KK + + + + d7 + PP_KK_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k7 + PP_KK_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a8 + P_K + KPase + + + + d8 + KPase_P_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k8 + KPase_P_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a9 + P_K + PP_KK + + + + d9 + PP_KK_P_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k9 + PP_KK_P_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + + + + + a10 + PP_K + KPase + + + + d10 + KPase_PP_K + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + compartment + k10 + KPase_PP_K + + + + + + + + + + \ No newline at end of file diff --git a/data/sbml/BIOMD0000000010.xml b/data/sbml/BIOMD0000000010.xml new file mode 100644 index 0000000..b50d6c5 --- /dev/null +++ b/data/sbml/BIOMD0000000010.xml @@ -0,0 +1,733 @@ + + + + + +
Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade
+
+

The combination of ultrasensitivity and negative feedback bring sustained oscillations in the mitogen-activated protein kinase cascades.

+
+
+

This model is described in the article:

+
+ Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. +
+
Kholodenko BN
+
Eur. J. Biochem. 2000; 267(6):1583-8
+

Abstract:

+
+

Functional organization of signal transduction into protein phosphorylation cascades, such as the mitogen-activated protein kinase (MAPK) cascades, greatly enhances the sensitivity of cellular targets to external stimuli. The sensitivity increases multiplicatively with the number of cascade levels, so that a tiny change in a stimulus results in a large change in the response, the phenomenon referred to as ultrasensitivity. In a variety of cell types, the MAPK cascades are imbedded in long feedback loops, positive or negative, depending on whether the terminal kinase stimulates or inhibits the activation of the initial level. Here we demonstrate that a negative feedback loop combined with intrinsic ultrasensitivity of the MAPK cascade can bring about sustained oscillations in MAPK phosphorylation. Based on recent kinetic data on the MAPK cascades, we predict that the period of oscillations can range from minutes to hours. The phosphorylation level can vary between the base level and almost 100% of the total protein. The oscillations of the phosphorylation cascades and slow protein diffusion in the cytoplasm can lead to intracellular waves of phospho-proteins.

+
+
+
+

This model is hosted on BioModels Database + and identified by: BIOMD0000000010 + .

+

To cite BioModels Database, please use: BioModels Database: An enhanced, curated and annotated resource for published quantitative kinetic models + .

+
+
+

To the extent possible under law, all copyright and related or neighbouring rights to this encoded model have been dedicated to the public domain worldwide. Please refer to CC0 Public Domain Dedication + for more information.

+
+ + + + + + + + + + Sauro + Herbert + + Herbert_Sauro@kgi.edu + + Keck Graduate Institute + + + + + + 2005-02-12T00:18:12Z + + + 2015-06-02T12:04:33Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V1 + MKKK + + + + + + 1 + + + + + MAPK_PP + Ki + + n + + + + + K1 + MKKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V2 + MKKK_P + + + + KK2 + MKKK_P + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + k3 + MKKK_P + MKK + + + + KK3 + MKK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + k4 + MKKK_P + MKK_P + + + + KK4 + MKK_P + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V5 + MKK_PP + + + + KK5 + MKK_PP + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V6 + MKK_P + + + + KK6 + MKK_P + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + k7 + MKK_PP + MAPK + + + + KK7 + MAPK + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + k8 + MKK_PP + MAPK_P + + + + KK8 + MAPK_P + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V9 + MAPK_PP + + + + KK9 + MAPK_PP + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + uVol + V10 + MAPK_P + + + + KK10 + MAPK_P + + + + + + + + + + + + \ No newline at end of file diff --git a/data/sbml/BIOMD0000000012.xml b/data/sbml/BIOMD0000000012.xml new file mode 100644 index 0000000..fe19368 --- /dev/null +++ b/data/sbml/BIOMD0000000012.xml @@ -0,0 +1,1052 @@ + + + + + +
Elowitz2000 - Repressilator
+
+

This model describes the deterministic version of the repressilator system.

+

The authors of this model (see reference) use three transcriptional repressor systems that are not part of any natural biological clock to build an oscillating network that they called the repressilator. The model system was induced in Escherichia coli.

+

In this system, LacI (variable X is the mRNA, variable PX is the protein) inhibits the tetracycline-resistance transposon tetR (Y, PY describe mRNA and protein). Protein tetR inhibits the gene Cl from phage Lambda (Z, PZ: mRNA, protein),and protein Cl inhibits lacI expression. With the appropriate parameter values this system oscillates.

+
+
+

This model is described in the article:

+
+ A synthetic oscillatory network of transcriptional regulators. +
+
Elowitz MB, Leibler S.
+
Nature. 2000 Jan; 403(6767):335-338
+

Abstract:

+
+

Networks of interacting biomolecules carry out many essential functions in living cells, but the 'design principles' underlying the functioning of such intracellular networks remain poorly understood, despite intensive efforts including quantitative analysis of relatively simple systems. Here we present a complementary approach to this problem: the design and construction of a synthetic network to implement a particular function. We used three transcriptional repressor systems that are not part of any natural biological clock to build an oscillating network, termed the repressilator, in Escherichia coli. The network periodically induces the synthesis of green fluorescent protein as a readout of its state in individual cells. The resulting oscillations, with typical periods of hours, are slower than the cell-division cycle, so the state of the oscillator has to be transmitted from generation to generation. This artificial clock displays noisy behaviour, possibly because of stochastic fluctuations of its components. Such 'rational network design may lead both to the engineering of new cellular behaviours and to an improved understanding of naturally occurring networks.

+
+
+
+

The model is based upon the equations in Box 1 of the paper; however, these equations as printed are dimensionless, and the correct dimensions have been returned to the equations, and the parameters set to reproduce Figure 1C (left).

+
+
+

The original model was generated by B.E. Shapiro using Cellerator version 1.0 update 2.1127 using Mathematica 4.2 for Mac OS X (June 4, 2002), November 27, 2002 12:15:32, using (PowerMac,PowerPC, Mac OS X,MacOSX,Darwin).

+

Nicolas Le Novere provided a corrected version generated by SBMLeditor on Sun Aug 20 00:44:05 BST 2006. This removed the EmptySet species. Ran fine on COPASI 4.0 build 18.

+

Bruce Shapiro revised the model with SBMLeditor on 23 October 2006 20:39 PST. This defines default units and correct reactions. The original Cellerator reactions while being mathematically correct did not accurately reflect the intent of the authors. The original notes were mostly removed because they were mostly incorrect in the revised version. Tested with MathSBML 2.6.0.

+

Nicolas Le Novere changed the volume to 1 cubic micrometre, to allow for stochastic simulation.

+

Changed by Lukas Endler to use the average livetime of mRNA instead of its halflife and a corrected value of alpha and alpha0.

+

Moreover, the equations used in this model were clarified, cf. below.

+

The equations given in box 1 + of the original publication are rescaled in three respects (lowercase letters denote the rescaled, uppercase letters the unscaled number of molecules per cell):

+
    +
  • the time is rescaled to the average mRNA lifetime, t_ave: τ = t/t_ave
    • +
  • the mRNA concentration is rescaled to the translation efficiency eff: m = M/eff
    • +
  • the protein concentration is rescaled to Km: p = P/Km
    • +
+

+ α + in the equations should be in units of rescaled proteins per promotor and cell, and β + is the ratio of the protein to the mRNA decay rates or the ratio of the mRNA to the protein halflife.

+

In this version of the model α + and β + are calculated correspondingly to the article, while p + and m + where just replaced by P/Km + resp. M/eff + and all equations multiplied by 1/t_ave + . Also, to make the equations easier to read, commonly used variables derived from the parameters given in the article by simple rules were introduced.

+

The parameters given in the article were:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
promotor strength (repressed) ( tps_repr + ): 5*10 -4transcripts/(promotor*s)
promotor strength (full) ( tps_active + ): 0.5transcripts/(promotor*s)
mRNA half life, τ 1/2,mRNA + : 2min
protein half life, τ 1/2,prot + : 10min
K M + : 40monomers/cell
Hill coefficient n:2 +
+

From these the following constants can be derived:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
average mRNA lifetime ( t_ave + ): + τ 1/2,mRNA + /ln(2) + = 2.89 min
mRNA decay rate ( kd_mRNA + ): + ln(2)/ τ 1/2,mRNA + = 0.347 min -1
protein decay rate ( kd_prot + ): + ln(2)/ τ 1/2,prot +
transcription rate ( a_tr + ): + tps_active*60 + = 29.97 transcripts/min
transcription rate (repressed) ( a0_tr + ): + tps_repr*60 + = 0.03 transcripts/min
translation rate ( k_tl + ): + eff*kd_mRNA + = 6.93 proteins/(mRNA*min)
α : + a_tr*eff*τ 1/2,prot + /(ln(2)*K M + ) + = 216.4 proteins/(promotor*cell*Km)
α 0 + : + a0_tr*eff*τ 1/2,prot + /(ln(2)*K M + ) + = 0.2164 proteins/(promotor*cell*Km)
β : + k_dp/k_dm + = 0.2
+
+

Annotation by the Kinetic Simulation Algorithm Ontology (KiSAO):

+

To reproduce the simulations run published by the authors, the model has to be simulated with any of two different approaches. First, one could use a deterministic method ( KISAO_0000035 + ) with continuous variables ( KISAO_0000018 + ). One sample algorithm to use is the CVODE solver ( KISAO_0000019 + ). Second, one could simulate the system using Gillespie's direct method ( KISAO_0000029 + ), which is a stochastic method ( KISAO_0000036 + ) supporting adaptive timesteps ( KISAO_0000041 + ) and using discrete variables ( KISAO_0000016 + ).

+
+
+

This model is hosted on BioModels Database + and identified by: BIOMD0000000012 + .

+

To cite BioModels Database, please use: BioModels Database: An enhanced, curated and annotated resource for published quantitative kinetic models + .

+
+
+

To the extent possible under law, all copyright and related or neighbouring rights to this encoded model have been dedicated to the public domain worldwide. Please refer to CC0 Public Domain Dedication + for more information.

+
+ + + + + + + + + + Le Novère + Nicolas + + lenov@ebi.ac.uk + + EMBL-EBI + + + + + + + + Shapiro + Bruce + + bshapiro@caltech.edu + + Jet Propulsion Laboratory + + + + + Juty + Nick + + juty@ebi.ac.uk + + EMBL-EBI + + + + + Endler + Lukas + + lukas@ebi.ac.uk + + EMBL-EBI + + + + + Chelliah + Vijayalakshmi + + viji@ebi.ac.uk + + EMBL-EBI + + + + + + 2009-01-20T14:03:56Z + + + 2013-07-10T10:59:30Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ lacI inhibitor

+ + + + + + + + + + + + + + + + +

+ Tet repressor protein

+ + + + + + + + + + + + + + + + +

+ lambda repressor

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ ratio of protein to mRNA decay rates

+ + + + +

+ Leakiness in protein copies per promoter and cell

+ + + + +

+ Protein copies per promoter and cell

+ + + + +

+ Average number of proteins per transcript

+ + + + +

+ Hill coefficient

+ + + + +

+ Number of repressor molecules per cell giving half maximal repression, in monomers per cell

+ + + + + + + +

+ mRNA decay rate constant

+ + + + +

+ Protein decay rate costant

+ + + + +

+ Translation rate constant

+ + + + +

+ Transcription rate from free promotor minus a0_tr

+ + + + +

+ Transcrition from free promotor in transcripts per second and promotor

+ + + + +

+ Transcrition from fully repressed promotor in transcripts per second and promotor

+ + + + +

+ Transcription rate from fully repressed promotor

+ + + + + + + + + tau_mRNA + + + 2 + + + + + + + + + tau_mRNA + tau_prot + + + + + + + + eff + t_ave + + + + + + + + + + ps_a + ps_0 + + 60 + + + + + + + + ps_0 + 60 + + + + + + + + + + 2 + + tau_prot + + + + + + + + + + 2 + + tau_mRNA + + + + + + + + + + a_tr + eff + tau_prot + + + + + + 2 + + KM + + + + + + + + + + + a0_tr + eff + tau_prot + + + + + + 2 + + KM + + + + + + + + + + + + + + + + + + + + + + + + + + + kd_mRNA + X + + + + + + + + + + + + + + + + + + + + + + + + + kd_mRNA + Y + + + + + + + + + + + + + + + + + + + + + + + + + kd_mRNA + Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + k_tl + X + + + + + + + + + + + + + + + + + + + + + + + + + + + + k_tl + Y + + + + + + + + + + + + + + + + + + + + + + + + + + + + k_tl + Z + + + + + + + + + + + + + + + + + + + + + + + + + kd_prot + PX + + + + + + + + + + + + + + + + + + + + + + + + + kd_prot + PY + + + + + + + + + + + + + + + + + + + + + + + + + kd_prot + PZ + + + + + + + + + + + + + + + + + + + + + + + + + + + + a0_tr + + + + + a_tr + + + KM + n + + + + + + + KM + n + + + + PZ + n + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + a0_tr + + + + + a_tr + + + KM + n + + + + + + + KM + n + + + + PX + n + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + a0_tr + + + + + a_tr + + + KM + n + + + + + + + KM + n + + + + PY + n + + + + + + + + + + \ No newline at end of file diff --git a/data/sbml/sources.json b/data/sbml/sources.json new file mode 100644 index 0000000..34e3fd7 --- /dev/null +++ b/data/sbml/sources.json @@ -0,0 +1,50 @@ +{ + "description": "Curated BioModels used as the SBML import regression set. Files are unmodified downloads.", + "retrieved": "2026-11-05", + "download_url_pattern": "https://www.ebi.ac.uk/biomodels/model/download/{id}?filename={id}_url.xml", + "models": [ + { + "id": "BIOMD0000000009", + "file": "BIOMD0000000009.xml", + "name": "Huang1996 - Ultrasensitivity in MAPK cascade", + "sbml": "L2V4", + "bytes": 54521, + "sha256": "1f95d793517d47dbf98321a89cdb14fcff9f461427b4ea062d398734c962ab98", + "url": "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000009?filename=BIOMD0000000009_url.xml", + "license": "CC0 1.0 (stated in the model's notes)", + "pubmed": [ + "20587024", + "8816754" + ] + }, + { + "id": "BIOMD0000000010", + "file": "BIOMD0000000010.xml", + "name": "Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade", + "sbml": "L2V4", + "bytes": 29295, + "sha256": "906964e07ffd570ad94a03c1222e230b54b913cd536b008392bfb6155930a966", + "url": "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000010?filename=BIOMD0000000010_url.xml", + "license": "CC0 1.0 (stated in the model's notes)", + "pubmed": [ + "10712587", + "20587024" + ] + }, + { + "id": "BIOMD0000000012", + "file": "BIOMD0000000012.xml", + "name": "Elowitz2000 - Repressilator", + "sbml": "L2V3", + "bytes": 43901, + "sha256": "d535c84bc5fe68b01728cac5c53f8fed971176e3f0f2367a9ba1e9c043b11a2d", + "url": "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000012?filename=BIOMD0000000012_url.xml", + "license": "CC0 1.0 (stated in the model's notes)", + "pubmed": [ + "10659856", + "20587024" + ] + } + ], + "note": "Only models whose own notes state the CC0 1.0 public-domain dedication are vendored. Survey-only models are identified in the auto-research report by id and SHA-256, not copied here." +} diff --git a/docs/graph_formalization.md b/docs/graph_formalization.md index 995fefa..efac275 100644 --- a/docs/graph_formalization.md +++ b/docs/graph_formalization.md @@ -86,6 +86,58 @@ flowchart LR The full diagrams remain reproducible from the linked Mermaid sources. Old raster/vector images with mixed reference values have been removed. +## Reaction networks from SBML + +`vh.graph.sbml` reads a declared subset of SBML (Levels 2 and 3 core) from +scratch: compartments, species, parameters, function definitions, initial +assignments, assignment rules and reactions with MathML kinetic laws. Anything +outside the subset — events, rate or algebraic rules, delays, fast reactions, +stoichiometryMath, conversion factors, required packages — is refused by name, +never silently dropped. A model becomes + +- an **executable network**: reaction rates, `d amount / dt = N v`, + trajectories, steady states and the Jacobian restricted to the stoichiometric + class; and +- a **typed graph**: `species`, `reaction`, `compartment` and `moiety` nodes; + `reactant`/`product` edges carrying exact stoichiometry, `modifier` edges, and + `member` edges from each species to the conserved moieties it belongs to. + +The graph schema adds four invariants, each recomputed from the graph's own +edges rather than trusted from the importer: + +- **stoichiometry** — every coefficient is a positive exact rational; +- **located** — every species is in exactly one compartment; +- **conservation** — each moiety is an exact left-null vector of `N` + (`l . N = 0` in rational arithmetic), the moieties are independent and span + the whole left nullspace, and each total equals the weighted initial amounts; +- **stoichiometric consistency** — every reacting species admits a strictly + positive mass conserved by all internal reactions (Gevorgyan et al. 2008). + Each answer carries an exact certificate: a mass vector, or a reaction + combination that provably creates the unconserved species from nothing. + +Validated on three CC0 BioModels in [`data/sbml/`](../data/sbml/sources.json) +against exact algebra and libRoadRunner as an independent referee (details in +the [auto-research report](../auto-research/2026-11-05-E-sbml-import/report.md)): + +| Model | Exact conservation laws | Referee agreement | +|---|---|---| +| Kholodenko 2000 MAPK oscillator (BIOMD0000000010) | 3 kinase totals (100, 300, 300) | trajectories 7e-10; unstable fixed point and eigenvalues match; period 1299.77 to 9e-12 | +| Huang & Ferrell 1996 MAPK cascade (BIOMD0000000009) | 7 enzyme/kinase/phosphatase totals | trajectories 6e-10; stable steady state 2e-12; eigenvalues 8e-10 | +| Elowitz & Leibler 2000 repressilator (BIOMD0000000012) | none (open system) | trajectories 7e-10; unstable fixed point; period 125.354 to 2e-12 | + +```mermaid +flowchart LR + MKKK -->|J0| MKKK_P -->|J1| MKKK + MKK -->|J2| MKK_P -->|J3| MKK_PP -->|J4| MKK_P -->|J5| MKK + MAPK -->|J6| MAPK_P -->|J7| MAPK_PP -->|J8| MAPK_P -->|J9| MAPK + MKKK_P -.->|activates J2, J3| MKK + MKK_PP -.->|activates J6, J7| MAPK + MAPK_PP -.->|inhibits J0| MKKK +``` + +The full imported graph, including reaction nodes and moieties, is in +[`img/kholodenko_mapk_sbml.mmd`](img/kholodenko_mapk_sbml.mmd). + ## The graph verifying the physics engine The sharpest demonstration that the formalization is faithful: a **pure graph @@ -130,5 +182,6 @@ physics a formal, inspectable, verifiable skeleton to hang on. ```bash python scripts/validate_graph.py # runs every algorithm above, renders the diagrams -pytest -q tests/test_graph.py tests/test_graph_clearance.py +python scripts/validate_sbml.py # SBML import against exact algebra and libRoadRunner +pytest -q tests/test_graph.py tests/test_graph_clearance.py tests/test_sbml.py ``` diff --git a/docs/img/kholodenko_mapk_sbml.mmd b/docs/img/kholodenko_mapk_sbml.mmd new file mode 100644 index 0000000..550465d --- /dev/null +++ b/docs/img/kholodenko_mapk_sbml.mmd @@ -0,0 +1,62 @@ +%% Kholodenko 2000 MAPK cascade imported from SBML (BIOMD0000000010) +flowchart LR + MKKK["Mos"] + MKKK_P["Mos-P"] + MKK["Mek1"] + MKK_P["Mek1-P"] + MKK_PP["Mek1-PP"] + MAPK["Erk2"] + MAPK_P["Erk2-P"] + MAPK_PP["Erk2-PP"] + J0["MAPKKK activation"] + J1["MAPKKK inactivation"] + J2["phosphorylation of MAPKK"] + J3["phosphorylation of MAPKK-P"] + J4["dephosphorylation of MAPKK-PP"] + J5["dephosphorylation of MAPKK-P"] + J6["phosphorylation of MAPK"] + J7["phosphorylation of MAPK-P"] + J8["dephosphorylation of MAPK-PP"] + J9["dephosphorylation of MAPK-P"] + moiety_1["MKKK + MKKK_P"] + moiety_2["MKK + MKK_P + MKK_PP"] + moiety_3["MAPK + MAPK_P + MAPK_PP"] + MKKK --> J0 + J0 --> MKKK_P + MAPK_PP -->|"modifies"| J0 + MKKK_P --> J1 + J1 --> MKKK + MKK --> J2 + J2 --> MKK_P + MKKK_P -->|"modifies"| J2 + MKK_P --> J3 + J3 --> MKK_PP + MKKK_P -->|"modifies"| J3 + MKK_PP --> J4 + J4 --> MKK_P + MKK_P --> J5 + J5 --> MKK + MAPK --> J6 + J6 --> MAPK_P + MKK_PP -->|"modifies"| J6 + MAPK_P --> J7 + J7 --> MAPK_PP + MKK_PP -->|"modifies"| J7 + MAPK_PP --> J8 + J8 --> MAPK_P + MAPK_P --> J9 + J9 --> MAPK + MKKK --> moiety_1 + MKKK_P --> moiety_1 + MKK --> moiety_2 + MKK_P --> moiety_2 + MKK_PP --> moiety_2 + MAPK --> moiety_3 + MAPK_P --> moiety_3 + MAPK_PP --> moiety_3 + classDef c_species fill:#2B44C4,stroke:#222,color:#ffffff; + class MKKK,MKKK_P,MKK,MKK_P,MKK_PP,MAPK,MAPK_P,MAPK_PP c_species; + classDef c_reaction fill:#586379,stroke:#222,color:#ffffff; + class J0,J1,J2,J3,J4,J5,J6,J7,J8,J9 c_reaction; + classDef c_moiety fill:#1A8A54,stroke:#222,color:#ffffff; + class moiety_1,moiety_2,moiety_3 c_moiety; diff --git a/docs/roadmap.md b/docs/roadmap.md index 71a02d8..6b78fe3 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -75,6 +75,9 @@ A cell is ~10¹⁴ atoms, so pure CCSD(T) will never reach it. The path up: from the molecular layer → action-potential prolongation → torsade risk category. *Done when* the 28 CiPA reference drugs are classified into their published risk categories. 10. **Cells and pathways.** SBML import (BioModels) so target occupancy drives signaling and cellular response. + *Started 2026-11-05:* a from-scratch SBML core reader with exact conservation and positive-mass + invariants, validated against libRoadRunner ([report](../auto-research/2026-11-05-E-sbml-import/report.md)). + Still open: events, rate rules, delays, units checking, and wiring `vh.bridge` occupancy into a pathway. 11. **Uncertainty end to end.** Propagate parameter uncertainty across scales; Bayesian calibration; because every layer is differentiable (JAX), gradients of a clinical endpoint with respect to molecular parameters. diff --git a/pyproject.toml b/pyproject.toml index cd26967..be4238e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,7 +28,9 @@ dependencies = [ prep = ["openmm>=8.1", "pdbfixer>=1.9"] # the quantum accuracy ceiling (rungs 4-5): ab initio electronic structure + geometry optimization quantum = ["pyscf>=2.5", "geometric>=1.0"] -dev = ["pytest>=8", "matplotlib>=3.8", "openmm>=8.1", "pdbfixer>=1.9", "pyscf>=2.5", "geometric>=1.0"] +# libroadrunner is the independent SBML referee for vh.graph.sbml (tests only) +dev = ["pytest>=8", "matplotlib>=3.8", "openmm>=8.1", "pdbfixer>=1.9", "pyscf>=2.5", "geometric>=1.0", + "libroadrunner>=2.7"] [tool.setuptools.packages.find] include = ["vh*"] diff --git a/scripts/validate_sbml.py b/scripts/validate_sbml.py new file mode 100644 index 0000000..033bd48 --- /dev/null +++ b/scripts/validate_sbml.py @@ -0,0 +1,251 @@ +"""SBML import validation: exact structure, an independent referee, and a survey. + +For each vendored BioModels file (data/sbml, CC0) this records the exact +conservation laws and their initial totals, positive-mass consistency with +exact certificates, and agreement with libRoadRunner (libSBML + CVODE) on +trajectories, the conservation space, steady states, reduced Jacobian +eigenvalues and, for oscillators, the limit-cycle period. + +With ``--survey DIR`` it also reads every BIOMD*.xml in DIR and records which +are refused (and why) and which fail positive-mass consistency, with the exact +certificate. Survey files are identified by SHA-256; they are not vendored. + + PYTHONPATH=. python scripts/validate_sbml.py [--survey DIR] +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +from pathlib import Path + +import numpy as np + +from vh.graph import Graph +from vh.graph import stoichiometry as S +from vh.graph.render import to_mermaid +from vh.graph.sbml import SBMLError, UnsupportedSBML, read_sbml + +ROOT = Path(__file__).resolve().parents[1] +DATA = ROOT / "data" / "sbml" +OUT = ROOT / "auto-research" / "2026-11-05-E-sbml-import" +IMG = ROOT / "docs" / "img" + +CASES = { # trajectory horizon, presimulation before the steady-state solve, oscillator? + "BIOMD0000000009": (200.0, 2000.0, False), + "BIOMD0000000010": (4000.0, 4000.0, True), + "BIOMD0000000012": (1000.0, 1000.0, True), +} + + +def _period(t, x): + """Mean interval between upward mean-crossings over the second half.""" + half = len(t) // 2 + tt, xx = t[half:], x[half:] - x[half:].mean() + up = np.nonzero((xx[:-1] < 0) & (xx[1:] >= 0))[0] + cross = tt[up] - xx[up] * (tt[up + 1] - tt[up]) / (xx[up + 1] - xx[up]) + gaps = np.diff(cross) + return float(gaps.mean()), float(gaps.std()), int(len(gaps)) + + +def _referee(path): + import roadrunner + + rr = roadrunner.RoadRunner(str(path)) + rr.integrator.relative_tolerance = 1e-12 + rr.integrator.absolute_tolerance = 1e-16 + rr.integrator.maximum_num_steps = 100000 + return rr, roadrunner.__version__ + + +def validate_model(mid: str) -> dict: + from scipy.optimize import linear_sum_assignment + + path = DATA / f"{mid}.xml" + m = read_sbml(path) + horizon, settle, oscillator = CASES[mid] + ids = m.dynamic_species + n = m.stoichiometry_matrix() + laws = m.conservation_laws() + x0 = dict(zip(ids, m.initial_amounts())) + cons = m.stoichiometric_consistency() + g, schema = m.to_graph() + row = { + "id": mid, "name": m.name, "sha256": m.source["sha256"], "sbml": f"L{m.level}V{m.version}", + "species": len(m.species), "dynamic_species": len(ids), "reactions": len(m.reactions), + "assignment_rules": len(m.rules), + "conservation_laws": [" + ".join(s if c == 1 else f"{c}*{s}" for s, c in law.items()) for law in laws], + "conserved_totals_initial_amount": [float(sum(c * x0[s] for s, c in law.items())) for law in laws], + "conserved_totals_initial_concentration": [ + float(sum(c * x0[s] for s, c in law.items()) / m.compartments[m.species[next(iter(law))].compartment].size) + if len({m.species[s].compartment for s in law}) == 1 else None for law in laws], + "exact_residual_l_dot_N": str(S.conservation_residual(S.conservation_laws(n), n)), + "rank_plus_laws_equals_dynamic_species": S.rank(n) + len(laws) == len(ids), + "positive_mass_consistent": cons["consistent"], "consistency_certified": cons["certified"], + "graph_invariant_problems": g.validate(schema), + } + + rr, version = _referee(path) + t = np.linspace(0, horizon, 101) + ref = np.array(rr.simulate(0, horizon, 101, selections=["time"] + ids)) + ours = m.simulate(ref[:, 0]) + traj = ours["amount"] + row["referee"] = {"engine": f"libRoadRunner {version}", + "trajectory_max_error_rel_to_species_max": float( + np.max(np.abs(traj - ref[:, 1:]) / np.abs(ref[:, 1:]).max(axis=0)))} + totals = np.array([[sum(c * traj[k, ids.index(s)] for s, c in law.items()) for law in laws] + for k in range(len(t))]) if laws else np.zeros((len(t), 0)) + row["max_relative_drift_of_conserved_totals"] = float( + np.max(np.abs(totals - totals[0]) / np.abs(totals[0]))) if laws else 0.0 + + rr2, _ = _referee(path) + rr2.conservedMoietyAnalysis = True + cm = rr2.getConservationMatrix() + arr = np.array(cm) + refm = arr[:, [list(cm.colnames).index(s) for s in ids]] if arr.size else np.zeros((0, len(ids))) + refm = refm[np.abs(refm).sum(axis=1) > 0] + oursm = np.array([[law.get(s, 0) for s in ids] for law in laws], float).reshape(-1, len(ids)) + rk = np.linalg.matrix_rank + row["referee"]["conservation_space_equal"] = bool( + rk(oursm) == rk(refm) == (rk(np.vstack([oursm, refm])) if len(oursm) else 0)) + + tail = m.simulate(np.linspace(0, settle, 401))["amount"][200:] + ss = m.steady_state(tail.mean(axis=0)) + rr2.simulate(0, settle, 3) + rr2.steadyState() + rss = np.array([rr2[s] for s in ids]) + ev = m.reduced_eigenvalues(ss["amount"]) + rev = np.array(rr2.getReducedEigenValues()) + cost = np.abs(ev[:, None] - rev[None, :]) + r_i, c_i = linear_sum_assignment(cost) + row["steady_state"] = { + "solver_success": ss["success"], "max_abs_rate_of_change": ss["max_abs_derivative"], + "max_diff_vs_referee_rel_to_max": float(np.max(np.abs(ss["amount"] - rss)) / np.max(np.abs(rss))), + "eigenvalue_max_diff_vs_referee_rel": float(np.max(cost[r_i, c_i]) / np.max(np.abs(rev))), + "max_real_eigenvalue": float(ev.real.max()), "stable": bool(ev.real.max() < 0), + } + if oscillator: + t_long = np.linspace(0, 1.5 * horizon if mid == "BIOMD0000000010" else 2 * horizon, 60001) + k = 0 + ours_p = _period(t_long, m.simulate(t_long)["amount"][:, k]) + rr3, _ = _referee(path) + ref_p = _period(t_long, np.array(rr3.simulate(0, t_long[-1], len(t_long), selections=["time", ids[k]]))[:, 1]) + row["oscillation"] = {"species": ids[k], "period": ours_p[0], "referee_period": ref_p[0], + "relative_difference": abs(ours_p[0] - ref_p[0]) / ref_p[0], + "cycles_averaged": ours_p[2], "model_time_units": "as written in the model"} + return row + + +SURVEY_HORIZON = 100.0 # model time units; the same window for every surveyed model + + +SURVEY_TOL = (1e-10, 1e-14) # (rtol, atol) in both engines +TIGHT_TOL = (1e-12, 1e-16) # re-run for any model above 1e-6, to separate integration error from semantics + + +def _survey_referee(path, m, tol=SURVEY_TOL): + """Worst |ours - libRoadRunner| / species range over [0, SURVEY_HORIZON], 51 points.""" + import roadrunner + + ids = m.dynamic_species + if not ids: + return None + roadrunner.Logger.setLevel(roadrunner.Logger.LOG_ERROR) + rr = roadrunner.RoadRunner(str(path)) + rr.integrator.relative_tolerance, rr.integrator.absolute_tolerance = tol + rr.integrator.maximum_num_steps = 500000 + ref = np.array(rr.simulate(0, SURVEY_HORIZON, 51, selections=["time"] + ids)) + ours = m.simulate(ref[:, 0], rtol=tol[0], atol=tol[1])["amount"] + scale = np.maximum(np.abs(ref[:, 1:]).max(axis=0), 1e-300) + return float(np.max(np.abs(ours - ref[:, 1:]) / scale)) + + +def survey(directory: Path) -> dict: + rows, refused = [], {} + for path in sorted(directory.glob("BIOMD*.xml")): + data = path.read_bytes() + entry = {"id": path.stem, "sha256": hashlib.sha256(data).hexdigest()} + try: + m = read_sbml(data) + except UnsupportedSBML as exc: + entry["status"] = "refused" + entry["reason"] = str(exc) + refused[str(exc).split(" '")[0].split(" http")[0]] = refused.get(str(exc).split(" '")[0].split(" http")[0], 0) + 1 + rows.append(entry) + continue + except SBMLError as exc: + entry["status"] = "error" + entry["reason"] = str(exc) + rows.append(entry) + continue + c = m.stoichiometric_consistency() + g, schema = m.to_graph() + graph_flag = any("[stoichiometric-consistency]" in p for p in g.validate(schema)) + entry.update({"status": "read", "name": m.name, "species": len(m.species), "reactions": len(m.reactions), + "conservation_laws": len(m.conservation_laws()), "positive_mass_consistent": c["consistent"], + "certified": c["certified"], "graph_invariant_agrees": graph_flag == (not c["consistent"])}) + if not c["consistent"]: + entry["unconserved"] = c["unconserved"] + entry["leak_certificate"] = c["leak_combination"] + entry["referee_trajectory_max_error"] = _survey_referee(path, m) + if entry["referee_trajectory_max_error"] is not None and entry["referee_trajectory_max_error"] > 1e-6: + entry["referee_trajectory_max_error_tight_tolerance"] = _survey_referee(path, m, TIGHT_TOL) + rows.append(entry) + read = [r for r in rows if r["status"] == "read"] + errs = [r["referee_trajectory_max_error"] for r in read if r["referee_trajectory_max_error"] is not None] + return {"files": len(rows), "read": len(read), "refused_by_construct": refused, + "referee": {"horizon": SURVEY_HORIZON, "rtol": SURVEY_TOL[0], "atol": SURVEY_TOL[1], + "compared": len(errs), "max_error_rel_to_species_range": max(errs), + "median_error": float(np.median(errs)), "within_1e-6": sum(e <= 1e-6 for e in errs), + "rerun_at_tight_tolerance": {r["id"]: [r["referee_trajectory_max_error"], + r["referee_trajectory_max_error_tight_tolerance"]] + for r in read if "referee_trajectory_max_error_tight_tolerance" in r}, + "tight_tolerance": {"rtol": TIGHT_TOL[0], "atol": TIGHT_TOL[1]}}, + "errors": [r for r in rows if r["status"] == "error"], + "inconsistent": [r["id"] for r in read if not r["positive_mass_consistent"]], + "all_classifications_certified": all(r["certified"] for r in read), + "graph_invariant_agrees_everywhere": all(r["graph_invariant_agrees"] for r in read), + "models": rows} + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("--survey", type=Path, help="directory of BIOMD*.xml files to survey") + args = ap.parse_args(argv) + results = {"authority": "exact rational algebra; libRoadRunner as independent SBML referee", + "models": [validate_model(mid) for mid in CASES]} + if args.survey: + results["survey"] = survey(args.survey) + OUT.mkdir(parents=True, exist_ok=True) + (OUT / "results.json").write_text(json.dumps(results, indent=1) + "\n") + + g, _ = read_sbml(DATA / "BIOMD0000000010.xml").to_graph() + view = Graph(g.name) # built in document order, so the output is reproducible + for node in g.nodes.values(): + if node.type in ("species", "reaction", "moiety"): + view.add_node(node.id, node.type, node.scale, node.label, **node.attrs) + for e in g.edges: + if e.type in ("reactant", "product", "modifier", "member"): + view.add_edge(e.src, e.dst, e.type, e.directed, e.model, **e.attrs) + (IMG / "kholodenko_mapk_sbml.mmd").write_text( + to_mermaid(view, "LR", title="Kholodenko 2000 MAPK cascade imported from SBML (BIOMD0000000010)") + "\n") + + for row in results["models"]: + print(f"{row['id']} {row['name']}: {len(row['conservation_laws'])} exact laws, " + f"consistent={row['positive_mass_consistent']} (certified {row['consistency_certified']}); " + f"trajectory err {row['referee']['trajectory_max_error_rel_to_species_max']:.1e}; " + f"steady state diff {row['steady_state']['max_diff_vs_referee_rel_to_max']:.1e}, " + f"eig diff {row['steady_state']['eigenvalue_max_diff_vs_referee_rel']:.1e}, " + f"stable={row['steady_state']['stable']}") + if args.survey: + sv = results["survey"] + print(f"survey: {sv['read']}/{sv['files']} read; refused {sv['refused_by_construct']}; " + f"inconsistent {sv['inconsistent']}; certified {sv['all_classifications_certified']}; " + f"referee {sv['referee']['within_1e-6']}/{sv['referee']['compared']} within 1e-6, " + f"worst {sv['referee']['max_error_rel_to_species_range']:.1e}") + print(f"wrote {OUT.relative_to(ROOT)}/results.json and docs/img/kholodenko_mapk_sbml.mmd") + + +if __name__ == "__main__": + main() diff --git a/tests/test_sbml.py b/tests/test_sbml.py new file mode 100644 index 0000000..069e173 --- /dev/null +++ b/tests/test_sbml.py @@ -0,0 +1,406 @@ +"""SBML import: analytic kinetics, exact stoichiometry, refusals, and a referee. + +Synthetic networks below are artificial fixtures with closed-form solutions. +The vendored BioModels (data/sbml, CC0) are checked against exact algebra and, +when installed, against libRoadRunner -- an independent SBML implementation +(libSBML + CVODE) used here only as a referee. +""" + +import hashlib +import json +import math +from fractions import Fraction +from pathlib import Path + +import numpy as np +import pytest + +from vh.graph import Graph +from vh.graph import stoichiometry as S +from vh.graph.sbml import SBMLError, UnsupportedSBML, read_sbml + +ROOT = Path(__file__).resolve().parents[1] +DATA = ROOT / "data" / "sbml" +M = 'xmlns="http://www.w3.org/1998/Math/MathML"' + + +def l3(body: str, version: int = 2, extra: str = "") -> bytes: + return (f'\n{body}' + ).encode() + + +def l2(body: str) -> bytes: + return ('\n{body}').encode() + + +def species(sid, comp="c", conc=None, amount=None, hosu="false", boundary="false", constant="false"): + init = f'initialConcentration="{conc}"' if conc is not None else f'initialAmount="{amount}"' + return (f'') + + +def reaction(rid, reactants, products, math, local="", reversible="false", modifiers=()): + def refs(tag, items): + if not items: + return "" + inner = "".join(f'' for s, c in items) + return f"<{tag}>{inner}" + mods = ("" + "".join(f'' for m in modifiers) + + "") if modifiers else "" + lp = f"{local}" if local else "" + return (f'{refs("listOfReactants", reactants)}' + f'{refs("listOfProducts", products)}{mods}{math}{lp}') + + +def apply(op, *args): + return f"<{op}/>{''.join(args)}" + + +def ci(x): + return f"{x}" + + +def cn(x): + return f"{x}" + + +COMP = '' + + +# --------------------------------------------------------------------------- +# kinetics against closed-form solutions +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("hosu", ["false", "true"]) +def test_first_order_decay_respects_concentration_versus_amount(hosu): + # V = 2.5, A0 = 10 (amount). Kinetic laws are substance/time: with + # hasOnlySubstanceUnits=false the symbol A is a concentration, so the law + # k*c*A and the amount-based law k*A describe the same amount decay. + size, k = 2.5, 0.3 + law = apply("times", ci("k"), ci("c"), ci("A")) if hosu == "false" else apply("times", ci("k"), ci("A")) + doc = l3(COMP.format(size=size) + + "" + species("A", amount=10, hosu=hosu) + species("B", amount=0, hosu=hosu) + + "" + + reaction("decay", [("A", 1)], [("B", 1)], law, f'') + + "") + m = read_sbml(doc) + t = np.linspace(0, 20, 41) + out = m.simulate(t) + exact = 10 * np.exp(-k * t) + np.testing.assert_allclose(out["amount"][:, 0], exact, rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(out["amount"].sum(axis=1), 10, rtol=1e-12) + expected_value = exact / size if hosu == "false" else exact + np.testing.assert_allclose(out["value"][:, 0], expected_value, rtol=1e-8, atol=1e-10) + assert m.conservation_laws() == [{"A": 1, "B": 1}] + + +def test_reversible_mass_action_relaxation_steady_state_and_eigenvalue(): + kf, kr, a0 = 2.0, 0.5, 3.0 + law = apply("minus", apply("times", ci("kf"), ci("A")), apply("times", ci("kr"), ci("B"))) + doc = l3(COMP.format(size=1) + '' + '' + + species("A", conc=a0) + species("B", conc=0) + "" + + reaction("r", [("A", 1)], [("B", 1)], law, reversible="true") + "") + m = read_sbml(doc) + a_eq = a0 * kr / (kf + kr) + t = np.linspace(0, 5, 26) + out = m.simulate(t) + np.testing.assert_allclose(out["amount"][:, 0], a_eq + (a0 - a_eq) * np.exp(-(kf + kr) * t), rtol=1e-8) + ss = m.steady_state() + assert ss["success"] + np.testing.assert_allclose(ss["amount"], [a_eq, a0 - a_eq], rtol=1e-12) + np.testing.assert_allclose(m.reduced_eigenvalues(ss["amount"]), [-(kf + kr)], rtol=1e-6) + + +def test_functions_piecewise_time_rules_initial_assignment_and_local_shadowing(): + body = ( + f'' + f'{ci("x")}{ci("K")}{apply("divide", ci("x"), apply("plus", ci("K"), ci("x")))}' + '' + + COMP.format(size=2) + + '' + '' + + "" + species("A", conc=1) + species("B", conc=0) + "" + + f'' + f'{apply("times", cn(3), ci("base"))}' + + f'' + f'{cn(2)}{apply("lt", "t", cn(5))}' + f'{cn(0)}' + + "" + + reaction("r", [("A", 1)], [("B", 1)], + apply("times", ci("c"), ci("k"), ci("g"), f"{ci('sat')}{ci('A')}{cn(0.5)}"), + '') + + "") + m = read_sbml(l3(body)) + x0 = m.initial_amounts() + np.testing.assert_allclose(x0, [3 * 1.5 * 2, 0]) # initial assignment is a concentration + a = 4.5 + expected = 2 * 0.25 * 2 * a / (0.5 + a) # local k shadows the global k = 100 + assert m.rates(x0, t=1.0)[0] == pytest.approx(expected, rel=1e-15) + assert m.rates(x0, t=6.0)[0] == 0.0 # piecewise on csymbol time + + +# --------------------------------------------------------------------------- +# refusals: unsupported constructs are named, never silently dropped +# --------------------------------------------------------------------------- +BASE = COMP.format(size=1) + "" + species("A", conc=1) + species("B", conc=0) + "" +DECAY = reaction("r", [("A", 1)], [("B", 1)], apply("times", cn(0.1), ci("A"))) + + +@pytest.mark.parametrize("doc,message", [ + (l3(BASE + "" + DECAY + ""), "listOfEvents"), + (l3(BASE + f'{cn(1)}'), "rateRule"), + (l3(BASE + f'{ci("A")}'), "algebraicRule"), + (l3(BASE + "" + reaction("r", [("A", 1)], [("B", 1)], + "d" + + ci("A") + cn(1) + "") + ""), "delay"), + (l3(BASE + "" + DECAY.replace('reversible="false"', 'reversible="false" fast="true"') + + "", version=1), "fast reaction"), + (l2('' + '' + '' + f'{cn(2)}' + f'{ci("A")}'), "stoichiometryMath"), + (l3(BASE, extra=' xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" fbc:required="true"'), + "required SBML package"), + (l3(BASE.replace('" + reaction("r", [("A", 1)], [("B", 1)], apply("factorial", ci("A"))) + + ""), "factorial"), + (l3(BASE + "" + reaction("r", [("A", 1)], [("B", 1)], ci("r")) + ""), + "reaction id"), + (l3(BASE.replace('constant="true"/>', 'constant="false"/>')), + "non-constant compartment"), + (l3(BASE + ""), "listOfConstraints"), +]) +def test_unsupported_constructs_are_refused_by_name(doc, message): + with pytest.raises(UnsupportedSBML, match=message): + read_sbml(doc) + + +@pytest.mark.parametrize("prefix", [ + b']>', + b'', +]) +def test_doctype_and_entity_declarations_are_refused(prefix): + with pytest.raises(SBMLError, match="DOCTYPE"): + read_sbml(prefix + l3(BASE).split(b"\n", 1)[1]) + + +@pytest.mark.parametrize("doc,message", [ + (l3(COMP.format(size=1) + "" + species("A", conc=1, constant="true") + species("B", conc=0) + + "" + DECAY + ""), "constant"), + (l3(BASE + '' + f'{ci("q")}' + f'{ci("p")}'), "algebraic loop"), + (l3(BASE + "" + reaction("r", [("A", 1)], [("B", 1)], ci("nope")) + ""), + "unknown symbols"), + (l3(COMP.format(size=1) + ''), "no initial value"), + (l3(COMP.format(size=1) + ''), + "needs attribute"), + (l3(BASE + "" + reaction("r", [("A", 0)], [("B", 1)], ci("A")) + ""), + "positive"), +]) +def test_semantic_errors_are_reported(doc, message): + with pytest.raises(SBMLError, match=message): + read_sbml(doc) + + +# --------------------------------------------------------------------------- +# exact stoichiometric algebra +# --------------------------------------------------------------------------- +def test_exact_rref_nullspace_and_conservation_laws(): + # A + E <-> AE -> B + E : conserved E + AE and A + AE + B. + n = [[-1, 1, 0], [-1, 1, 1], [1, -1, -1], [0, 0, 1]] # rows A, E, AE, B + laws = S.conservation_laws(n) + assert laws == [[1, 0, 1, 1], [0, 1, 1, 0]] + assert S.conservation_residual(laws, n) == 0 + assert S.rank(n) + len(laws) == 4 + for v in S.nullspace(n): + assert all(sum(Fraction(row[j]) * v[j] for j in range(3)) == 0 for row in n) + assert S.integer_row([Fraction(-1, 2), Fraction(1, 3)]) == [3, -2] + with pytest.raises(TypeError): + S.as_fraction(0.1) + + +def _check_certificates(n, res): + """Independent exact check of the returned certificates.""" + internal = res["internal_reactions"] + sub = [[Fraction(row[j]) for j in internal] for row in n] + mass = res["mass"] + assert all(m >= 0 for m in mass) and all(mass[i] > 0 for i in res["conserved"]) + assert all(sum(mass[i] * sub[i][k] for i in range(len(n))) == 0 for k in range(len(internal))) + if res["unconserved"]: + v = res["leak_combination"] + produced = [sum(sub[i][k] * v[k] for k in range(len(internal))) for i in range(len(n))] + assert all(p >= 0 for p in produced) and all(produced[i] > 0 for i in res["unconserved"]) + + +@pytest.mark.parametrize("n,internal,unconserved", [ + ([[-1, 0], [-1, 0], [1, -1], [0, 1]], None, []), # A + B -> C -> D + # A -> A + B: the net column hides the catalytic A, so the caller (as + # SBMLModel does from reactant/product lists) marks the reaction internal. + ([[0], [1]], [True], [1]), + ([[0], [1]], None, []), # same column read as pure synthesis + ([[-1], [-1], [1]], None, []), # A + B -> C + ([[-1, 0], [0, 1]], None, []), # A -> (out), (in) -> B: exchanges only + ([[-1, 1], [1, -2]], None, [0, 1]), # A -> B, 2B -> A: mass cannot balance + ([[-2, 1], [1, 0], [0, -1]], None, []), # 2A -> B, C -> A (fractional masses) +]) +def test_stoichiometric_consistency_with_exact_certificates(n, internal, unconserved): + res = S.stoichiometric_consistency(n, internal) + assert res["unconserved"] == unconserved + assert res["consistent"] == (not unconserved) + assert res["certified"] + _check_certificates(n, res) + + +# --------------------------------------------------------------------------- +# vendored BioModels: provenance, exact structure and graph invariants +# --------------------------------------------------------------------------- +SOURCES = json.loads((DATA / "sources.json").read_text()) +EXPECTED_LAWS = { + # Kinase totals of the three-tier cascade (Kholodenko 2000). + "BIOMD0000000010": [{"MKKK": 1, "MKKK_P": 1}, {"MKK": 1, "MKK_P": 1, "MKK_PP": 1}, + {"MAPK": 1, "MAPK_P": 1, "MAPK_PP": 1}], + # Two enzymes, three kinases and two phosphatases, each free plus bound (Huang & Ferrell 1996). + "BIOMD0000000009": [ + {"E1": 1, "E1_KKK": 1}, {"E2": 1, "E2_P_KKK": 1}, + {"KKK": 1, "P_KKK": 1, "E1_KKK": 1, "E2_P_KKK": 1, "P_KKK_KK": 1, "P_KKK_P_KK": 1}, + {"KK": 1, "P_KK": 1, "PP_KK": 1, "P_KKK_KK": 1, "P_KKK_P_KK": 1, "PP_KK_K": 1, "PP_KK_P_K": 1, + "KKPase_PP_KK": 1, "KKPase_P_KK": 1}, + {"K": 1, "P_K": 1, "PP_K": 1, "PP_KK_K": 1, "PP_KK_P_K": 1, "KPase_PP_K": 1, "KPase_P_K": 1}, + {"KPase": 1, "KPase_PP_K": 1, "KPase_P_K": 1}, {"KKPase": 1, "KKPase_PP_KK": 1, "KKPase_P_KK": 1}], + # Synthesis and degradation of every species: nothing is conserved (Elowitz & Leibler 2000). + "BIOMD0000000012": [], +} + + +def load(mid): + return read_sbml(DATA / f"{mid}.xml") + + +def test_vendored_models_match_recorded_provenance(): + assert {m["id"] for m in SOURCES["models"]} == set(EXPECTED_LAWS) + for entry in SOURCES["models"]: + data = (DATA / entry["file"]).read_bytes() + assert hashlib.sha256(data).hexdigest() == entry["sha256"] and len(data) == entry["bytes"] + assert entry["license"].startswith("CC0") and b"CC0" in data + + +@pytest.mark.parametrize("mid", sorted(EXPECTED_LAWS)) +def test_biomodels_exact_conservation_consistency_and_graph(mid): + m = load(mid) + assert m.conservation_laws() == EXPECTED_LAWS[mid] + n = m.stoichiometry_matrix() + laws = S.conservation_laws(n) + assert S.conservation_residual(laws, n) == 0 and S.rank(n) + len(laws) == len(m.dynamic_species) + res = m.stoichiometric_consistency() + assert res["consistent"] and res["certified"] + g, schema = m.to_graph() + assert g.validate(schema) == [] + restored = Graph.from_dict(json.loads(json.dumps(g.to_dict()))) + assert restored.validate(schema) == [] + + +def test_kholodenko_totals_are_conserved_along_the_trajectory(): + m = load("BIOMD0000000010") + out = m.simulate(np.linspace(0, 3000, 61)) + idx = {s: i for i, s in enumerate(out["species"])} + for law, total in zip(EXPECTED_LAWS["BIOMD0000000010"], (100, 300, 300)): + series = sum(c * out["amount"][:, idx[s]] for s, c in law.items()) + np.testing.assert_allclose(series, total, rtol=1e-9) + + +def test_graph_invariants_catch_tampered_structure(): + g, schema = load("BIOMD0000000009").to_graph() + member = next(e for e in g.edges_of_type("member") if e.src == "E1_KKK") + member.attrs["coefficient"] = "2" + assert any("[conservation] moiety_1: not conserved" in p for p in g.validate(schema)) + + g, schema = load("BIOMD0000000009").to_graph() + g = g.subgraph([n for n in g.nodes if n != "moiety_7"]) + assert any("6 moieties recorded" in p for p in g.validate(schema)) + + g, schema = load("BIOMD0000000010").to_graph() + next(iter(g.edges_of_type("product"))).attrs["stoichiometry"] = "0" + assert any("[stoichiometry]" in p for p in g.validate(schema)) + + # A product edge that creates MAPK from a conversion of MKKK breaks both + # the recorded moieties and positive-mass consistency. + g, schema = load("BIOMD0000000010").to_graph() + g.add_edge("J0", "MAPK", "product", stoichiometry="1") + problems = g.validate(schema) + assert any("[conservation]" in p for p in problems) + assert any("[stoichiometric-consistency]" in p and "exact certificate" in p for p in problems) + + +# --------------------------------------------------------------------------- +# independent referee: libRoadRunner +# --------------------------------------------------------------------------- +TIMES = {"BIOMD0000000010": 4000.0, "BIOMD0000000009": 200.0, "BIOMD0000000012": 1000.0} +# Presimulation before a steady-state solve: long enough to settle (stable +# model) or to sample a whole orbit (oscillators). Huang-Ferrell amounts are +# ~1e-13, so the referee's absolute steady-state threshold needs a start near +# the fixed point. +SETTLE = {"BIOMD0000000010": 4000.0, "BIOMD0000000009": 2000.0, "BIOMD0000000012": 1000.0} + + +def _referee(mid): + roadrunner = pytest.importorskip("roadrunner") + rr = roadrunner.RoadRunner(str(DATA / f"{mid}.xml")) + rr.integrator.relative_tolerance = 1e-12 + rr.integrator.absolute_tolerance = 1e-16 + rr.integrator.maximum_num_steps = 100000 + return rr + + +@pytest.mark.parametrize("mid", sorted(TIMES)) +def test_trajectories_match_libroadrunner(mid): + rr = _referee(mid) + m = load(mid) + ids = m.dynamic_species + ref = np.array(rr.simulate(0, TIMES[mid], 101, selections=["time"] + ids)) + ours = m.simulate(ref[:, 0]) + scale = np.abs(ref[:, 1:]).max(axis=0) + assert np.max(np.abs(ours["amount"] - ref[:, 1:]) / scale) < 1e-8 + + +@pytest.mark.parametrize("mid", sorted(TIMES)) +def test_conservation_space_matches_libroadrunner(mid): + rr = _referee(mid) + rr.conservedMoietyAnalysis = True + m = load(mid) + ids = m.dynamic_species + cm = rr.getConservationMatrix() + ref = np.array(cm)[:, [list(cm.colnames).index(s) for s in ids]] if np.array(cm).size else np.zeros((0, len(ids))) + ref = ref[np.abs(ref).sum(axis=1) > 0] # rows on rule-set species vanish here + ours = np.array([[law.get(s, 0) for s in ids] for law in m.conservation_laws()], float).reshape(-1, len(ids)) + rank = np.linalg.matrix_rank + assert rank(ours) == rank(ref) == (rank(np.vstack([ours, ref])) if len(ours) else 0) + + +@pytest.mark.parametrize("mid,stable", [("BIOMD0000000009", True), ("BIOMD0000000010", False), + ("BIOMD0000000012", False)]) +def test_steady_state_and_stability_match_libroadrunner(mid, stable): + from scipy.optimize import linear_sum_assignment + + rr = _referee(mid) + rr.conservedMoietyAnalysis = True + m = load(mid) + ids = m.dynamic_species + tail = m.simulate(np.linspace(0, SETTLE[mid], 401))["amount"][200:] + ss = m.steady_state(tail.mean(axis=0)) # an orbit's mean starts Newton near the fixed point + assert ss["success"] and ss["max_abs_derivative"] < 1e-9 + rr.simulate(0, SETTLE[mid], 3) + rr.steadyState() + ref = np.array([rr[s] for s in ids]) + assert np.max(np.abs(ss["amount"] - ref)) / np.max(np.abs(ref)) < 1e-7 + ours_ev = m.reduced_eigenvalues(ss["amount"]) + ref_ev = np.array(rr.getReducedEigenValues()) + cost = np.abs(ours_ev[:, None] - ref_ev[None, :]) + rows, cols = linear_sum_assignment(cost) + assert len(ours_ev) == len(ref_ev) + assert np.max(cost[rows, cols]) / np.max(np.abs(ref_ev)) < 1e-6 + assert (ours_ev.real.max() < 0) == stable diff --git a/vh/graph/render.py b/vh/graph/render.py index e8006e6..29adffe 100644 --- a/vh/graph/render.py +++ b/vh/graph/render.py @@ -22,6 +22,9 @@ "lumen": ("#B26A12", "#ffffff"), "target": ("#1A8A54", "#ffffff"), "atom": ("#2B44C4", "#ffffff"), + "species": ("#2B44C4", "#ffffff"), + "reaction": ("#586379", "#ffffff"), + "moiety": ("#1A8A54", "#ffffff"), } _EDGE_LABEL = { "flow": lambda e: f"{e.attrs.get('flow', 0):.1f} L/h", @@ -30,6 +33,11 @@ "absorbs": lambda e: f"Fa·Fg={e.attrs.get('fa',0)*e.attrs.get('fg',1):.2f}", "binds": lambda e: f"Kd {e.attrs.get('kd_nm','?')} nM", "bond": lambda e: "", + # SBML reaction networks: show a coefficient only when it is not 1 + "reactant": lambda e: "" if str(e.attrs.get("stoichiometry", "1")) == "1" else str(e.attrs["stoichiometry"]), + "product": lambda e: "" if str(e.attrs.get("stoichiometry", "1")) == "1" else str(e.attrs["stoichiometry"]), + "modifier": lambda e: "modifies", + "member": lambda e: "" if str(e.attrs.get("coefficient", "1")) == "1" else str(e.attrs["coefficient"]), } diff --git a/vh/graph/sbml.py b/vh/graph/sbml.py new file mode 100644 index 0000000..6968de4 --- /dev/null +++ b/vh/graph/sbml.py @@ -0,0 +1,1045 @@ +"""Read a declared subset of SBML into an executable reaction network and a typed graph. + +Written from scratch on the standard-library XML parser; libSBML/libRoadRunner +are used only as independent referees in the tests. Supported: SBML Level 2 +(versions 1-5) and Level 3 (versions 1-2) core models built from + + compartments (constant size), species, global and local parameters, + function definitions, initial assignments, assignment rules, and reactions + with reactants, products, modifiers and a MathML kinetic law. + +Anything outside that subset -- events, rate and algebraic rules, constraints, +delays, fast reactions, stoichiometryMath, conversion factors, variable +compartments, required packages, or a MathML element not listed in +``_OPERATORS`` -- raises ``UnsupportedSBML`` naming the construct. Nothing is +silently dropped except notes, annotations, units and SBO terms. Units are not +checked: the numbers are used as written, as an SBML simulator would. + +SBML semantics implemented here: a species symbol in math means its +concentration (amount / compartment size) unless ``hasOnlySubstanceUnits``; a +kinetic law is in substance per time, so ``d amount_i / dt = sum_j N_ij v_j``; +boundary-condition species are not changed by reactions. + +Documents containing a DOCTYPE are refused (no entity expansion); MathML is +compiled to closures over a fixed operator table, never to ``eval``. +""" + +from __future__ import annotations + +import hashlib +import math +import xml.etree.ElementTree as ET +from dataclasses import dataclass, field +from fractions import Fraction +from pathlib import Path + +import numpy as np + +from . import stoichiometry as S +from .core import EdgeRule, Graph, Schema + +MATHML = "http://www.w3.org/1998/Math/MathML" +_TIME = "http://www.sbml.org/sbml/symbols/time" +_AVOGADRO = "http://www.sbml.org/sbml/symbols/avogadro" +_DELAY = "http://www.sbml.org/sbml/symbols/delay" + + +class SBMLError(ValueError): + """The document is malformed or violates SBML semantics this reader relies on.""" + + +class UnsupportedSBML(SBMLError): + """The document is valid SBML but uses a construct outside the supported subset.""" + + +# --------------------------------------------------------------------------- +# MathML -> expression tree +# --------------------------------------------------------------------------- +# A node is a tuple: ("num", float) | ("sym", id) | ("time",) | ("op", name, args) +# | ("piecewise", ((value, cond), ...), otherwise_or_None) | ("call", fname, args) + +def _local(tag: str) -> str: + return tag.rsplit("}", 1)[-1] + + +def _ieee_div(a, b): + if b == 0: + if a == 0 or math.isnan(a): + return math.nan + return math.copysign(math.inf, a) * math.copysign(1.0, b) + return a / b + + +def _ieee_pow(a, b): + try: + r = a ** b + except ZeroDivisionError: + return math.inf + except OverflowError: + return math.inf + return math.nan if isinstance(r, complex) else r + + +def _ieee(fn): + """Domain errors become NaN and overflow becomes inf, as in IEEE arithmetic, + so a kinetic law evaluated off its domain does not abort an integration.""" + def wrapped(args): + try: + return fn(args) + except (ValueError, ZeroDivisionError): + return math.nan + except OverflowError: + return math.inf + return wrapped + + +def _log(args): + return math.log10(args[0]) if len(args) == 1 else math.log(args[1], args[0]) + + +def _chain(cmp): + return lambda args: all(cmp(a, b) for a, b in zip(args, args[1:])) + + +_OPERATORS = { + "plus": lambda a: sum(a) if a else 0.0, + "minus": lambda a: -a[0] if len(a) == 1 else a[0] - a[1], + "times": lambda a: math.prod(a) if a else 1.0, + "divide": lambda a: _ieee_div(a[0], a[1]), + "power": lambda a: _ieee_pow(a[0], a[1]), + "root": _ieee(lambda a: _ieee_pow(a[1], 1.0 / a[0]) if len(a) == 2 else math.sqrt(a[0])), + "exp": _ieee(lambda a: math.exp(a[0])), + "ln": _ieee(lambda a: math.log(a[0])), + "log": _ieee(_log), + "abs": lambda a: abs(a[0]), + "floor": lambda a: float(math.floor(a[0])), + "ceiling": lambda a: float(math.ceil(a[0])), + "sin": lambda a: math.sin(a[0]), "cos": lambda a: math.cos(a[0]), + "tan": lambda a: math.tan(a[0]), "arcsin": _ieee(lambda a: math.asin(a[0])), + "arccos": _ieee(lambda a: math.acos(a[0])), "arctan": lambda a: math.atan(a[0]), + "sinh": _ieee(lambda a: math.sinh(a[0])), "cosh": _ieee(lambda a: math.cosh(a[0])), + "tanh": lambda a: math.tanh(a[0]), + "eq": _chain(lambda x, y: x == y), "neq": lambda a: a[0] != a[1], + "gt": _chain(lambda x, y: x > y), "lt": _chain(lambda x, y: x < y), + "geq": _chain(lambda x, y: x >= y), "leq": _chain(lambda x, y: x <= y), + "and": lambda a: all(a), "or": lambda a: any(a), "xor": lambda a: sum(map(bool, a)) % 2 == 1, + "not": lambda a: not a[0], +} +_ARITY = {"minus": (1, 2), "divide": (2, 2), "power": (2, 2), "neq": (2, 2), "not": (1, 1)} +_UNARY = {"exp", "ln", "abs", "floor", "ceiling", "sin", "cos", "tan", "arcsin", "arccos", + "arctan", "sinh", "cosh", "tanh"} +_CONSTANTS = {"pi": math.pi, "exponentiale": math.e, "true": 1.0, "false": 0.0, + "infinity": math.inf, "notanumber": math.nan} +_INFIX = {"plus": "+", "minus": "-", "times": "*", "divide": "/", "power": "^"} + + +def _cn(el) -> tuple: + kind = el.attrib.get("type", "real") + parts = [p.strip() for p in "".join(el.itertext()).split()] if kind not in ("e-notation", "rational") else None + if kind in ("e-notation", "rational"): + seps = [c for c in el if _local(c.tag) == "sep"] + if len(seps) != 1: + raise SBMLError(f" needs exactly one ") + a = (el.text or "").strip() + b = (seps[0].tail or "").strip() + return ("num", float(a) * 10.0 ** float(b) if kind == "e-notation" else float(Fraction(a) / Fraction(b))) + if kind not in ("real", "integer", "double"): + raise UnsupportedSBML(f"") + if len(parts) != 1: + raise SBMLError(" must contain one number") + return ("num", float(parts[0])) + + +def parse_mathml(el) -> tuple: + """Expression tree of one MathML element (```` or a child).""" + tag = _local(el.tag) + if tag in ("math", "semantics"): + kids = [c for c in el if _local(c.tag) not in ("annotation", "annotation-xml")] + if len(kids) != 1: + raise SBMLError(f"<{tag}> must contain exactly one expression") + return parse_mathml(kids[0]) + if tag == "cn": + return _cn(el) + if tag == "ci": + return ("sym", "".join(el.itertext()).strip()) + if tag == "csymbol": + url = el.attrib.get("definitionURL", "").strip() + if url == _TIME: + return ("time",) + if url == _AVOGADRO: + return ("num", 6.02214076e23) + raise UnsupportedSBML(f"csymbol {url or '(no definitionURL)'}") + if tag in _CONSTANTS: + return ("num", _CONSTANTS[tag]) + if tag == "piecewise": + pieces, otherwise = [], None + for c in el: + kids = list(c) + if _local(c.tag) == "piece" and len(kids) == 2: + pieces.append((parse_mathml(kids[0]), parse_mathml(kids[1]))) + elif _local(c.tag) == "otherwise" and len(kids) == 1: + otherwise = parse_mathml(kids[0]) + else: + raise SBMLError("malformed ") + return ("piecewise", tuple(pieces), otherwise) + if tag == "apply": + kids = list(el) + if not kids: + raise SBMLError("empty ") + head, rest = kids[0], kids[1:] + h = _local(head.tag) + if h == "ci": + return ("call", "".join(head.itertext()).strip(), tuple(parse_mathml(k) for k in rest)) + if h == "csymbol": + raise UnsupportedSBML(f"csymbol function {head.attrib.get('definitionURL', '')}") + if h not in _OPERATORS: + raise UnsupportedSBML(f"MathML operator <{h}>") + args = [] + for k in rest: + q = _local(k.tag) + if q in ("degree", "logbase"): + inner = list(k) + if len(inner) != 1: + raise SBMLError(f"<{q}> must contain one expression") + args.insert(0, parse_mathml(inner[0])) + else: + args.append(parse_mathml(k)) + lo, hi = _ARITY.get(h, (1, 1) if h in _UNARY else (0, None)) + if h in ("root", "log"): + lo, hi = 1, 2 + if len(args) < lo or (hi is not None and len(args) > hi): + raise SBMLError(f"<{h}> applied to {len(args)} arguments") + return ("op", h, tuple(args)) + raise UnsupportedSBML(f"MathML element <{tag}>") + + +def symbols(node) -> set[str]: + kind = node[0] + if kind == "sym": + return {node[1]} + if kind == "op" or kind == "call": + out = set() + for a in node[2]: + out |= symbols(a) + return out + if kind == "piecewise": + out = set() + for v, c in node[1]: + out |= symbols(v) | symbols(c) + if node[2] is not None: + out |= symbols(node[2]) + return out + return set() + + +def infix(node) -> str: + """Human-readable form for graph attributes and reports (not re-parsed).""" + kind = node[0] + if kind == "num": + return repr(node[1]) if not float(node[1]).is_integer() else str(int(node[1])) + if kind == "sym": + return node[1] + if kind == "time": + return "time" + if kind == "call": + return f"{node[1]}({', '.join(infix(a) for a in node[2])})" + if kind == "piecewise": + body = ", ".join(f"{infix(v)} if {infix(c)}" for v, c in node[1]) + return f"piecewise({body}{', otherwise ' + infix(node[2]) if node[2] else ''})" + name, args = node[1], node[2] + if name in _INFIX and len(args) >= 2: + return "(" + f" {_INFIX[name]} ".join(infix(a) for a in args) + ")" + if name == "minus" and len(args) == 1: + return f"-{infix(args[0])}" + return f"{name}({', '.join(infix(a) for a in args)})" + + +def _substitute(node, mapping): + kind = node[0] + if kind == "sym": + return mapping.get(node[1], node) + if kind in ("op", "call"): + return (kind, node[1], tuple(_substitute(a, mapping) for a in node[2])) + if kind == "piecewise": + return ("piecewise", tuple((_substitute(v, mapping), _substitute(c, mapping)) for v, c in node[1]), + None if node[2] is None else _substitute(node[2], mapping)) + return node + + +def _inline(node, functions, depth=0): + """Replace calls to SBML function definitions by their bodies.""" + if depth > 64: + raise SBMLError("function definitions nest too deeply (recursive?)") + kind = node[0] + if kind == "call": + if node[1] not in functions: + raise SBMLError(f"call to undefined function {node[1]!r}") + params, body = functions[node[1]] + args = [_inline(a, functions, depth) for a in node[2]] + if len(args) != len(params): + raise SBMLError(f"function {node[1]!r} takes {len(params)} arguments, got {len(args)}") + return _inline(_substitute(body, dict(zip(params, args))), functions, depth + 1) + if kind == "op": + return ("op", node[1], tuple(_inline(a, functions, depth) for a in node[2])) + if kind == "piecewise": + return ("piecewise", tuple((_inline(v, functions, depth), _inline(c, functions, depth)) + for v, c in node[1]), + None if node[2] is None else _inline(node[2], functions, depth)) + return node + + +def compile_expr(node, index: dict[str, int], constants: dict[str, float] | None = None): + """Closure ``f(env, t) -> float``; ``index`` maps symbols to ``env`` slots and + ``constants`` binds local parameters (which shadow global symbols).""" + constants = constants or {} + kind = node[0] + if kind == "num": + v = float(node[1]) + return lambda env, t: v + if kind == "time": + return lambda env, t: t + if kind == "sym": + name = node[1] + if name in constants: + v = float(constants[name]) + return lambda env, t: v + if name not in index: + raise SBMLError(f"math references unknown or unsupported symbol {name!r}") + i = index[name] + return lambda env, t: env[i] + if kind == "piecewise": + pieces = [(compile_expr(v, index, constants), compile_expr(c, index, constants)) for v, c in node[1]] + other = compile_expr(node[2], index, constants) if node[2] is not None else (lambda env, t: math.nan) + + def piecewise(env, t): + for value, cond in pieces: + if cond(env, t): + return value(env, t) + return other(env, t) + return piecewise + if kind == "call": + raise SBMLError("function calls must be inlined before compilation") + fn = _OPERATORS[node[1]] + args = [compile_expr(a, index, constants) for a in node[2]] + if len(args) == 1: + a0 = args[0] + return lambda env, t: fn([a0(env, t)]) + if len(args) == 2: + a0, a1 = args + return lambda env, t: fn([a0(env, t), a1(env, t)]) + return lambda env, t: fn([a(env, t) for a in args]) + + +# --------------------------------------------------------------------------- +# model components +# --------------------------------------------------------------------------- +@dataclass(frozen=True) +class Compartment: + id: str + size: float | None + name: str = "" + + +@dataclass(frozen=True) +class Species: + id: str + compartment: str + initial: tuple[str, float] | None # ("amount" | "concentration", value) as written + has_only_substance_units: bool + boundary_condition: bool + constant: bool + name: str = "" + + +@dataclass(frozen=True) +class Reaction: + id: str + reactants: tuple[tuple[str, Fraction], ...] + products: tuple[tuple[str, Fraction], ...] + modifiers: tuple[str, ...] + rate: tuple | None # inlined expression tree, or None (no kinetic law) + local_parameters: dict = field(default_factory=dict) + reversible: bool = True + name: str = "" + + +def _bool(el, attr, default, level): + v = el.attrib.get(attr) + if v is None: + if level >= 3 and default is None: + raise SBMLError(f"<{_local(el.tag)} id={el.attrib.get('id')!r}> needs attribute {attr!r} in Level 3") + return bool(default) + if v not in ("true", "false", "1", "0"): + raise SBMLError(f"attribute {attr}={v!r} is not a boolean") + return v in ("true", "1") + + +def _float(el, attr): + v = el.attrib.get(attr) + if v is None: + return None + x = float(v) + if not math.isfinite(x): + raise SBMLError(f"attribute {attr}={v!r} is not finite") + return x + + +def _stoichiometry(ref, level) -> Fraction: + for c in ref: + if _local(c.tag) == "stoichiometryMath": + raise UnsupportedSBML("stoichiometryMath") + raw = ref.attrib.get("stoichiometry") + if raw is None: + if level >= 3: + raise UnsupportedSBML(f"Level 3 speciesReference to {ref.attrib.get('species')!r} without stoichiometry") + raw = "1" + try: + value = Fraction(raw.strip()) + except ValueError as exc: + raise SBMLError(f"stoichiometry {raw!r} is not a finite decimal") from exc + if value <= 0: + raise SBMLError(f"stoichiometry {raw!r} must be positive") + return value + + +_IGNORED = {"notes", "annotation", "listOfUnitDefinitions", "listOfCompartmentTypes", "listOfSpeciesTypes"} +_HANDLED = {"listOfFunctionDefinitions", "listOfCompartments", "listOfSpecies", "listOfParameters", + "listOfInitialAssignments", "listOfRules", "listOfReactions"} + + +class SBMLModel: + """An SBML reaction network: structure (exact stoichiometry), executable + kinetics, and conversion to a typed :class:`Graph`.""" + + def __init__(self, *, id, name, level, version, compartments, species, parameters, + functions, initial_assignments, rules, reactions, source=None): + self.id, self.name, self.level, self.version = id, name, level, version + self.compartments: dict[str, Compartment] = compartments + self.species: dict[str, Species] = species + self.parameters: dict[str, tuple[float | None, bool]] = parameters + self.functions = functions + self.initial_assignments: dict[str, tuple] = initial_assignments + self.rules: dict[str, tuple] = rules + self.reactions: dict[str, Reaction] = reactions + self.source = source or {} + self._check_semantics() + self._compile() + + # --- structure ----------------------------------------------------------- + def _check_semantics(self): + ids = list(self.compartments) + list(self.species) + list(self.parameters) + list(self.reactions) + dup = {i for i in ids if ids.count(i) > 1} + if dup: + raise SBMLError(f"duplicate ids {sorted(dup)}") + for s in self.species.values(): + if s.compartment not in self.compartments: + raise SBMLError(f"species {s.id!r} in unknown compartment {s.compartment!r}") + for target in self.rules: + if target in self.compartments: + raise UnsupportedSBML(f"assignment rule for compartment {target!r} (variable volume)") + if target not in self.species and target not in self.parameters: + raise SBMLError(f"assignment rule for unknown or reaction id {target!r}") + if target in self.initial_assignments: + raise SBMLError(f"{target!r} has both an assignment rule and an initial assignment") + if (target in self.parameters and self.parameters[target][1]) or ( + target in self.species and self.species[target].constant): + raise SBMLError(f"constant {target!r} cannot have an assignment rule") + for target in self.initial_assignments: + if target not in self.compartments and target not in self.species and target not in self.parameters: + raise UnsupportedSBML(f"initial assignment to {target!r} (not a compartment, species or parameter)") + for c in self.compartments.values(): + if c.size is None and c.id not in self.initial_assignments: + raise SBMLError(f"compartment {c.id!r} has no size") + for s in self.species.values(): + if s.initial is None and s.id not in self.initial_assignments and s.id not in self.rules: + raise SBMLError(f"species {s.id!r} has no initial value") + for p, (value, _) in self.parameters.items(): + if value is None and p not in self.initial_assignments and p not in self.rules: + raise SBMLError(f"parameter {p!r} has no value") + for r in self.reactions.values(): + for sid, _ in r.reactants + r.products: + if sid not in self.species: + raise SBMLError(f"reaction {r.id!r} references unknown species {sid!r}") + sp = self.species[sid] + if not sp.boundary_condition and (sp.constant or sid in self.rules): + raise SBMLError(f"species {sid!r} is changed by reaction {r.id!r} but is " + f"{'constant' if sp.constant else 'set by an assignment rule'} " + "and not a boundary condition") + for sid in r.modifiers: + if sid not in self.species: + raise SBMLError(f"reaction {r.id!r} modifier {sid!r} is not a species") + + @property + def dynamic_species(self) -> list[str]: + """Species whose amounts are state variables (changed only by reactions).""" + return [s.id for s in self.species.values() + if not s.boundary_condition and not s.constant and s.id not in self.rules] + + @property + def reaction_ids(self) -> list[str]: + return list(self.reactions) + + def stoichiometry_matrix(self, species_ids: list[str] | None = None) -> list[list[Fraction]]: + """Exact net stoichiometry N (rows: species, columns: reactions).""" + species_ids = self.dynamic_species if species_ids is None else species_ids + row = {s: i for i, s in enumerate(species_ids)} + n = [[Fraction(0)] * len(self.reactions) for _ in species_ids] + for j, r in enumerate(self.reactions.values()): + for sid, c in r.reactants: + if sid in row: + n[row[sid]][j] -= c + for sid, c in r.products: + if sid in row: + n[row[sid]][j] += c + return n + + def conservation_laws(self) -> list[dict[str, int]]: + """Exact conserved moieties over dynamic species amounts (see stoichiometry).""" + ids = self.dynamic_species + laws = S.conservation_laws(self.stoichiometry_matrix(ids)) + return [{ids[i]: c for i, c in enumerate(law) if c} for law in laws] + + def reacting_species(self) -> list[str]: + """Every species consumed or produced by some reaction, in document order.""" + used = {sid for r in self.reactions.values() for sid, _ in r.reactants + r.products} + return [s for s in self.species if s in used] + + def stoichiometric_consistency(self) -> dict: + """Positive-mass check over all reacting species. + + Boundary and constant species keep their rows (being held fixed does + not change a species' mass), but a reaction counts as internal only if + each side has at least one species that is not held fixed. Reactions + whose reactants or products are all fixed pools (for example + ``A -> EmptySet``, or synthesis from a clamped precursor) exchange + matter with the environment, like pure synthesis or degradation.""" + ids = self.reacting_species() + free = {s.id for s in self.species.values() if not s.boundary_condition and not s.constant} + internal = [any(sid in free for sid, _ in r.reactants) and any(sid in free for sid, _ in r.products) + for r in self.reactions.values()] + res = S.stoichiometric_consistency(self.stoichiometry_matrix(ids), internal) + return {**res, + "conserved": [ids[i] for i in res["conserved"]], + "unconserved": [ids[i] for i in res["unconserved"]], + "internal_reactions": [self.reaction_ids[j] for j in res["internal_reactions"]], + "mass": None if res["mass"] is None else {ids[i]: str(m) for i, m in enumerate(res["mass"])}, + "leak_combination": None if res["leak_combination"] is None else { + self.reaction_ids[j]: str(v) for j, v in zip(res["internal_reactions"], res["leak_combination"]) + if v}} + + # --- kinetics ------------------------------------------------------------ + def _compile(self): + self._symbols = list(self.compartments) + list(self.species) + list(self.parameters) + self._index = {s: i for i, s in enumerate(self._symbols)} + allowed = set(self._index) + + def check_refs(expr, where, local=()): + unknown = symbols(expr) - allowed - set(local) + if unknown & set(self.reactions): + raise UnsupportedSBML(f"{where} uses a reaction id as a rate symbol") + if unknown: + raise SBMLError(f"{where} references unknown symbols {sorted(unknown)}") + + for target, expr in {**self.rules, **self.initial_assignments}.items(): + check_refs(expr, f"rule/assignment for {target!r}") + self._rule_order = _topological(self.rules, set(self.rules)) + self._rule_fns = [(self._index[t], compile_expr(self.rules[t], self._index)) for t in self._rule_order] + self._rate_fns = [] + for r in self.reactions.values(): + if r.rate is None: + self._rate_fns.append(None) + continue + check_refs(r.rate, f"kinetic law of {r.id!r}", r.local_parameters) + self._rate_fns.append(compile_expr(r.rate, self._index, r.local_parameters)) + self._missing_laws = [r for r, fn in zip(self.reactions, self._rate_fns) if fn is None] + dyn = self.dynamic_species + self._dyn_env = np.array([self._index[s] for s in dyn], dtype=int) + self._dyn_size_ix = np.array([self._index[self.species[s].compartment] for s in dyn], dtype=int) + self._dyn_conc = np.array([not self.species[s].has_only_substance_units for s in dyn], dtype=bool) + n = self.stoichiometry_matrix(dyn) + self._n = np.array([[float(x) for x in row] for row in n]).reshape(len(dyn), len(self.reactions)) + self._env0 = self._initial_env() + + def _initial_env(self) -> np.ndarray: + """Symbol values at t = 0: attributes, then initial assignments and rules in + dependency order. Species values are concentrations unless + hasOnlySubstanceUnits.""" + env = np.full(len(self._symbols), np.nan) + defs: dict[str, tuple] = {} + for c in self.compartments.values(): + if c.size is not None: + env[self._index[c.id]] = c.size + for p, (value, _) in self.parameters.items(): + if value is not None: + env[self._index[p]] = value + raw_species = {} + for s in self.species.values(): + if s.id in self.initial_assignments or s.id in self.rules: + continue + raw_species[s.id] = s + defs[s.id] = ("sym", s.compartment) # dependency marker only + defs.update(self.initial_assignments) + defs.update(self.rules) + order = _topological(defs, set(defs)) + for target in order: + i = self._index[target] + if target in raw_species: + s = raw_species[target] + kind, value = s.initial + size = env[self._index[s.compartment]] + as_amount = value if kind == "amount" else value * size + env[i] = as_amount if s.has_only_substance_units else _ieee_div(as_amount, size) + else: + env[i] = compile_expr(_inline(defs[target], self.functions), self._index)(env, 0.0) + if np.isnan(env).any(): + missing = [s for s, v in zip(self._symbols, env) if np.isnan(v)] + raise SBMLError(f"no initial value for {missing}") + return env + + def sizes(self) -> np.ndarray: + """Compartment size of each dynamic species (state scaling).""" + return self._env0[self._dyn_size_ix].copy() + + def initial_amounts(self) -> np.ndarray: + """Initial amount of each dynamic species, in model substance units.""" + v = self._env0[self._dyn_env] + return np.where(self._dyn_conc, v * self.sizes(), v) + + def _env(self, amounts, t) -> np.ndarray: + env = self._env0.copy() + sizes = env[self._dyn_size_ix] + env[self._dyn_env] = np.where(self._dyn_conc, amounts / sizes, amounts) + for i, fn in self._rule_fns: + env[i] = fn(env, t) + return env + + def rates(self, amounts, t: float = 0.0) -> np.ndarray: + """Reaction rates (substance / time) at the given dynamic species amounts.""" + if self._missing_laws: + raise SBMLError(f"reactions without kinetic laws cannot be simulated: {self._missing_laws}") + env = self._env(np.asarray(amounts, dtype=float), t) + return np.array([fn(env, t) for fn in self._rate_fns], dtype=float) + + def amount_derivative(self, amounts, t: float = 0.0) -> np.ndarray: + return self._n @ self.rates(amounts, t) if len(self.reactions) else np.zeros(len(amounts)) + + def simulate(self, times, rtol=1e-10, atol=1e-12, method="LSODA", amounts0=None) -> dict: + """Integrate from t = times[0]. The integrated state is amount / size + (concentration-scaled) so tolerances are comparable across compartments.""" + from scipy.integrate import solve_ivp + + times = np.asarray(times, dtype=float) + sizes = self.sizes() + y0 = (self.initial_amounts() if amounts0 is None else np.asarray(amounts0, float)) / sizes + sol = solve_ivp(lambda t, y: self.amount_derivative(y * sizes, t) / sizes, + (times[0], times[-1]), y0, t_eval=times, method=method, rtol=rtol, atol=atol) + if not sol.success: + raise RuntimeError(f"integration failed: {sol.message}") + amounts = (sol.y.T * sizes) + conc = np.where(self._dyn_conc, amounts / sizes, amounts) + return {"time": sol.t, "species": self.dynamic_species, "amount": amounts, + "value": conc, "nfev": int(sol.nfev)} + + def _reduction(self): + """Exact conservation laws in RREF on amounts: x_P = T - M x_Q.""" + ids = self.dynamic_species + laws = S.conservation_laws(self.stoichiometry_matrix(ids)) + if not laws: + return [], [], list(range(len(ids))), np.zeros((0, len(ids))) + r, pivots = S.rref(laws) + free = [c for c in range(len(ids)) if c not in pivots] + return r, pivots, free, np.array([[float(x) for x in row] for row in r]) + + def jacobian(self, amounts, t: float = 0.0, rel_step=1e-7) -> np.ndarray: + """d(amount derivative)/d(amount) by central differences.""" + x = np.asarray(amounts, dtype=float) + scale = np.maximum(np.abs(x), self.sizes() * 1e-9) + j = np.empty((len(x), len(x))) + for k in range(len(x)): + h = rel_step * scale[k] + xp, xm = x.copy(), x.copy() + xp[k] += h + xm[k] -= h + j[:, k] = (self.amount_derivative(xp, t) - self.amount_derivative(xm, t)) / (2 * h) + return j + + def steady_state(self, amounts_guess=None, tol=1e-13) -> dict: + """Newton/hybrid solve of N v(x) = 0 on the conservation-reduced system. + + The conserved totals are fixed at those of ``amounts_guess`` (default: + the initial amounts), so the answer is the steady state of that + stoichiometric class, nearest the guess when several exist. + """ + from scipy.optimize import root + + x0 = self.initial_amounts() if amounts_guess is None else np.asarray(amounts_guess, float) + sizes = self.sizes() + _, pivots, free, lr = self._reduction() + totals = lr @ x0 if len(pivots) else np.zeros(0) + + def residual(y): + x = y * sizes + f = self.amount_derivative(x) / sizes + return np.concatenate([f[free], (lr @ x - totals) / np.maximum(np.abs(totals), 1e-300)]) + + sol = root(residual, x0 / sizes, method="hybr", tol=tol) + x = sol.x * sizes + f = self.amount_derivative(x) + return {"amount": x, "species": self.dynamic_species, "success": bool(sol.success), + "max_abs_derivative": float(np.max(np.abs(f / sizes))) if len(f) else 0.0, + "conservation_residual": float(np.max(np.abs(lr @ x - totals))) if len(pivots) else 0.0} + + def reduced_eigenvalues(self, amounts) -> np.ndarray: + """Eigenvalues of the Jacobian restricted to the stoichiometric class + (the conservation laws' zero modes removed via the link matrix).""" + j = self.jacobian(amounts) + _, pivots, free, lr = self._reduction() + link = np.zeros((len(amounts), len(free))) + for col, q in enumerate(free): + link[q, col] = 1.0 + for row, p in enumerate(pivots): + link[p, :] = -lr[row, free] + return np.linalg.eigvals(j[np.ix_(free, range(len(amounts)))] @ link) + + # --- graph --------------------------------------------------------------- + def to_graph(self) -> tuple[Graph, Schema]: + """Typed graph: compartments, species, reactions and exact conserved + moieties, with the schema and invariants in :func:`sbml_schema`.""" + g = Graph(self.id or "sbml") + env0 = self._env0 + for c in self.compartments.values(): + g.add_node(c.id, "compartment", "cell", c.name or c.id, size=float(env0[self._index[c.id]])) + x0 = dict(zip(self.dynamic_species, self.initial_amounts())) + for s in self.species.values(): + v = env0[self._index[s.id]] + amount = v if s.has_only_substance_units else v * env0[self._index[s.compartment]] + g.add_node(s.id, "species", "molecule", s.name or s.id, + compartment=s.compartment, initial_amount=float(x0.get(s.id, amount)), + boundary_condition=s.boundary_condition, constant=s.constant, + has_only_substance_units=s.has_only_substance_units, + assignment_rule=infix(self.rules[s.id]) if s.id in self.rules else None) + g.add_edge(s.id, s.compartment, "located_in") + for r in self.reactions.values(): + g.add_node(r.id, "reaction", "cell", r.name or r.id, reversible=r.reversible, + kinetic_law=None if r.rate is None else infix(r.rate), + local_parameters={k: float(v) for k, v in r.local_parameters.items()}) + key = f"sbml:{self.id}/{r.id}" + for sid, c in r.reactants: + g.add_edge(sid, r.id, "reactant", model=key, stoichiometry=str(c)) + for sid, c in r.products: + g.add_edge(r.id, sid, "product", model=key, stoichiometry=str(c)) + for sid in r.modifiers: + g.add_edge(sid, r.id, "modifier", model=key) + for k, law in enumerate(self.conservation_laws()): + total = sum(c * x0[s] for s, c in law.items()) + mid = f"moiety_{k + 1}" + g.add_node(mid, "moiety", "molecule", " + ".join( + (f"{c}*{s}" if c != 1 else s) for s, c in law.items()), total=float(total)) + for s, c in law.items(): + g.add_edge(s, mid, "member", coefficient=str(c)) + return g, sbml_schema() + + +def _topological(defs: dict[str, tuple], targets: set[str]) -> list[str]: + """Order definitions so each follows the targets it references.""" + deps = {t: symbols(e) & targets - {t} for t, e in defs.items()} + for t, e in defs.items(): + if t in symbols(e) and e[0] != "sym": + raise SBMLError(f"{t!r} is defined in terms of itself") + order, done, active = [], set(), set() + + def visit(t): + if t in done: + return + if t in active: + raise SBMLError(f"algebraic loop through {t!r}") + active.add(t) + for d in sorted(deps[t]): + visit(d) + active.discard(t) + done.add(t) + order.append(t) + + for t in defs: + visit(t) + return order + + +# --------------------------------------------------------------------------- +# reading +# --------------------------------------------------------------------------- +def read_sbml(source) -> SBMLModel: + """Read an SBML document from a path or bytes. A path's SHA-256 is recorded.""" + if isinstance(source, (bytes, bytearray)): + data, origin = bytes(source), None + else: + origin = Path(source) + data = origin.read_bytes() + head = data[:4096].lower() + if b"") + ns = root.tag[1:].split("}")[0] if root.tag.startswith("{") else "" + level, version = int(root.attrib.get("level", 0)), int(root.attrib.get("version", 0)) + if (level, version) not in {(2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (3, 1), (3, 2)}: + raise UnsupportedSBML(f"SBML Level {level} Version {version}") + for attr, value in root.attrib.items(): + if _local(attr) == "required" and value == "true": + raise UnsupportedSBML(f"required SBML package ({attr})") + + def q(tag): + return f"{{{ns}}}{tag}" if ns else tag + + model = root.find(q("model")) + if model is None: + raise SBMLError("no ") + if level >= 3 and "conversionFactor" in model.attrib: + raise UnsupportedSBML("model conversionFactor") + for child in model: + name = _local(child.tag) + if name in _IGNORED or (name in _HANDLED and child.tag == q(name)): + continue + if len(child) == 0 and name.startswith("listOf"): + continue + raise UnsupportedSBML(f"model element <{name}>") + + def items(list_tag, item_tag): + lst = model.find(q(list_tag)) + if lst is None: + return [] + out = [] + for el in lst: + if _local(el.tag) in ("notes", "annotation"): + continue + if el.tag != q(item_tag): + raise UnsupportedSBML(f"<{_local(el.tag)}> in <{list_tag}>") + out.append(el) + return out + + def math_of(el, where): + m = el.find(f"{{{MATHML}}}math") + if m is None: + raise SBMLError(f"{where} has no MathML") + return parse_mathml(m) + + functions = {} + for fd in items("listOfFunctionDefinitions", "functionDefinition"): + m = fd.find(f"{{{MATHML}}}math") + if m is None: + raise SBMLError(f"function {fd.attrib.get('id')!r} has no MathML") + lam_el = [c for c in m if _local(c.tag) not in ("annotation", "annotation-xml")] + if len(lam_el) != 1 or _local(lam_el[0].tag) != "lambda": + raise SBMLError(f"function {fd.attrib.get('id')!r} is not a lambda") + kids = list(lam_el[0]) + params = [] + for k in kids[:-1]: + if _local(k.tag) != "bvar" or len(k) != 1: + raise SBMLError(f"function {fd.attrib.get('id')!r}: malformed ") + params.append("".join(k[0].itertext()).strip()) + functions[fd.attrib["id"]] = (tuple(params), parse_mathml(kids[-1])) + + compartments = {} + for c in items("listOfCompartments", "compartment"): + if not _bool(c, "constant", True if level < 3 else None, level): + raise UnsupportedSBML(f"non-constant compartment {c.attrib['id']!r}") + size = _float(c, "size") + compartments[c.attrib["id"]] = Compartment(c.attrib["id"], size, c.attrib.get("name", "")) + + species = {} + for s in items("listOfSpecies", "species"): + if "conversionFactor" in s.attrib: + raise UnsupportedSBML(f"species conversionFactor on {s.attrib['id']!r}") + amount, conc = _float(s, "initialAmount"), _float(s, "initialConcentration") + if amount is not None and conc is not None: + raise SBMLError(f"species {s.attrib['id']!r} has both initialAmount and initialConcentration") + initial = ("amount", amount) if amount is not None else ("concentration", conc) if conc is not None else None + species[s.attrib["id"]] = Species( + s.attrib["id"], s.attrib.get("compartment", ""), initial, + _bool(s, "hasOnlySubstanceUnits", False if level < 3 else None, level), + _bool(s, "boundaryCondition", False if level < 3 else None, level), + _bool(s, "constant", False if level < 3 else None, level), + s.attrib.get("name", "")) + + parameters = {} + for p in items("listOfParameters", "parameter"): + parameters[p.attrib["id"]] = (_float(p, "value"), _bool(p, "constant", True if level < 3 else None, level)) + + initial_assignments = {} + for ia in items("listOfInitialAssignments", "initialAssignment"): + initial_assignments[ia.attrib["symbol"]] = _inline(math_of(ia, "initialAssignment"), functions) + + rules = {} + lst = model.find(q("listOfRules")) + for rule in (lst if lst is not None else []): + kind = _local(rule.tag) + if kind in ("notes", "annotation"): + continue + if kind != "assignmentRule": + raise UnsupportedSBML(f"<{kind}>") + var = rule.attrib["variable"] + if var in rules: + raise SBMLError(f"two assignment rules for {var!r}") + rules[var] = _inline(math_of(rule, f"rule for {var!r}"), functions) + + reactions = {} + for r in items("listOfReactions", "reaction"): + rid = r.attrib["id"] + if _bool(r, "fast", False if (level < 3 or version >= 2) else None, level): + raise UnsupportedSBML(f"fast reaction {rid!r}") + def refs(tag, item): + lst = r.find(q(tag)) + out = [] + for el in (lst if lst is not None else []): + if _local(el.tag) in ("notes", "annotation"): + continue + if el.tag != q(item): + raise UnsupportedSBML(f"<{_local(el.tag)}> in <{tag}>") + out.append(el) + return out + + reactants = tuple((x.attrib["species"], _stoichiometry(x, level)) + for x in refs("listOfReactants", "speciesReference")) + products = tuple((x.attrib["species"], _stoichiometry(x, level)) + for x in refs("listOfProducts", "speciesReference")) + modifiers = tuple(x.attrib["species"] for x in refs("listOfModifiers", "modifierSpeciesReference")) + kl = r.find(q("kineticLaw")) + rate, local = None, {} + if kl is not None: + for tag, item in (("listOfParameters", "parameter"), ("listOfLocalParameters", "localParameter")): + lst = kl.find(q(tag)) + for p in (lst if lst is not None else []): + if _local(p.tag) in ("notes", "annotation"): + continue + if p.tag != q(item): + raise UnsupportedSBML(f"<{_local(p.tag)}> in kinetic law of {rid!r}") + v = _float(p, "value") + if v is None: + raise SBMLError(f"local parameter {p.attrib['id']!r} of {rid!r} has no value") + local[p.attrib["id"]] = v + rate = _inline(math_of(kl, f"kinetic law of {rid!r}"), functions) + reactions[rid] = Reaction(rid, reactants, products, modifiers, rate, local, + _bool(r, "reversible", True if level < 3 else None, level), + r.attrib.get("name", "")) + + src = {"sha256": hashlib.sha256(data).hexdigest()} + if origin is not None: + src["path"] = str(origin) + return SBMLModel(id=model.attrib.get("id", ""), name=model.attrib.get("name", ""), + level=level, version=version, compartments=compartments, species=species, + parameters=parameters, functions=functions, + initial_assignments=initial_assignments, rules=rules, + reactions=reactions, source=src) + + +# --------------------------------------------------------------------------- +# schema and invariants for an SBML-derived graph +# --------------------------------------------------------------------------- +def _graph_stoichiometry(g: Graph, dynamic_only: bool = True): + """N recomputed from graph edges alone (independent of SBMLModel): rows are + the dynamic species, or with ``dynamic_only=False`` every reacting species.""" + if dynamic_only: + dyn = [n.id for n in g.nodes_of_type("species") + if not n.attrs.get("boundary_condition") and not n.attrs.get("constant") + and n.attrs.get("assignment_rule") is None] + else: + used = {e.src for e in g.edges_of_type("reactant")} | {e.dst for e in g.edges_of_type("product")} + dyn = [n.id for n in g.nodes_of_type("species") if n.id in used] + rx = [n.id for n in g.nodes_of_type("reaction")] + row, col = {s: i for i, s in enumerate(dyn)}, {r: j for j, r in enumerate(rx)} + n = [[Fraction(0)] * len(rx) for _ in dyn] + for e in g.edges_of_type("reactant"): + if e.src in row: + n[row[e.src]][col[e.dst]] -= Fraction(e.attrs["stoichiometry"]) + for e in g.edges_of_type("product"): + if e.dst in row: + n[row[e.dst]][col[e.src]] += Fraction(e.attrs["stoichiometry"]) + return dyn, rx, n + + +def sbml_schema() -> Schema: + schema = Schema( + node_types={"compartment": {"size"}, "species": {"compartment", "initial_amount", "boundary_condition"}, + "reaction": {"reversible"}, "moiety": {"total"}}, + edge_rules=[ + EdgeRule("reactant", {"species"}, {"reaction"}, {"stoichiometry"}), + EdgeRule("product", {"reaction"}, {"species"}, {"stoichiometry"}), + EdgeRule("modifier", {"species"}, {"reaction"}), + EdgeRule("located_in", {"species"}, {"compartment"}), + EdgeRule("member", {"species"}, {"moiety"}, {"coefficient"}), + ]) + + @schema.invariant("stoichiometry") + def _stoich(g: Graph) -> list[str]: + out = [] + for e in g.edges_of_type("reactant", "product"): + try: + c = Fraction(str(e.attrs.get("stoichiometry"))) + ok = c > 0 + except (TypeError, ValueError, ZeroDivisionError): + ok = False + if not ok: + out.append(f"edge {e.key}: stoichiometry must be a positive exact rational") + return out + + @schema.invariant("located") + def _located(g: Graph) -> list[str]: + return [f"species {n.id}: needs exactly one compartment" + for n in g.nodes_of_type("species") if len(g.out_edges(n.id, ("located_in",))) != 1] + + @schema.invariant("conservation") + def _conservation(g: Graph) -> list[str]: + """Each moiety node must be an exact left-null vector of N, the moieties + must be independent and span the whole left nullspace, and each total + must equal the initial amounts it weights.""" + try: + dyn, _, n = _graph_stoichiometry(g) + except (KeyError, ValueError, ZeroDivisionError) as exc: + return [f"cannot rebuild stoichiometry: {exc}"] + out, laws = [], [] + idx = {s: i for i, s in enumerate(dyn)} + for m in g.nodes_of_type("moiety"): + law = [Fraction(0)] * len(dyn) + for e in g.in_edges(m.id, ("member",)): + if e.src not in idx: + out.append(f"{m.id}: member {e.src} is not a dynamic species") + continue + law[idx[e.src]] = Fraction(e.attrs["coefficient"]) + if S.conservation_residual([law], n) != 0: + out.append(f"{m.id}: not conserved by every reaction (l.N != 0)") + total = sum(float(c) * g.nodes[s].attrs["initial_amount"] for s, c in zip(dyn, law) if c) + if not math.isclose(total, m.attrs["total"], rel_tol=1e-12, abs_tol=1e-300): + out.append(f"{m.id}: total {m.attrs['total']} != weighted initial amounts {total}") + laws.append(law) + expected = len(dyn) - S.rank(n) if n and n[0] else len(dyn) + if laws and S.rank(laws) != len(laws): + out.append("moieties are linearly dependent") + if len(laws) != expected: + out.append(f"{len(laws)} moieties recorded, but the stoichiometry has {expected} " + "independent conservation laws") + return out + + @schema.invariant("stoichiometric-consistency") + def _consistency(g: Graph) -> list[str]: + """Every dynamic species must admit a positive mass conserved by all + internal reactions (Gevorgyan et al. 2008).""" + dyn, rx, n = _graph_stoichiometry(g, dynamic_only=False) + fixed = {v.id for v in g.nodes_of_type("species") + if v.attrs.get("boundary_condition") or v.attrs.get("constant")} + internal = [any(e.src not in fixed for e in g.in_edges(r, ("reactant",))) + and any(e.dst not in fixed for e in g.out_edges(r, ("product",))) for r in rx] + res = S.stoichiometric_consistency(n, internal) + out = [] + if res["unconserved"]: + names = [dyn[i] for i in res["unconserved"]] + leak = res["leak_combination"] + how = ("exact certificate: net production by " + ", ".join( + f"{v}*{rx[j]}" for j, v in zip(res["internal_reactions"], leak) if v) + if leak else "LP only; no exact certificate") + out.append(f"species {names} cannot carry a positive conserved mass ({how})") + elif not res["certified"]: + out.append("consistency found by LP but no exact mass vector could be certified") + return out + + return schema diff --git a/vh/graph/stoichiometry.py b/vh/graph/stoichiometry.py new file mode 100644 index 0000000..1847175 --- /dev/null +++ b/vh/graph/stoichiometry.py @@ -0,0 +1,259 @@ +"""Exact stoichiometric analysis: conservation laws and stoichiometric consistency. + +Stoichiometric coefficients are exact rationals (``fractions.Fraction``), so the +linear algebra here has no rounding: a conservation law ``l`` is accepted only if +``l . N`` is exactly zero, and the number of laws is exactly ``n - rank(N)``. + +Stoichiometric consistency (Gevorgyan, Poolman & Fell, Bioinformatics 24:2245, +2008) asks whether every species can be given a strictly positive mass such that +every internal reaction conserves mass. It is a linear program; the LP is solved +in floating point (SciPy/HiGHS), then each answer is turned into an exact +rational certificate and checked exactly: + + - a conserved species gets a nonnegative mass vector ``m`` with + ``m^T N_internal = 0`` exactly and ``m_i > 0``; + - an unconserved species gets a reaction combination ``v`` with + ``N_internal v >= 0`` exactly and ``(N_internal v)_i > 0``. For any + nonnegative ``m`` with ``m^T N_internal = 0``, ``0 = m^T N v >= m_i (N v)_i``, + so ``m_i = 0``: no positive mass exists (Farkas' lemma). + +If rounding prevents an exact certificate, the result says so instead of +reporting the floating-point answer as proven. +""" + +from __future__ import annotations + +from fractions import Fraction +from math import gcd, lcm +from numbers import Rational + +Matrix = list[list[Fraction]] + + +def as_fraction(value) -> Fraction: + """Exact rational from an int, Fraction or decimal string; floats are rejected + because their binary value is rarely the coefficient the author meant.""" + if isinstance(value, bool): + raise TypeError("a boolean is not a stoichiometric coefficient") + if isinstance(value, (Rational, str)): + return Fraction(value) + raise TypeError(f"exact rational required, got {type(value).__name__}") + + +def _matrix(rows) -> Matrix: + out = [[as_fraction(x) for x in row] for row in rows] + if out and any(len(r) != len(out[0]) for r in out): + raise ValueError("matrix rows must have equal length") + return out + + +def transpose(rows: Matrix, ncols: int | None = None) -> Matrix: + if not rows: + return [[] for _ in range(ncols or 0)] + return [list(col) for col in zip(*rows)] + + +def rref(rows) -> tuple[Matrix, list[int]]: + """Reduced row echelon form over the rationals, and its pivot columns.""" + a = [list(r) for r in _matrix(rows)] + if not a: + return [], [] + nrows, ncols = len(a), len(a[0]) + pivots: list[int] = [] + r = 0 + for c in range(ncols): + p = next((i for i in range(r, nrows) if a[i][c] != 0), None) + if p is None: + continue + a[r], a[p] = a[p], a[r] + inv = 1 / a[r][c] + a[r] = [x * inv for x in a[r]] + for i in range(nrows): + if i != r and a[i][c] != 0: + f = a[i][c] + a[i] = [x - f * y for x, y in zip(a[i], a[r])] + pivots.append(c) + r += 1 + if r == nrows: + break + return a[:r], pivots + + +def rank(rows) -> int: + return len(rref(rows)[1]) + + +def nullspace(rows, ncols: int | None = None) -> Matrix: + """Basis of {v : A v = 0} (one vector per free column of the RREF).""" + a = _matrix(rows) + n = len(a[0]) if a else (ncols or 0) + r, pivots = rref(a) + free = [c for c in range(n) if c not in pivots] + basis = [] + for f in free: + v = [Fraction(0)] * n + v[f] = Fraction(1) + for row, p in zip(r, pivots): + v[p] = -row[f] + basis.append(v) + return basis + + +def integer_row(row) -> list[int]: + """Scale a rational row to coprime integers with a positive first nonzero.""" + row = [as_fraction(x) for x in row] + den = lcm(*(x.denominator for x in row)) if row else 1 + ints = [int(x * den) for x in row] + g = gcd(*ints) if any(ints) else 1 + ints = [x // g for x in ints] + first = next((x for x in ints if x), 0) + return [-x for x in ints] if first < 0 else ints + + +def conservation_laws(stoichiometry) -> list[list[int]]: + """Canonical integer basis of the left nullspace of N (species x reactions). + + Each row ``l`` satisfies ``l . N = 0`` exactly, so ``sum_i l_i x_i`` is + constant along every trajectory of ``dx/dt = N v(x)``. The basis is the RREF + of the left nullspace (unique for the space), scaled to coprime integers, so + two implementations can be compared row for row. + """ + n = _matrix(stoichiometry) + nspecies = len(n) + nreactions = len(n[0]) if n else 0 + if nreactions == 0: + return [integer_row([int(i == j) for j in range(nspecies)]) for i in range(nspecies)] + left = nullspace(transpose(n), nspecies) + if not left: + return [] + canonical, _ = rref(left) + return [integer_row(r) for r in canonical] + + +def conservation_residual(laws, stoichiometry) -> Fraction: + """max |l . N_j| over laws and reactions, computed exactly (0 means proven).""" + n = _matrix(stoichiometry) + worst = Fraction(0) + for law in laws: + law = [as_fraction(x) for x in law] + for j in range(len(n[0]) if n else 0): + worst = max(worst, abs(sum(law[i] * n[i][j] for i in range(len(n))))) + return worst + + +def _matmul(a: Matrix, b: Matrix) -> Matrix: + return [[sum(x * y for x, y in zip(row, col)) for col in zip(*b)] for row in a] if b and b[0] else [[] for _ in a] + + +def _certify(b: Matrix, positive: list[int], max_rounds: int = 64): + """Exact y with w = B y >= 0 and w_i > 0 for every i in ``positive``. + + Solves an LP in floating point (maximise the number of positive entries + requested), rationalises the answer and checks it exactly. Entries that + the LP leaves numerically at zero are then constrained to be exactly zero + by restricting y to the exact nullspace of those rows, and the LP is + solved again; each round fixes at least one more row, so the loop ends. + Returns (y, w) as Fractions, or None when no exact certificate is found. + """ + from scipy.optimize import linprog + + nrows = len(b) + ncols = len(b[0]) if b else 0 + zero: set[int] = set() + for _ in range(max_rounds): + if zero: + k = nullspace([b[i] for i in sorted(zero)], ncols) + if not k: + return None + basis = transpose(k) # ncols x dim + else: + basis = [[Fraction(int(i == j)) for j in range(ncols)] for i in range(ncols)] + m = _matmul(b, basis) # nrows x dim + dim = len(basis[0]) if basis else 0 + free_rows = [i for i in range(nrows) if i not in zero] + want = [i for i in positive if i not in zero] + if len(want) < len(positive): + return None + mf = [[float(x) for x in m[i]] for i in range(nrows)] + # variables [y (dim, free), z (len(want), in [0,1])]: M y >= 0 on free rows; z_k <= (M y)_want_k + a_ub = [[-x for x in mf[i]] + [0.0] * len(want) for i in free_rows] + a_ub += [[-x for x in mf[i]] + [1.0 if c == k else 0.0 for c in range(len(want))] + for k, i in enumerate(want)] + lp = linprog(c=[0.0] * dim + [-1.0] * len(want), A_ub=a_ub or None, + b_ub=[0.0] * len(a_ub) or None, + bounds=[(None, None)] * dim + [(0, 1)] * len(want), method="highs") + if lp.status != 0: + return None + y = lp.x[:dim] + wf = [sum(mf[i][c] * y[c] for c in range(dim)) for i in range(nrows)] + scale = max((abs(x) for x in wf), default=0.0) or 1.0 + yq = [Fraction(float(v) / scale).limit_denominator(10**9) for v in y] + w = [sum(m[i][c] * yq[c] for c in range(dim)) for i in range(nrows)] + if all(x >= 0 for x in w) and all(w[i] > 0 for i in positive): + return [sum(basis[j][c] * yq[c] for c in range(dim)) for j in range(ncols)], w + tiny = {i for i in free_rows if i not in positive and abs(wf[i]) / scale < 1e-9} + if not tiny - zero: + return None + zero |= tiny + return None + + +def stoichiometric_consistency(stoichiometry, internal: list[bool] | None = None) -> dict: + """Which species can carry a strictly positive conserved mass? + + Rows of ``stoichiometry`` should be every species that takes part in a + reaction, including boundary or constant species: holding a species fixed + changes its dynamics, not its mass. ``internal[j]`` marks reactions that + should conserve mass; by default a reaction is internal when it has both a + consumed and a produced species (pure synthesis/degradation/exchange + reactions are open by construction). + + Returns ``consistent`` (every species conservable), the conserved and + unconserved species indices, the exact certificates described in the + module docstring, and ``certified`` (True only if every classification has + an exactly verified certificate). + """ + from scipy.optimize import linprog + + n = _matrix(stoichiometry) + ns = len(n) + nr = len(n[0]) if n else 0 + if internal is None: + internal = [any(n[i][j] < 0 for i in range(ns)) and any(n[i][j] > 0 for i in range(ns)) + for j in range(nr)] + cols = [j for j in range(nr) if internal[j]] + sub = [[n[i][j] for j in cols] for i in range(ns)] + result = {"internal_reactions": cols, "n_species": ns} + if ns == 0 or not cols: + return {**result, "consistent": True, "conserved": list(range(ns)), "unconserved": [], + "mass": [Fraction(1)] * ns, "leak_combination": [], "certified": True} + + # Classification LP: maximise the number of species with positive mass. + # variables [m (ns, >= 0), z (ns, in [0, 1])]; m^T N_int = 0; z <= m. + a_eq = [[float(sub[i][k]) for i in range(ns)] + [0.0] * ns for k in range(len(cols))] + a_ub = [[-1.0 if c == i else 0.0 for c in range(ns)] + [1.0 if c == i else 0.0 for c in range(ns)] + for i in range(ns)] + lp = linprog(c=[0.0] * ns + [-1.0] * ns, A_ub=a_ub, b_ub=[0.0] * ns, A_eq=a_eq, + b_eq=[0.0] * len(cols), bounds=[(0, None)] * ns + [(0, 1)] * ns, method="highs") + if lp.status != 0: + raise RuntimeError(f"consistency LP failed: {lp.message}") + conserved = [i for i in range(ns) if lp.x[ns + i] > 0.5] + unconserved = [i for i in range(ns) if i not in conserved] + + # Mass certificate: m = B c in the exact left nullspace of N_int, m >= 0, m_conserved > 0. + mass = None + left = nullspace(transpose(sub), ns) + if left and conserved: + cert = _certify(transpose(left), conserved) + if cert is not None: + mass = cert[1] + elif not conserved: + mass = [Fraction(0)] * ns + # Leak certificate: v with N_int v >= 0 and (N_int v)_unconserved > 0. + leak: list[Fraction] | None = [] + if unconserved: + cert = _certify(sub, unconserved) + leak = None if cert is None else cert[0] + return {**result, "consistent": not unconserved, "conserved": conserved, + "unconserved": unconserved, "mass": mass, "leak_combination": leak, + "certified": mass is not None and leak is not None} From 05dd3177b8dda43f71dbd0d0f9f209efd3c176a7 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 5 Nov 2026 16:20:31 +0000 Subject: [PATCH 2/2] Refuse DTDs in any encoding, not only ASCII-compatible bytes The DOCTYPE/entity refusal scanned raw bytes for " Claude-Session: https://claude.ai/code/session_01CoFhtFtE7wJgvHKAzNB3xD --- tests/test_sbml.py | 24 ++++++++++++++++++++++++ vh/graph/sbml.py | 27 ++++++++++++++++++++++++--- 2 files changed, 48 insertions(+), 3 deletions(-) diff --git a/tests/test_sbml.py b/tests/test_sbml.py index 069e173..3546379 100644 --- a/tests/test_sbml.py +++ b/tests/test_sbml.py @@ -187,6 +187,30 @@ def test_doctype_and_entity_declarations_are_refused(prefix): read_sbml(prefix + l3(BASE).split(b"\n", 1)[1]) +UTF16_ENTITY = ('\n]>\n' + '' + '') + + +@pytest.mark.parametrize("encoding", ["utf-16", "utf-16-le", "utf-16-be"]) +def test_entity_refusal_does_not_depend_on_the_byte_encoding(encoding): + # Regression: an ASCII byte scan missed these, and the entity was expanded. + with pytest.raises(SBMLError, match="DOCTYPE"): + read_sbml(UTF16_ENTITY.encode(encoding)) + + +def test_utf16_documents_without_a_doctype_still_load(): + doc = l3(BASE + "" + DECAY + "").decode() + m = read_sbml(doc.replace('encoding="UTF-8"', 'encoding="UTF-16"').encode("utf-16")) + assert m.dynamic_species == ["A", "B"] + + +@pytest.mark.parametrize("data", [b"", b"".encode("utf-32")]) +def test_unreadable_xml_is_an_sbml_error(data): + with pytest.raises(SBMLError): + read_sbml(data) + + @pytest.mark.parametrize("doc,message", [ (l3(COMP.format(size=1) + "" + species("A", conc=1, constant="true") + species("B", conc=0) + "" + DECAY + ""), "constant"), diff --git a/vh/graph/sbml.py b/vh/graph/sbml.py index 6968de4..1807b06 100644 --- a/vh/graph/sbml.py +++ b/vh/graph/sbml.py @@ -768,6 +768,29 @@ def visit(t): # --------------------------------------------------------------------------- # reading # --------------------------------------------------------------------------- +def _refuse_dtd(data: bytes) -> None: + """Refuse any DOCTYPE or entity declaration before the real parse. + + Expat decodes the document in whatever encoding it declares (UTF-8, + UTF-16 with or without a byte-order mark, Latin-1, ...), so the check + sees the declarations themselves. A scan of the raw bytes would not: in + UTF-16, `` SBMLModel: """Read an SBML document from a path or bytes. A path's SHA-256 is recorded.""" if isinstance(source, (bytes, bytearray)): @@ -775,9 +798,7 @@ def read_sbml(source) -> SBMLModel: else: origin = Path(source) data = origin.read_bytes() - head = data[:4096].lower() - if b"")