The gap
For fitting ordinary differential equation models there is an accepted set of benchmark problems, the PEtab benchmark collection of Hass and colleagues. Anyone proposing a new fitting method for those models is expected to run it on that collection and report the result.
Nothing equivalent exists for fitting stochastic, rule-based models. There is no agreed set of problems, no reference answers, and no scoring protocol. So no one can say whether one method fits these models better than another, including us.
We are better placed to build one than anyone else, and we need it for our own work regardless. Issues #659, #660, #661, and #662 all claim that some change will produce better fits, and not one of those claims is measurable today.
Why a stochastic benchmark can be better than the ODE one
The PEtab collection scores a method against the best answer anyone has found so far, because the true parameters of a real published model are unknown.
For a stochastic model we can do better. Pick a model, choose parameter values, simulate it to produce data, and then ask whether a fitting run recovers the values we chose. The right answer is known exactly, because we chose it. That makes the score meaningful in a way "did you match the best previous attempt" is not.
It also means we are not limited by how many published fitting problems exist. Any model we can simulate becomes a benchmark problem. What is scarce is careful selection and a scoring protocol, not models.
What we already have
The machinery exists in tests/recovery_harness.py. Its simulate_truth function already takes a model and a set of true parameter values, simulates it, and writes out a data file. tests/test_recovery.py already drives fits against that data and checks whether the true values come back. There are seven models in tests/recovery_models/. This is a working parameter recovery benchmark. It is small, it is aimed mostly at differential equation models, and it is wired up as a test tier rather than as something publishable.
For models, we have a great deal to draw on:
| Source |
BNGL files |
Declaring a stochastic method |
bngsim/parity_checks/bng_parity |
895 |
341 |
RuleMonkey |
222 |
194 |
BNGL-Models/models |
129 |
49 |
bngsim/benchmarks |
113 |
26 |
BNGL_library/bngl_models |
56 |
55 |
BNGL-Models/pybnf-jobs |
81, plus 102 conf and 327 exp files |
0 |
PyBNF/tests/recovery_models |
7 |
3 |
Two notes on that table. The parity check models exist to confirm that bngsim reproduces BNG2.pl, so many of them are likely small and synthetic rather than real published models, and they will need filtering. The pybnf-jobs collection is the opposite. It is entirely differential equation work, but it is the only source that contains complete fitting jobs with configuration files and data, so it is the best template for what a finished benchmark entry should look like.
The curated published models in BNGL-Models/models and BNGL_library/bngl_models are the most promising starting pool, because a benchmark built from real published biology is far more persuasive than one built from toy models.
What to build
-
Choose the problems. Aim for something like twenty to thirty stochastic models, chosen deliberately for a range of size, number of free parameters, and how noisy their output is. Prefer published models with a citation. Quality matters much more than quantity here.
-
Fix the problem definitions. For each model, record the true parameter values, the simulation settings used to generate the data, the observables measured, the sampling times, the number of replicates, and the search bounds. Freeze and version these so that a result reported today can be compared with one reported in two years.
-
Write down the scoring rules. At minimum: how close the recovered parameters are to the true ones, how often a run gets there at all across repeated attempts, and how many simulations it took. Repeated attempts matter more here than for differential equation models, since both the fitting method and the simulation are random.
-
Produce baseline results. Run our own methods across the suite. Then, and this is the part that decides whether anyone believes it, run methods that are not ours. bngsim/parity_checks already contains working setups for AMICI, COPASI and RoadRunner, which is an unusual head start on exactly this.
-
Publish it. A public repository with the models, the generated data, the protocol and the baseline results, given a permanent identifier so it can be cited. See the note below.
Where it should live
The work starts here, because the harness is here. The finished suite probably should not live in the PyBNF repository, since a benchmark that lives inside one of the tools it scores is easy to dismiss. BNGL-Models or a new repository would be a better long term home. Worth deciding early, because it affects how the harness is factored.
Relationship to other issues
This is what makes #659, #660, #661 and #662 measurable instead of asserted. #659 in particular has a direct connection, since the size of the noise bias it describes is exactly the kind of thing this suite would quantify across many models rather than on one example.
Doing a small version of this first, with perhaps five models, would be enough to start scoring the algorithm work while the full suite is assembled.
The gap
For fitting ordinary differential equation models there is an accepted set of benchmark problems, the PEtab benchmark collection of Hass and colleagues. Anyone proposing a new fitting method for those models is expected to run it on that collection and report the result.
Nothing equivalent exists for fitting stochastic, rule-based models. There is no agreed set of problems, no reference answers, and no scoring protocol. So no one can say whether one method fits these models better than another, including us.
We are better placed to build one than anyone else, and we need it for our own work regardless. Issues #659, #660, #661, and #662 all claim that some change will produce better fits, and not one of those claims is measurable today.
Why a stochastic benchmark can be better than the ODE one
The PEtab collection scores a method against the best answer anyone has found so far, because the true parameters of a real published model are unknown.
For a stochastic model we can do better. Pick a model, choose parameter values, simulate it to produce data, and then ask whether a fitting run recovers the values we chose. The right answer is known exactly, because we chose it. That makes the score meaningful in a way "did you match the best previous attempt" is not.
It also means we are not limited by how many published fitting problems exist. Any model we can simulate becomes a benchmark problem. What is scarce is careful selection and a scoring protocol, not models.
What we already have
The machinery exists in
tests/recovery_harness.py. Itssimulate_truthfunction already takes a model and a set of true parameter values, simulates it, and writes out a data file.tests/test_recovery.pyalready drives fits against that data and checks whether the true values come back. There are seven models intests/recovery_models/. This is a working parameter recovery benchmark. It is small, it is aimed mostly at differential equation models, and it is wired up as a test tier rather than as something publishable.For models, we have a great deal to draw on:
bngsim/parity_checks/bng_parityRuleMonkeyBNGL-Models/modelsbngsim/benchmarksBNGL_library/bngl_modelsBNGL-Models/pybnf-jobsPyBNF/tests/recovery_modelsTwo notes on that table. The parity check models exist to confirm that bngsim reproduces BNG2.pl, so many of them are likely small and synthetic rather than real published models, and they will need filtering. The
pybnf-jobscollection is the opposite. It is entirely differential equation work, but it is the only source that contains complete fitting jobs with configuration files and data, so it is the best template for what a finished benchmark entry should look like.The curated published models in
BNGL-Models/modelsandBNGL_library/bngl_modelsare the most promising starting pool, because a benchmark built from real published biology is far more persuasive than one built from toy models.What to build
Choose the problems. Aim for something like twenty to thirty stochastic models, chosen deliberately for a range of size, number of free parameters, and how noisy their output is. Prefer published models with a citation. Quality matters much more than quantity here.
Fix the problem definitions. For each model, record the true parameter values, the simulation settings used to generate the data, the observables measured, the sampling times, the number of replicates, and the search bounds. Freeze and version these so that a result reported today can be compared with one reported in two years.
Write down the scoring rules. At minimum: how close the recovered parameters are to the true ones, how often a run gets there at all across repeated attempts, and how many simulations it took. Repeated attempts matter more here than for differential equation models, since both the fitting method and the simulation are random.
Produce baseline results. Run our own methods across the suite. Then, and this is the part that decides whether anyone believes it, run methods that are not ours.
bngsim/parity_checksalready contains working setups for AMICI, COPASI and RoadRunner, which is an unusual head start on exactly this.Publish it. A public repository with the models, the generated data, the protocol and the baseline results, given a permanent identifier so it can be cited. See the note below.
Where it should live
The work starts here, because the harness is here. The finished suite probably should not live in the PyBNF repository, since a benchmark that lives inside one of the tools it scores is easy to dismiss.
BNGL-Modelsor a new repository would be a better long term home. Worth deciding early, because it affects how the harness is factored.Relationship to other issues
This is what makes #659, #660, #661 and #662 measurable instead of asserted. #659 in particular has a direct connection, since the size of the noise bias it describes is exactly the kind of thing this suite would quantify across many models rather than on one example.
Doing a small version of this first, with perhaps five models, would be enough to start scoring the algorithm work while the full suite is assembled.