benchmarks.md

Benchmarks & validation

BNGsim's SSA (Gillespie direct method with dependency graph + 4-ary sum tree) and PSA (partial-scaling approximation) were benchmarked against run_network 3.0 (BioNetGen 2.9.3) across 10 SSA models and 3 PSA models (6 configurations). BNGsim timing covers sim.run() only; run_network timing includes full subprocess overhead. Protocol: 2 warmup + 5 timed runs, median reported.

SSA results (BNGsim wins every model)

Model	Species	Reactions	BNGsim	run_network	Speedup
gene_expression_hill	2	4	0.0001 s	0.046 s	547×
simple_system	4	4	0.002 s	0.045 s	20×
flagellar_motor	4	4	0.0008 s	0.045 s	54×
gene_expr_3stage	6	6	0.50 s	1.02 s	2.1×
oscillatory_system	5	8	0.25 s	0.34 s	1.3×
gene_expression	10	14	0.004 s	0.053 s	15×
tcr_signaling	37	97	0.68 s	1.69 s	2.5×
egfr_net	356	3749	0.086 s	0.44 s	5.2×
multisite_phos	1026	7680	0.36 s	1.75 s	4.9×
fceri_gamma	3744	58276	2.24 s	64.5 s	28.8×

Geometric mean speedup: 8.7×

PSA results (BNGsim wins every configuration)

Model	Nc	BNGsim	run_network	Speedup
tcr_signaling	10	0.015 s	0.098 s	6.5×
tcr_signaling	100	0.039 s	0.173 s	4.4×
prion_aggregation	10	0.46 s	1.33 s	2.9×
prion_aggregation	100	0.45 s	1.35 s	3.0×
erk_activation	10	0.14 s	6.38 s	44.8×
erk_activation	100	1.70 s	10.9 s	6.4×

Geometric mean PSA speedup: 6.8×

BNGsim's advantage comes from: (1) zero subprocess overhead (in-process execution), (2) O(log N) sum-tree reaction selection (vs O(N) linear scan in run_network), and (3) pre-computed propensity and dependency-graph data — including the per-reaction affected-set used after each fire — with zero heap allocation, sort, or dedup in the SSA hot loop. The advantage grows with network size — 28.8× on fceri_gamma (58K reactions).

PSA on SBML / Antimony models

Simulator(model, method="psa", poplevel=N_c) accepts SBML or Antimony models loaded via Model.from_sbml(...) / Model.from_antimony(...) — the same dispatch as .net, sharing the validate_for_ssa gate at construct time (reversible_non_mass_action, non_integer_stoichiometry, etc. raise SsaValidationError before the simulator is built). PSA semantics for N_c are identical to the .net path; choose per Lin, Feng, Hlavacek (2019) — accuracy degrades when N_c is too small relative to per-species populations (e.g. erk_activation requires N_c ≳ 100 to track the network reference).

PSA's leaping advantage carries fully through the SBML loader. On tcr_signaling loaded via Model.from_sbml(...) (3 timed runs, median; t_end=300, n_steps=1000), PSA beats SSA by 8.6× at N_c=300, 15.3× at N_c=100, 34.8× at N_c=30, and 46.7× at N_c=10, matching the SSA→PSA speedup ranges measured on the .net path.

ODE parity vs RoadRunner (results ship in-repo)

bngsim's SBML ODE engine is cross-checked against libRoadRunner over the full BioModels ODE corpus (1,323 models) by the parity_checks/rr_parity suite: both engines integrate each model at a shared tolerance and their trajectories are compared cell-by-cell. The latest released result is committed in this repository (so it's there the moment you clone or install) as a gzipped report under parity_checks/rr_parity/snapshots/ — see that directory's README for the convention. Regenerate the interactive HTML matrix (per-model verdicts, three-tier timing breakdown, BNGsim-vs-RoadRunner cost scatter, and win-by-species tables) from a snapshot with:

cd parity_checks/rr_parity
gunzip -c snapshots/report_ode_<version>.json.gz > runs/report_ode.json
python generate_matrix.py            # writes runs/parity_matrix.html

As of the latest snapshot: 1,237/1,323 models agree to ~1e-4 relative error; on warm per-integration cost BNGsim is ~1.65× faster by geometric mean and wins the majority of models, with RoadRunner pulling ahead on the high-species tail (the crossover is around 20–35 species).

SBML semantic test suite — capability boundary and intentional deviations

