Add Cook_AIChE2022 PEtab problem, fixes #229#316
Conversation
Adds the bone-remodeling model of Cook et al. (2022), "Mathematical Modeling of the Effects of Wnt-10b on Bone Metabolism", AIChE Journal 68(12):e17809 (doi:10.1002/aic.17809), as a new benchmark problem. Source model and code: https://github.com/ashleefv/Wnt10bBoneCompartment Model - Based on the published SimBiology/SBML export (single remodeling cycle). - SimBiology parameter/compartment IDs renamed to readable names; model id/name set to the problem ID; reference DOI added as an is-described-by annotation. - The multiple remodeling cycles of the original MATLAB implementation (state reset every 100 days: osteocytes reduced by 20, cell populations below 1 set to 0) are encoded as 11 time-triggered SBML events at t = 100, 200, ..., 1100. Bone volume z is continuous across the resets, so the bone-volume observable is unambiguous at the measurement times. PEtab problem - 3 conditions for the Wnt-10b fold changes -1, 5 and 50. - 4 measurements: the Bennett 2005 / 2007 BV/TV data used for the fit in the paper, taken after 4, 6, 6 and 12 remodeling cycles. - Observable: relative bone volume change (z - 100), unweighted least squares reproduced with fixed unit noise (normal). - 4 estimated parameters (beta1adj, alpha3adj, beta2adj, K); nominal values are the published final fit. - Dose-response visualization (BV/TV change vs Wnt-10b fold change). - Validated with petablint; imports and simulates with AMICI using default solver settings. Documentation and figures - README.md documents the model, data, event-based cycle encoding, nominal-parameter source and differences from the publication. - make_figures.py reproduces three paper/thesis figures via petab.v1.visualize (dose-response, bone volume over cycles, cell dynamics); the rendered figures are included. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01RJLtEaShERvSBHRvwCEGSA
- Add the Roser-Page (2014) CTLA-4Ig validation data (Wnt-10b fold change 1.8) as a new condition Wnt_1_8 with two BV/TV endpoints (26.6 % at 600 d / 6 cycles, 36.6 % at 1200 d / 12 cycles). These were held out from the fit in the original study; they carry their reported SDs (19.2, 40.6) as noise so they contribute negligibly to the objective. - Update the measurement, condition, simulated-data and visualization tables accordingly; the dose-response visualization now also overlays the validation points. - Extend make_figures.py to reproduce, using the publication's figure numbering, Figure 5 (validation vs Roser-Page: bone volume over 12 cycles at fold change 1.8 with the 1.2-2.4 envelope) and Figure 6 (cell-population response to Wnt-10b: max pre-osteoblast/osteoblast/osteoclast counts vs fold change). - Document the validation data and the two new figures in README.md. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Match the numbering of the publication: - Figure 4: Roser-Page validation (renamed from fig5_validation_roserpage). - Figure 5: activated cell-population dynamics over a single remodeling cycle for Wnt-10b fold changes -1, 5, 50 (osteocytes ~unchanged, pre-osteoblasts slightly up, osteoblasts up, osteoclasts down), matching the paper caption. - Figure 6: pre-osteoblast:osteoblast and osteoclast:osteoblast AUC ratios vs Wnt-10b fold change (replaces the earlier max-cell-count figure). Regenerate figures, use a numpy>=2 compatible trapezoid, and update README. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Convert the problem to PEtab v2 (format_version 2.0.0) and move the remodeling-cycle resets out of the SBML model into the PEtab experiment/condition tables: - Model: remove the 11 time-triggered SBML events; the model is now a plain event-free ODE system. - Conditions (conditions_Cook_AIChE2022.tsv): Wnt-10b dose conditions plus a cycle_reset condition (Osteocytes__S -> S - 20; P/B/C set to 0 if < 1, via piecewise), expressed with self-referential target values. - Experiments (experiments_Cook_AIChE2022.tsv): one timecourse per dose that applies cycle_reset at each 100-day boundary up to the required cycles. - Measurements now reference experimentId; observable table drops observableTransformation; parameter table drops parameterScale (v2). - Drop the v1 condition and visualization tables (v2 config has no visualization files). The period-wise experiment formulation is numerically identical to the former SBML-event encoding (max |diff| ~1e-7 at tight solver tolerances). Rewrite make_figures.py to read the v2 problem and simulate the event-free model cycle-by-cycle (no petab.v1.visualize); regenerate all figures and the simulated-data table. Update README for the v2 formulation. Validated with petablint and petab.v2 (no issues). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Use AMICI (the engine behind src/python/simulate.py) as the simulator: - Regenerate simulatedData with AMICI via the PEtab v2 importer/simulator (amici PetabImporter/PetabSimulator), which encodes the experiment periods as events and applies the cycle_reset condition natively. Total log-likelihood at the nominal parameters is ~-152.69. - Rewrite make_figures.py to compile the event-free SBML with AMICI and integrate it cycle-by-cycle (AMICI as the ODE integrator), dropping the libroadrunner dependency. Trajectories at fold changes outside the problem (the 1.2-2.4 validation envelope, the Figure 6 sweep) use the same routine. - Loosen the AMICI absolute tolerance from the default 1e-16 to 1e-12: with the default the solver reports a too-small step right after the cycle reset. - Document the AMICI simulation and the tolerance requirement in the README. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
make_figures.py now simulates the PEtab v2 problem itself: it imports the problem with AMICI's PEtab importer and simulates each of the four dose experiments with dense output, so AMICI applies the cycle_reset condition at each 100-day boundary. All figures are derived from those experiments instead of a hand-rolled per-cycle integration. Consequently two figures are reduced (as noted in the README): Figure 4 drops the 1.2-2.4 envelope and Figure 6 shows the AUC ratios at the four available fold changes (-1, 1.8, 5, 50) as discrete points, since those extra fold changes are not experiments in the problem. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
The collection tooling loaded every problem with the PEtab v1 reader, which fails on a v2 problem YAML and broke `bmp-create-overview`, `bmp-check-sbml-metadata` and the pytest suite. - base.py: `get_problem` now dispatches on the YAML `format_version` and returns a `petab.v1` or `petab.v2` problem accordingly. - overview.py: add a version-tolerant `get_sbml_model` helper; count conditions via experiments for v2; count estimated parameters via `x_free_ids`; read model files from either the v1 (`problems`/`sbml_files`) or v2 (`model_files`) YAML layout. - check_sbml_metadata.py: use the shared `get_sbml_model` helper. - Refresh the README overview row for Cook_AIChE2022 and drop a stale note from its README. petablint, overview, sbml-metadata checks and pytest all pass for the full collection (v1 and v2 problems). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Move the AMICI PEtab-v2 simulation code out of the Cook_AIChE2022 example and into src/python/simulate.py, and add an explicit v1/v2 switch: - `simulate(problem)` dispatches on the problem type: v1 problems use `import_petab_problem` + `simulate_petab`; v2 problems use the AMICI PEtab importer (`PetabImporter.create_simulator`), which encodes the experiment periods as events. - `create_v2_simulator(problem)` is a reusable helper (loosens AMICI's very tight default absolute tolerance, which otherwise fails at experiment-period state resets). Cook_AIChE2022/make_figures.py now imports it instead of building its own importer. - simulate() reports per-condition/experiment llh, chi2 and status and the total log-likelihood for both formats. Update the Simulate workflow to install the released `amici[petab]` (plus the swig/BLAS/HDF5 build dependencies and the collection package) instead of the pinned development build. Verified locally: `simulate.py Boehm_JProteomeRes2014` (v1) and `simulate.py Cook_AIChE2022` (v2) both run and report likelihoods; the Cook figures still reproduce. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
AMICI 1.0 is released, so target its public import paths and stop catering to the development build: - simulate.py imports `PetabImporter` from `amici.importers.petab` and the v1 helpers from `amici.sim.sundials.petab.v1`, dropping the try/except fallback. - Declare the simulation dependency as an unpinned `simulate` extra (`amici[petab]`) in pyproject; the Simulate workflow installs `src/python[simulate]`. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
Reproduction figures (Cook_AIChE2022)Generated with Figure 1 — dose-responseSimulation (o) vs Bennett 2005/2007 fitting data and Roser-Page 2014 validation data (x). Figure 2 — relative bone volume over remodeling cyclesSawtooth trajectory for Wnt-10b fold changes -1, 5 and 50, with literature endpoints. Figure 3 — cell-population dynamics (fold change 50)Figure 4 — validation vs Roser-Page 2014Bone volume at fold change 1.8 with the two held-out Roser-Page endpoints. Figure 5 — cell dynamics over a single cycleFigure 6 — cell-population AUC ratios vs Wnt-10b |
The collection tooling now imports `petab.v2` (in `get_problem`, and via `petablint` for the v2 problem), and `petab.v2` imports `petab.v1.distributions` which requires scipy. scipy was not a declared dependency, so a clean install (as in CI) failed with `ModuleNotFoundError: No module named 'scipy'` on the petablint / overview / metadata / pytest steps once a v2 problem was present. Verified in a fresh Python 3.12 venv (`pip install -e src/python`): petablint (v1 and v2), bmp-create-overview, bmp-check-sbml-metadata and pytest all pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014d2j46C3U5ft3yP6xrk1XZ
There was a problem hiding this comment.
Pull request overview
Adds a new PEtab v2 benchmark problem (Cook_AIChE2022) for the Cook et al. (2022) bone-remodeling model and updates the repository’s Python tooling/CI to better support PEtab v2 problems (simulation, metadata/overview utilities).
Changes:
- Add the Cook_AIChE2022 benchmark model + PEtab v2 problem definition (SBML, PEtab tables, README, figure reproduction script, and simulated/measurement data).
- Extend
src/pythontooling to handle both PEtab v1 and v2 when loading problems, simulating with AMICI, and generating overview/metadata. - Update CI simulation workflow and Python package extras/dependencies to install AMICI and PEtab v2 requirements.
Reviewed changes
Copilot reviewed 17 out of 23 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
src/python/simulate.py |
Adds PEtab v2 simulation path via AMICI PetabImporter, and a v1/v2 dispatch simulate() helper. |
src/python/pyproject.toml |
Adds scipy dependency (for PEtab v2 import chain) and a simulate extra for AMICI. |
src/python/benchmark_models_petab/overview.py |
Makes overview stats and SBML4Humans links work for PEtab v1 and v2 layouts. |
src/python/benchmark_models_petab/check_sbml_metadata.py |
Uses a v1/v2 SBML accessor to validate model metadata across formats. |
src/python/benchmark_models_petab/base.py |
Loads PEtab problems as v1 or v2 based on format_version; keeps simulation TSV loading via v1 helper. |
README.md |
Adds the Cook_AIChE2022 row to the repository-wide benchmark overview table. |
Benchmark-Models/Cook_AIChE2022/README.md |
Documents the new benchmark problem, data sources, v2 structure, cycle-reset encoding, and figure reproduction. |
Benchmark-Models/Cook_AIChE2022/model_Cook_AIChE2022.xml |
Adds the SBML model (no SBML events; intended to be used with PEtab v2 experiment-period resets). |
Benchmark-Models/Cook_AIChE2022/Cook_AIChE2022.yaml |
Declares the PEtab v2 problem and associated tables/files. |
Benchmark-Models/Cook_AIChE2022/conditions_Cook_AIChE2022.tsv |
Defines dose conditions and the cycle_reset state-change condition. |
Benchmark-Models/Cook_AIChE2022/experiments_Cook_AIChE2022.tsv |
Encodes repeated cycles via experiment periods applying cycle_reset at 100-day boundaries. |
Benchmark-Models/Cook_AIChE2022/measurementData_Cook_AIChE2022.tsv |
Adds fitting + validation measurements. |
Benchmark-Models/Cook_AIChE2022/simulatedData_Cook_AIChE2022.tsv |
Adds nominal-parameter simulated outputs for the measurements. |
Benchmark-Models/Cook_AIChE2022/parameters_Cook_AIChE2022.tsv |
Defines estimated parameters, bounds, and nominal values. |
Benchmark-Models/Cook_AIChE2022/observables_Cook_AIChE2022.tsv |
Defines the bone-volume observable and its noise placeholder wiring. |
Benchmark-Models/Cook_AIChE2022/make_figures.py |
Reproduces paper/thesis figures by simulating the PEtab v2 problem with AMICI. |
.github/workflows/simulate.yml |
Installs AMICI via the project’s simulate extra and adds required system deps for building AMICI. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
dilpath
left a comment
There was a problem hiding this comment.
Are you able to get a similarly-good or better fit from fitting this PEtab problem, as the original paper?
Since this is a bit different to the normal submission, I'll wait for feedback on some of my comments before approving.
There was a problem hiding this comment.
Regarding a few of these package files, it's a bit messy with the combined v1/v2 functions but not sure what else to do. Does it make sense to instead maintain one overview table for v1 problems, and a separate overview table for v2 problems?
There was a problem hiding this comment.
I think we first should decide on how we want to maintain this collection. Have a distinct set of v1 and v2 problems, or have them implemented in both versions? Due to the lack of resources I would just upconvert the old problems, but not backport anything. In that case, having different tables and modularizing the code makes more sense.
| * **Estimated parameters** (`parameters`): the four Wnt-10b-related parameters | ||
| `beta1adj`, `alpha3adj`, `beta2adj`, `K`. (PEtab v2 has no `parameterScale` | ||
| column; the bounds are given on linear scale.) | ||
|
|
There was a problem hiding this comment.
@dweindl do you think these parameters should rather be exp(parameter) in the model, or the log handled elsewhere, in some tool-specific manner? In general, for parameterScale functionality.
There was a problem hiding this comment.
In the original study, no log-scaling was used as far as I understand. It's the question how much close we want to be to the original. Not sure.
There was a problem hiding this comment.
Until now, this kind of information belongs in the GitHub repo issue. Fine to keep here but then add a link to this from the github issue.
There was a problem hiding this comment.
I think it's good to have it in the repository, this way, it will also be included in the zenodo archives.
There was a problem hiding this comment.
I think also the scripts are helpful for better understand what has been done initially. For other problems, this turned out to be very difficult... However, I would clearly separate the PEtab problem from the rest. The Python package should ideally only include the PEtab files, not all the other import/conversion/validation scripts.
| exp_Wnt_5 obs_BV 69.2 600 1 Bennett2007_6cycles_Wnt_5 | ||
| exp_Wnt_50 obs_BV 339.0 1200 1 Bennett2005_12cycles | ||
| exp_Wnt_1_8 obs_BV 26.6 600 19.2 RoserPage2014_6cycles | ||
| exp_Wnt_1_8 obs_BV 36.6 1200 40.6 RoserPage2014_12cycles |
There was a problem hiding this comment.
If you move the measurement data out, I would suggest moving the 1 noiseParameters to the training data observable table to simplify the specification of the training problem overall.
I'm also confused why standard deviation is considered for the validation data, but ignored for the training data.
| <parameter id="g_31" name="g_31" value="1" constant="true"/> | ||
| <parameter id="g_21" name="g_21" value="2" constant="true"/> | ||
| <parameter id="g_22" name="g_22" value="1" constant="true"/> | ||
| <parameter id="g_32" name="g_32" value="1" constant="true"/> | ||
| <parameter id="g_41" name="g_41" value="1" constant="true"/> | ||
| <parameter id="g_42" name="g_42" value="1" constant="true"/> | ||
| <parameter id="g_43" name="g_43" value="-1" constant="true"/> | ||
| <parameter id="g_44" name="g_44" value="1" constant="true"/> | ||
| <parameter id="f_12" name="f_12" value="1" constant="true"/> | ||
| <parameter id="f_14" name="f_14" value="1" constant="true"/> | ||
| <parameter id="f_23" name="f_23" value="1" constant="true"/> | ||
| <parameter id="f_34" name="f_34" value="1" constant="true"/> |
There was a problem hiding this comment.
Should these be parameters?
| def _relax_tolerances(solver) -> None: | ||
| """Loosen AMICI's very tight default absolute tolerance. | ||
|
|
||
| The default ``atol`` (1e-16) can trigger a too-small step size right after | ||
| a PEtab experiment-period state change; 1e-12 is robust while remaining | ||
| accurate. | ||
| """ | ||
| solver.set_relative_tolerance(1e-10) | ||
| solver.set_absolute_tolerance(1e-12) | ||
| solver.set_max_steps(10**6) |
There was a problem hiding this comment.
This seems like problem-specific information. Does it substantially change any of the computed LLHs?
Co-authored-by: Dilan Pathirana <59329744+dilpath@users.noreply.github.com>
- base.py: use petab.versions.get_major_version for the v1/v2 dispatch (also handles non-numeric versions such as the beta_* problems) - Split the Roser-Page validation data into a separate measurement table that is not referenced by the problem YAML, so the fitted problem holds only the four Bennett points; move the unit noise into the observable (noiseFormula = 1) and drop the per-measurement noiseParameters - Reparameterize to estimate k2adj (= alpha3adj - beta1adj, lower bound > 0) and derive alpha3adj = beta1adj + k2adj via the condition table, so the original's alpha3adj > beta1adj constraint holds during fitting - simulate.py: make the loosened AMICI solver tolerances opt-in per problem instead of forcing them on every PEtab v2 problem - make_figures.py: read validation points and SDs from the validation table, use np.trapezoid, and pass the solver settings explicitly - README: document the reparameterization/constraint, add the clamp fitting caveat, and fix the parameterization wording Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Dn5uzjifxYJnEmqEKh2cR8
| - name: Install system dependencies | ||
| run: | | ||
| sudo apt-get update | ||
| sudo apt-get install -y swig libopenblas-dev libhdf5-dev | ||
|
|
There was a problem hiding this comment.
I think it's good to have it in the repository, this way, it will also be included in the zenodo archives.
| beta1adj $\beta_1^{\mathrm{adj}}$ 0.0001 1.0 0.177617716487146 true | ||
| k2adj $k_2$ 0.0001 1.1 0.083313316273387 true | ||
| beta2adj $\beta_2^{\mathrm{adj}}$ 1e-07 1.0 0.000709650034656732 true | ||
| K $K$ 1.0 10.0 6.26349707992014 true |
There was a problem hiding this comment.
Isn't that supposed to be [1, 100]? https://github.com/ashleefv/Wnt10bBoneCompartment/blob/5a8a7da42968e713516eda82394ab1512a722bba/BoneCompartmentUp4Para.m#L50
There was a problem hiding this comment.
Ah, there are different configurations below. But the paper says
We assume that the saturation parameter KM could take on a value between 1 and 100, and we bound the parameter fitting accordingly.
| * **Estimated parameters** (`parameters`): the four Wnt-10b-related parameters | ||
| `beta1adj`, `alpha3adj`, `beta2adj`, `K`. (PEtab v2 has no `parameterScale` | ||
| column; the bounds are given on linear scale.) | ||
|
|
There was a problem hiding this comment.
In the original study, no log-scaling was used as far as I understand. It's the question how much close we want to be to the original. Not sure.
There was a problem hiding this comment.
I think we first should decide on how we want to maintain this collection. Have a distinct set of v1 and v2 problems, or have them implemented in both versions? Due to the lack of resources I would just upconvert the old problems, but not backport anything. In that case, having different tables and modularizing the code makes more sense.
There was a problem hiding this comment.
I think also the scripts are helpful for better understand what has been done initially. For other problems, this turned out to be very difficult... However, I would clearly separate the PEtab problem from the rest. The Python package should ideally only include the PEtab files, not all the other import/conversion/validation scripts.






Adds the bone-remodeling model of Cook et al. (2022), "Mathematical Modeling of the Effects of Wnt-10b on Bone Metabolism", AIChE Journal 68(12):e17809 (doi:10.1002/aic.17809), as a new benchmark problem.
Source model and code: https://github.com/ashleefv/Wnt10bBoneCompartment
Model
PEtab problem
Documentation and figures
Claude-Session: https://claude.ai/code/session_01RJLtEaShERvSBHRvwCEGSA
Checklist for the submission of new PEtab problems
{LAST_NAME_OF_FIRST_AUTHOR}_{ABBREVIATED_JOURNAL_NAME}{YEAR_OF_PUBLICATION}petablint -vy problem.yaml)parameterNamein the parameter table orxLabelin the visualization table) encode special characters with tex-style strings instead of unicode, e.g. useH$_2$Ofor water.bmp-create-overview --update(requirespip install -e src/pythonfrom the repository root)