Against the SBML semantic test suite (v3.3.0, 1823 cases), graded through a faithful local port of the official test-runner grading, bngsim scores 1577 Match (1577/1789 in-scope time-course cases = 88.1%). The full breakdown is 1577 Match, 242 Unsupported (the declared capability boundary below), 3 NoMatch (the intentional deviation below), 1 CannotSolve, and 0 Error — no case is ever a wrong-but-plausible answer. In a shared-tolerance head-to-head bngsim also passes more cases than RoadRunner (1578 vs 1529 under the in-repo fair harness). The cases bngsim does not match are all deliberate — a declared capability boundary and one intentional deviation, not open bugs.

That score is reproducible and reviewable: benchmarks/suites/sbml_test_suite/testrunner/ ships a fail-closed wrapper for the official runner, a committed unsupported-tags manifest (bngsim-unsupported-tags.txt, pinned by a unit test to the shipped source of truth bngsim._sbml_unsupported), and score.py — a line-for-line port of the runner's comparison, outcome enum, and tag matching. See that directory's README for the runbook and the full per-category breakdown.

Unsupported constructs (refused at load — never silently approximated). bngsim is an ODE/SSA kinetic engine with no differential-algebraic, delay, or fast-equilibrium solver. Where such a construct appears in the equations that define the system, bngsim fails closed with a ModelError naming it — rather than integrate a different mathematical system and return a confident-but-wrong trajectory (opt into the legacy silent approximation with BNGSIM_ALLOW_UNSUPPORTED_CONSTRUCTS=1):

Construct (tag)	Why unsupported	Trivial form still simulated
AlgebraicRule (non-empty)	DAE constraint — no algebraic/DAE solver	an empty <algebraicRule/> (no constraint)
csymbol delay / delay(x, τ≠0)	delay-differential equation — no DDE solver	delay(x, 0) = identity
FastReaction (fast="true") under ODE	fast-equilibrium constraint — no constraint solver	—

A whole SBML Level 3 package declared required="true" is refused the same way and for the same reason: required="true" means the package changes the model's mathematical meaning, so reading only the core layer of a multi, qual or spatial document integrates a different system than the one on disk. comp (flattened) and distrib (random draws refused individually by the math translator) are the two bngsim accounts for; a presentation-only package declares required="false" and loads untouched. No SBML Test Suite case is affected — the suite declares only comp:required="true" and fbc:required="false".

The suite's flux-balance test type (fbc / FluxBalanceSteadyState, 34 cases) is likewise declared unsupported — it is a steady-state constraint problem, not a time-course, and bngsim is a time-course ODE/SSA engine. All of the above are recorded in the committed manifest, so the official grading reports them as Unsupported (an honest capability boundary) rather than as failures.

This is a genuine capability boundary shared with RoadRunner, which declares the same constructs unsupported in the SBML Test Suite Database. bngsim's declared-unsupported set is narrower, though: it fully or partially supports several tags others declare unsupported — e.g. VolumeConcentrationRates (time-varying compartment-volume dilution, 15/19 suite cases) and AssignedVariableStoichiometry (named / rate-ruled <speciesReference> stoichiometry, 27/73 — ahead of RoadRunner's 20).

Intentional deviations (simulated, but deliberately differ). Three cases where bngsim runs the model but does not match the suite's reference on purpose:

• Avogadro constant (3 cases: 00960, 00961, 01323). bngsim uses the SI-exact Avogadro number Nₐ = 6.02214076e23 (fixed by the 2019 redefinition of the SI mole) everywhere the avogadro csymbol appears. The suite's analytic reference predates that redefinition and uses the older measured value, so bngsim differs in the last few digits. We match the current SI definition on purpose.

Reproducible random event ordering. When several events fire at the same instant with equal (or unset) priority, SBML L3v2 §4.11.6 chooses their order at random. bngsim implements this with a per-run seed (SolverOptions.event_seed, exposed as Simulator(method="ode").run(seed=…)): equal-priority ties are broken by a std::mt19937_64, so the RandomEventExecution suite family (00962, 01588, 01590, 01591, 01599, 01605, 01627 and the T0/no-delay variants) passes, while a fixed default seed keeps every run reproducible out of the box — the hard requirement for the parameter-fitting workloads bngsim is built for (PyBNF). The RNG is consumed only at a genuine equal-priority tie, so any model without such simultaneity is byte-identical regardless of the seed, and an event-free ODE run stays seed-less (Result.seed is None). Pass an explicit seed for an independent event-ordering realization.

All other suite gaps are tracked as ordinary compliance work.