diff --git a/.github/workflows/simulate.yml b/.github/workflows/simulate.yml index 1a228c4..87125bf 100644 --- a/.github/workflows/simulate.yml +++ b/.github/workflows/simulate.yml @@ -20,15 +20,17 @@ jobs: with: python-version: "3.13" + - name: Install system dependencies + run: | + sudo apt-get update + sudo apt-get install -y swig libopenblas-dev libhdf5-dev + - name: Install AMICI env: AMICI_PARALLEL_COMPILE: "" run: | pip install --upgrade pip - # TODO: once amici 1.0 is released, replace with: - # pip install amici[petab] - pip install petab - pip install "git+https://github.com/dweindl/amici.git@6a85422ad233fc268fb71272bba39f5b82db387a#egg=amici&subdirectory=python/sdist" + pip install -e "src/python[simulate]" - name: Simulate with nominal parameters env: diff --git a/Benchmark-Models/Cook_AIChE2022/Cook_AIChE2022.yaml b/Benchmark-Models/Cook_AIChE2022/Cook_AIChE2022.yaml new file mode 100644 index 0000000..94e11c3 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/Cook_AIChE2022.yaml @@ -0,0 +1,17 @@ +format_version: 2.0.0 +parameter_files: +- parameters_Cook_AIChE2022.tsv +model_files: + model_Cook_AIChE2022: + location: model_Cook_AIChE2022.xml + language: sbml +measurement_files: +- measurementData_Cook_AIChE2022.tsv +condition_files: +- conditions_Cook_AIChE2022.tsv +experiment_files: +- experiments_Cook_AIChE2022.tsv +observable_files: +- observables_Cook_AIChE2022.tsv +mapping_files: [] +extensions: {} diff --git a/Benchmark-Models/Cook_AIChE2022/README.md b/Benchmark-Models/Cook_AIChE2022/README.md new file mode 100644 index 0000000..b89adc5 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/README.md @@ -0,0 +1,223 @@ +# Cook_AIChE2022 + +PEtab implementation of the bone-remodeling model of + +> C. V. Cook, M. A. Islam, B. J. Smith, A. N. Ford Versypt. +> *Mathematical Modeling of the Effects of Wnt-10b on Bone Metabolism.* +> AIChE Journal, 2022, 68(12), e17809. doi:[10.1002/aic.17809](https://doi.org/10.1002/aic.17809) + +Original model, data and MATLAB code: https://github.com/ashleefv/Wnt10bBoneCompartment +Additional background: C. V. Cook, *Mathematical Modeling to Connect Bone +Responses to Systemic Mechanisms* (PhD thesis). + +This is a **PEtab v2** problem (`format_version: 2.0.0`). The repeated +remodeling cycles are encoded in the PEtab experiment/condition tables rather +than in the SBML model (see below). + +## Model + +A compartmental ODE model of a bone multicellular unit over repeated +remodeling cycles. Five species are tracked: osteocytes (`S`), +pre-osteoblasts (`P`), osteoblasts (`B`), osteoclasts (`C`) and relative +bone volume (`z`, in % of the normal value, starting at 100). Wnt-10b enters +the system through a Hill-type factor `piwnta = Wnt / (Wnt + K)` that modulates +pre-osteoblast proliferation, pre-osteoblast-to-osteoblast differentiation and +osteoblast apoptosis. + +## Data + +Four bone-volume-fraction (BV/TV) measurements compiled in the paper from +Wnt-10b transgenic / knockout mouse experiments (Bennett et al. 2005, 2007), +expressed as percentage change relative to normal Wnt-10b: + +| condition | Wnt-10b fold change | remodeling cycles | time [d] | BV/TV change [%] | source | +|-----------|--------------------:|------------------:|---------:|-----------------:|--------| +| `Wnt_m1` | -1 | 4 | 400 | -29.7 | Bennett 2005 | +| `Wnt_m1` | -1 | 6 | 600 | -41.9 | Bennett 2007 | +| `Wnt_5` | 5 | 6 | 600 | 69.2 | Bennett 2007 | +| `Wnt_50` | 50 | 12 | 1200 | 339.0 | Bennett 2005 | + +### Validation data (Roser-Page 2014) + +The paper additionally *validates* the fitted model against an independent +dataset that was **not** used for parameter estimation: the CTLA-4Ig +(abatacept) mouse experiments of Roser-Page et al. 2014, which correspond to a +Wnt-10b fold change of ~1.8 (with a 1.2–2.4 range). The two BV/TV endpoints are +reproduced here from the paper's validation figure (`ValidationResults` in the +original `GraphsforPaper.m`): + +| condition | Wnt-10b fold change | remodeling cycles | time [d] | BV/TV change [%] | SD | source | +|-----------|--------------------:|------------------:|---------:|-----------------:|----:|--------| +| `Wnt_1_8` | 1.8 | 6 | 600 | 26.6 | 19.2 | Roser-Page 2014 | +| `Wnt_1_8` | 1.8 | 12 | 1200 | 36.6 | 40.6 | Roser-Page 2014 | + +> Reference: S. Roser-Page et al., *CTLA-4Ig-induced T cell anergy promotes +> Wnt-10b production and bone formation in a mouse model*, Arthritis & +> Rheumatology, 2014, 66(4), 990–999. doi:[10.1002/art.38319](https://doi.org/10.1002/art.38319) + +These two points are kept in a **separate** measurement table +(`measurementData_validation_Cook_AIChE2022.tsv`, experiment `exp_Wnt_1_8`, +dataset ids `RoserPage2014_*`) that is **not** referenced by the problem YAML, +so they are held out of parameter estimation — matching the original study, +where Roser-Page was used only for validation. Their reported standard +deviations (19.2, 40.6) are carried in that table purely as error bars for the +validation figure; they are not used as PEtab noise. + +## PEtab problem + +PEtab v2 (`format_version: 2.0.0`). The files are: + +* **Conditions** (`conditions_Cook_AIChE2022.tsv`): the Wnt-10b dose conditions + `Wnt_m1`, `Wnt_5`, `Wnt_50`, `Wnt_1_8` (each sets the `Wnt` parameter and + `alpha3adj = beta1adj + k2adj`, see *Differences* below) and the `cycle_reset` + condition that performs the cycle-boundary state reset (below). +* **Experiments** (`experiments_Cook_AIChE2022.tsv`): one timecourse per Wnt-10b + dose (`exp_Wnt_m1`, `exp_Wnt_5`, `exp_Wnt_50`, `exp_Wnt_1_8`). Each starts at + `t = 0` with the dose condition and then applies `cycle_reset` at every 100-day + boundary (`t = 100, 200, ...`) up to the number of cycles required by that + experiment's measurements. +* **Observable** (`observables`): `obs_BV = Bone_volume__z - 100`, the relative + BV/TV change. The paper fits with unweighted least squares; this is + reproduced with a fixed unit noise (`noiseFormula = 1`, `normal`), so the + objective equals the residual sum of squares up to a constant. +* **Estimated parameters** (`parameters`): the four Wnt-10b-related parameters + `beta1adj`, `k2adj`, `beta2adj`, `K`, with `alpha3adj = beta1adj + k2adj` + derived via the condition table so that `alpha3adj > beta1adj` holds during + fitting (see *Differences from the original publication*). (PEtab v2 has no + `parameterScale` column; the bounds are given on linear scale.) + +### Multiple remodeling cycles as PEtab experiments + +The number of remodeling cycles is the key feature of the original model, and +the SBML export shipped with the paper explicitly *does not* support it. In the +MATLAB code each 100-day cycle is integrated separately and the state is reset +at the cycle boundary: cell populations that have decayed below 1 are set to 0 +and the osteocyte count is reduced by 20 to initiate the next cycle. + +Here this reset is expressed as the PEtab `cycle_reset` **condition** + +``` +cycle_reset Osteocytes__S Osteocytes__S - 20 +cycle_reset Pre_Osteoblasts__P piecewise(0, Pre_Osteoblasts__P < 1, Pre_Osteoblasts__P) +cycle_reset Osteoblasts__B piecewise(0, Osteoblasts__B < 1, Osteoblasts__B) +cycle_reset Osteoclasts__C piecewise(0, Osteoclasts__C < 1, Osteoclasts__C) +``` + +which each experiment applies at `t = 100, 200, ...` via the experiment table. +The SBML model contains **no events**; the periods of the experiment define the +cycle boundaries. Crucially, bone volume `z` is **not** reset, so it is +continuous across cycle boundaries and the `z`-based observable is unambiguous +at the measurement times. Keeping the reset in the experiment/condition tables +leaves the model itself a plain, reusable ODE system. + +### Simulation + +The problem imports and simulates with AMICI (`src/python/simulate.py`); AMICI +encodes the experiment periods as events internally. The cycle-boundary state +change makes the system stiff, so the AMICI **absolute tolerance must be +loosened from its very tight default (`1e-16`) to about `1e-12`** — with the +default `atol` the solver reports a too-small step after the reset. The +`simulatedData` table was generated with AMICI at `rtol = 1e-10`, `atol = 1e-12` +(total log-likelihood ≈ -143.38 over the four fitting measurements at the +nominal parameters). + +## Nominal parameters + +Taken from the publication's final fit (the parameter set used to generate the +paper figures in `GraphsforPaper.m`, "after the 4th data point was added"): + +| parameter | nominal | note | +|-----------|--------:|------| +| `beta1adj` | 0.177617716487146 | `= k1` | +| `k2adj` | 0.083313316273387 | `= k2`; `alpha3adj = beta1adj + k2adj = k1 + k2` | +| `beta2adj` | 0.000709650034656732 | `= k3` | +| `K` | 6.26349707992014 | `= k4` | + +## Differences from the original publication + +* **Estimated parameters and the `alpha3adj > beta1adj` constraint.** The MATLAB + code estimates `k1..k4` with `beta1adj = k1`, `beta2adj = k3`, `K = k4` and + `alpha3adj = k1 + k2`, so `alpha3adj > beta1adj` (because `k2 > 0`). To keep + that constraint during fitting, the estimated parameters here are `beta1adj`, + `k2adj` (`= k2`, lower bound > 0), `beta2adj` and `K`, and `alpha3adj` is + derived as `alpha3adj = beta1adj + k2adj` via the condition table rather than + estimated independently. This also corrects the shipped SBML export, whose + `beta1adj = 0.0833` is in fact the value of `k2` (the nominal `k2adj` here). +* **Multiple remodeling cycles** are applied through the experiment/condition + tables rather than the MATLAB code's sequential re-initialized integrations + (see *Multiple remodeling cycles as PEtab experiments* above). +* The reaction kinetics match the shipped SBML export, which omits the + `max(1 - S/K_S, 0)` clamp present in the MATLAB code. It has no effect on the + nominal trajectories (`S` stays below `K_S`); whether it matters under fitting, + where the parameters vary and `S` could exceed `K_S`, remains to be verified. + +## Fitting notes + +The nominal fit was obtained in the original study with MATLAB `lsqcurvefit` +(Levenberg-Marquardt) from Latin-hypercube / normally-distributed multistarts. +It is a least-squares compromise across the four data points (residual sum of +squares ~279 in BV/TV-% units); the `Wnt_50 / 12-cycle` point is reproduced +almost exactly. + +## Reproducing the figures + +`make_figures.py` simulates the PEtab problem itself with **AMICI**: the problem +is imported with AMICI's PEtab importer (which encodes the experiment periods as +events) and each of the four dose experiments is simulated with dense output so +that AMICI applies the `cycle_reset` condition at every 100-day boundary. Every +figure is derived from those four experiments. + +```bash +python make_figures.py # requires petab (v2), amici, matplotlib, numpy>=2 +``` + +Because the problem only contains the four fold changes -1, 1.8, 5 and 50, two +figures are reduced relative to the publication: Figure 4 omits the 1.2–2.4 +envelope, and Figure 6 shows the AUC ratios at the four available fold changes +as discrete points rather than a continuous sweep. + +### Figure 1 — dose-response +Simulation (filled circles) vs Bennett 2005/2007 fitting data and Roser-Page +2014 validation data (crosses). + +![dose-response validation](fig1_validation_doseresponse.png) + +### Figure 2 — relative bone volume over remodeling cycles +The sawtooth trajectory (resorption dip then formation rebound each cycle) for +Wnt-10b fold changes -1, 5 and 50, with the literature endpoints overlaid. + +![bone volume vs time](fig2_bone_volume_vs_time.png) + +### Figure 3 — cell-population dynamics (Wnt-10b fold change 50) +Osteocyte, pre-osteoblast, osteoblast and osteoclast counts across the 12 cycles. + +![cell populations](fig3_cell_populations_Wnt50.png) + +The following three figures use the **numbering of the original publication**. + +### Figure 4 — model validation against Roser-Page 2014 data +Relative bone volume over 12 remodeling cycles at a Wnt-10b fold change of 1.8, +compared with the two independent Roser-Page BV/TV endpoints (held out from the +fit). Both data points lie on the simulated trajectory. (The publication's +1.2–2.4 fold-change envelope is omitted: those fold changes are not experiments +in the PEtab problem.) + +![validation vs Roser-Page](fig4_validation_roserpage.png) + +### Figure 5 — activated cell-population dynamics over a single remodeling cycle +Osteocyte, pre-osteoblast, osteoblast and osteoclast time courses for Wnt-10b +fold changes -1, 5 and 50: (A) osteocytes change little with Wnt-10b, (B) +pre-osteoblasts increase slightly, (C) osteoblasts increase, and (D) osteoclasts +decrease with increasing Wnt-10b. Populations settle well before the 100-day +cycle boundary. + +![cell dynamics over a single cycle](fig5_cell_dynamics_single_cycle.png) + +### Figure 6 — cell-population AUC ratios vs Wnt-10b +Pre-osteoblast:osteoblast (A) and osteoclast:osteoblast (B) area-under-curve +ratios over a single remodeling cycle, at the four fold changes present in the +problem (-1, 1.8, 5, 50); both ratios decrease with increasing Wnt-10b. (The +publication shows a continuous fold-change sweep, which is not available from +the problem's discrete experiments.) + +![AUC ratios](fig6_auc_ratios.png) diff --git a/Benchmark-Models/Cook_AIChE2022/conditions_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/conditions_Cook_AIChE2022.tsv new file mode 100644 index 0000000..ede4efc --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/conditions_Cook_AIChE2022.tsv @@ -0,0 +1,13 @@ +conditionId targetId targetValue +Wnt_m1 Wnt -1 +Wnt_m1 alpha3adj beta1adj + k2adj +Wnt_5 Wnt 5 +Wnt_5 alpha3adj beta1adj + k2adj +Wnt_50 Wnt 50 +Wnt_50 alpha3adj beta1adj + k2adj +Wnt_1_8 Wnt 1.8 +Wnt_1_8 alpha3adj beta1adj + k2adj +cycle_reset Osteocytes__S Osteocytes__S - 20 +cycle_reset Pre_Osteoblasts__P piecewise(0, Pre_Osteoblasts__P < 1, Pre_Osteoblasts__P) +cycle_reset Osteoblasts__B piecewise(0, Osteoblasts__B < 1, Osteoblasts__B) +cycle_reset Osteoclasts__C piecewise(0, Osteoclasts__C < 1, Osteoclasts__C) diff --git a/Benchmark-Models/Cook_AIChE2022/experiments_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/experiments_Cook_AIChE2022.tsv new file mode 100644 index 0000000..2cecc65 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/experiments_Cook_AIChE2022.tsv @@ -0,0 +1,37 @@ +experimentId conditionId time +exp_Wnt_m1 Wnt_m1 0 +exp_Wnt_m1 cycle_reset 100 +exp_Wnt_m1 cycle_reset 200 +exp_Wnt_m1 cycle_reset 300 +exp_Wnt_m1 cycle_reset 400 +exp_Wnt_m1 cycle_reset 500 +exp_Wnt_5 Wnt_5 0 +exp_Wnt_5 cycle_reset 100 +exp_Wnt_5 cycle_reset 200 +exp_Wnt_5 cycle_reset 300 +exp_Wnt_5 cycle_reset 400 +exp_Wnt_5 cycle_reset 500 +exp_Wnt_50 Wnt_50 0 +exp_Wnt_50 cycle_reset 100 +exp_Wnt_50 cycle_reset 200 +exp_Wnt_50 cycle_reset 300 +exp_Wnt_50 cycle_reset 400 +exp_Wnt_50 cycle_reset 500 +exp_Wnt_50 cycle_reset 600 +exp_Wnt_50 cycle_reset 700 +exp_Wnt_50 cycle_reset 800 +exp_Wnt_50 cycle_reset 900 +exp_Wnt_50 cycle_reset 1000 +exp_Wnt_50 cycle_reset 1100 +exp_Wnt_1_8 Wnt_1_8 0 +exp_Wnt_1_8 cycle_reset 100 +exp_Wnt_1_8 cycle_reset 200 +exp_Wnt_1_8 cycle_reset 300 +exp_Wnt_1_8 cycle_reset 400 +exp_Wnt_1_8 cycle_reset 500 +exp_Wnt_1_8 cycle_reset 600 +exp_Wnt_1_8 cycle_reset 700 +exp_Wnt_1_8 cycle_reset 800 +exp_Wnt_1_8 cycle_reset 900 +exp_Wnt_1_8 cycle_reset 1000 +exp_Wnt_1_8 cycle_reset 1100 diff --git a/Benchmark-Models/Cook_AIChE2022/fig1_validation_doseresponse.png b/Benchmark-Models/Cook_AIChE2022/fig1_validation_doseresponse.png new file mode 100644 index 0000000..9b11907 Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig1_validation_doseresponse.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/fig2_bone_volume_vs_time.png b/Benchmark-Models/Cook_AIChE2022/fig2_bone_volume_vs_time.png new file mode 100644 index 0000000..dbe635f Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig2_bone_volume_vs_time.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/fig3_cell_populations_Wnt50.png b/Benchmark-Models/Cook_AIChE2022/fig3_cell_populations_Wnt50.png new file mode 100644 index 0000000..d96d6ed Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig3_cell_populations_Wnt50.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/fig4_validation_roserpage.png b/Benchmark-Models/Cook_AIChE2022/fig4_validation_roserpage.png new file mode 100644 index 0000000..e0cce80 Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig4_validation_roserpage.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/fig5_cell_dynamics_single_cycle.png b/Benchmark-Models/Cook_AIChE2022/fig5_cell_dynamics_single_cycle.png new file mode 100644 index 0000000..e72f3d3 Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig5_cell_dynamics_single_cycle.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/fig6_auc_ratios.png b/Benchmark-Models/Cook_AIChE2022/fig6_auc_ratios.png new file mode 100644 index 0000000..6835d35 Binary files /dev/null and b/Benchmark-Models/Cook_AIChE2022/fig6_auc_ratios.png differ diff --git a/Benchmark-Models/Cook_AIChE2022/make_figures.py b/Benchmark-Models/Cook_AIChE2022/make_figures.py new file mode 100644 index 0000000..2d617b7 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/make_figures.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python3 +"""Reproduce figures from Cook et al. (2022) from the PEtab v2 problem. + +Every trajectory here is obtained by simulating the PEtab v2 problem itself +with AMICI, reusing the collection's simulation code +(``src/python/simulate.py``:``create_v2_simulator``, which encodes the +experiment periods as events). Each experiment is simulated with dense output +timepoints so that AMICI applies the ``cycle_reset`` condition at every 100-day +boundary. All figures are derived from the four dose experiments +(``exp_Wnt_m1``, ``exp_Wnt_5``, ``exp_Wnt_50``, ``exp_Wnt_1_8``). + +Because only these four Wnt-10b fold changes exist in the problem, some figures +are necessarily reduced relative to the publication: + +* Figure 4 shows the fold-change-1.8 trajectory and the Roser-Page data, but + not the 1.2-2.4 envelope (those fold changes are not part of the problem). +* Figure 6 shows the AUC ratios at the four available fold changes as discrete + points rather than a continuous sweep. + +The cycle-boundary state change makes the system stiff, so the AMICI absolute +tolerance is loosened from its very tight default (1e-16) to 1e-12; these +solver settings are passed explicitly to the simulator (see ``COOK_SOLVER``). + +Requirements: petab (v2), amici, matplotlib, numpy>=2 (a C++ compiler is needed +the first time, to build the AMICI model). +Run from anywhere: ``python make_figures.py`` (figures are written next to it). +""" +import logging +import sys +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import amici.sim.sundials as ass +import petab.v2 as petab + +HERE = Path(__file__).resolve().parent +PID = "Cook_AIChE2022" +CYCLE = 100.0 +# the cycle-reset state change makes the system stiff; loosen AMICI's very +# tight default absolute tolerance (see the module docstring). +COOK_SOLVER = dict(rtol=1e-10, atol=1e-12, max_steps=10**6) + +# reuse the collection's AMICI PEtab simulation code +sys.path.insert(0, str(HERE.parents[1] / "src" / "python")) +from simulate import create_v2_simulator # noqa: E402 + +# experiment id -> (Wnt-10b fold change, number of remodeling cycles) +EXP = {"exp_Wnt_m1": (-1.0, 6), "exp_Wnt_5": (5.0, 6), + "exp_Wnt_50": (50.0, 12), "exp_Wnt_1_8": (1.8, 12)} + +problem = petab.Problem.from_yaml(str(HERE / f"{PID}.yaml")) +nominal = problem.get_x_nominal_dict() +meas_df = problem.measurement_df # Bennett fitting data (in the problem YAML) +# Roser-Page validation data, held out of the fit (separate table, not in YAML) +val_df = pd.read_csv(HERE / f"measurementData_validation_{PID}.tsv", sep="\t") + +# --- import & compile the PEtab problem with AMICI ------------------------- +_sim = create_v2_simulator(problem, verbose=logging.WARNING, + solver_settings=COOK_SOLVER) +_model, _em, _solver = _sim.model, _sim.exp_man, _sim.solver +# states are needed for the trajectory figures +_solver.set_return_data_reporting_mode(ass.RDataReporting.full) +STATES = list(_model.get_state_ids()) # S, P, B, C, z +Z = STATES.index("Bone_volume__z") + + +def trajectory(exp_id, n_per_cycle=401): + """Dense state trajectory of a PEtab experiment simulated with AMICI. + + Returns ``(t, Y)`` with the state columns ordered as ``STATES``. AMICI + applies the ``cycle_reset`` condition at each 100-day period boundary.""" + _, ncyc = EXP[exp_id] + edata = _em.create_edata(exp_id, problem_parameters=nominal) + edata.set_timepoints(np.linspace(0, ncyc * CYCLE, ncyc * (n_per_cycle - 1) + 1)) + rdata = ass.run_simulation(_model, _solver, edata) + assert rdata.status == 0, f"{exp_id} failed with status {rdata.status}" + return np.asarray(rdata.ts), np.asarray(rdata.x) + + +TRAJ = {exp_id: trajectory(exp_id) for exp_id in EXP} + + +def first_cycle(exp_id): + """The first remodeling cycle (t < 100, before the first reset).""" + t, y = TRAJ[exp_id] + m = t < CYCLE + return t[m], y[m] + + +def save(name): + plt.tight_layout() + plt.savefig(HERE / name, dpi=120, bbox_inches="tight") + plt.close("all") + print("wrote", name) + + +# --- Figure 1: dose-response, simulation vs data ---------------------------- +# Simulated observables at the measurement timepoints, straight from the PEtab +# simulator. +_res = _sim.simulate() +print("total llh:", _res.llh) +sim_obs = {} +for rd in _res.rdatas: + if rd.y is None or np.asarray(rd.y).size == 0: + continue # e.g. the validation experiment, which has no measurements + for t, y in zip(np.atleast_1d(np.asarray(rd.ts)), + np.atleast_1d(np.asarray(rd.y).ravel())): + sim_obs[(rd.id, round(float(t)))] = float(y) + + +def _sim_obs(exp_id, time): + """Simulated observable (z - 100) for an experiment at a given time. + + Fitting experiments come straight from the PEtab simulator; the held-out + validation experiment is not in the fitted problem, so its value is taken + from the dense trajectory (z is continuous across the cycle resets).""" + key = (exp_id, round(float(time))) + if key in sim_obs: + return sim_obs[key] + t, y = TRAJ[exp_id] + return float(np.interp(float(time), t, y[:, Z])) - 100.0 + + +fig, ax = plt.subplots(figsize=(7, 6)) +colors = plt.cm.tab10.colors +i = 0 +for _, row in pd.concat([meas_df, val_df], ignore_index=True).iterrows(): + wnt = EXP[row["experimentId"]][0] + ax.plot(wnt, row["measurement"], "x", color=colors[i % 10], ms=9, mew=2) + ax.plot(wnt, _sim_obs(row["experimentId"], row["time"]), "o", + color=colors[i % 10], ms=9, label=str(row["datasetId"]).replace("_", " ")) + i += 1 +ax.set_xlabel("Wnt-10b (fold change)") +ax.set_ylabel("BV/TV (% change from normal)") +ax.set_title("Cook et al. 2022 - BV/TV change vs Wnt-10b (x: data, o: simulation)") +ax.legend(loc="best", fontsize=8) +save("fig1_validation_doseresponse.png") + +# --- Figure 2: relative bone volume over remodeling cycles ------------------ +data_points = { # relative bone volume = 100 + BV/TV % change + "exp_Wnt_m1": [(400, 100 - 29.7), (600, 100 - 41.9)], + "exp_Wnt_5": [(600, 100 + 69.2)], + "exp_Wnt_50": [(1200, 100 + 339)], +} +fig, ax = plt.subplots(figsize=(8, 6)) +for i, exp_id in enumerate(["exp_Wnt_m1", "exp_Wnt_5", "exp_Wnt_50"]): + t, y = TRAJ[exp_id] + ax.plot(t, y[:, Z], lw=1.6, color=f"C{i}", + label=f"Wnt-10b fold change {EXP[exp_id][0]:g}") + for tt, zz in data_points[exp_id]: + ax.plot(tt, zz, "x", color=f"C{i}", ms=9, mew=2) +ax.set_xlabel("time (days)") +ax.set_ylabel("relative bone volume (%)") +ax.set_title("Relative bone volume over remodeling cycles") +ax.legend(loc="best") +save("fig2_bone_volume_vs_time.png") + +# --- Figure 3: cell-population dynamics for Wnt-10b fold change 50 ---------- +t, y = TRAJ["exp_Wnt_50"] +labels = ["osteocytes (S)", "pre-osteoblasts (P)", "osteoblasts (B)", + "osteoclasts (C)"] +fig, axes = plt.subplots(2, 2, figsize=(10, 7)) +for ax, idx, lbl in zip(axes.flat, range(4), labels): + ax.plot(t, y[:, idx], color="C0", lw=1.2) + ax.set_xlabel("time (days)") + ax.set_ylabel(lbl + " (cells)") + ax.set_title(lbl) +fig.suptitle("Cell-population dynamics (Wnt-10b fold change 50)") +save("fig3_cell_populations_Wnt50.png") + +# --- Figure 4 (paper): model validation against Roser-Page 2014 data -------- +# The 1.2-2.4 envelope of the publication is omitted: those fold changes are +# not experiments in the PEtab problem. +t, y = TRAJ["exp_Wnt_1_8"] +fig, ax = plt.subplots(figsize=(8, 6)) +ax.plot(t, y[:, Z], color="C0", lw=1.8, + label="simulation (Wnt-10b fold change 1.8)") +# validation points and their reported SDs come from the validation table; +# measurements are BV/TV % change, plotted as relative bone volume (100 + %). +val_styles = {600: ("o", "C1"), 1200: ("s", "C3")} +for _, row in val_df.iterrows(): + tt = float(row["time"]) + fmt, col = val_styles.get(round(tt), ("o", "C1")) + ax.errorbar(tt, 100 + row["measurement"], row["noiseParameters"], + fmt=fmt, color=col, capsize=4, lw=2, + label=f"Roser-Page 2014 ({round(tt / CYCLE)} remodeling cycles)") +ax.set_xlabel("time (days)") +ax.set_ylabel("relative bone volume (%)") +ax.set_xlim(-10, 1210) +ax.set_title("Cook et al. 2022 Fig. 4 - validation vs Roser-Page data") +ax.legend(loc="best") +save("fig4_validation_roserpage.png") + +# --- Figure 5 (paper): activated cell-population dynamics, single cycle ------ +panels = [(0, "osteocytes (S)", "A"), (1, "pre-osteoblasts (P)", "B"), + (2, "osteoblasts (B)", "C"), (3, "osteoclasts (C)", "D")] +fig, axes = plt.subplots(2, 2, figsize=(10, 7)) +for ax, (idx, lbl, ttl) in zip(axes.flat, panels): + for exp_id in ["exp_Wnt_m1", "exp_Wnt_5", "exp_Wnt_50"]: + tc, yc = first_cycle(exp_id) + ax.plot(tc, yc[:, idx], lw=2, label=f"Wnt-10b fold change {EXP[exp_id][0]:g}") + ax.set_xlabel("time (days)") + ax.set_ylabel(lbl + " (cells)") + ax.set_title(ttl) +axes.flat[0].legend(loc="best", fontsize=8) +fig.suptitle("Cook et al. 2022 Fig. 5 - cell-population dynamics over a single cycle") +save("fig5_cell_dynamics_single_cycle.png") + +# --- Figure 6 (paper): AUC ratios vs Wnt-10b (at the available fold changes) - +# A continuous sweep is not available from the problem; the ratios are shown at +# the four experiment fold changes. +rows = [] +for exp_id, (wnt, _) in EXP.items(): + tc, yc = first_cycle(exp_id) + a_po, a_ob, a_oc = (np.trapezoid(yc[:, k], tc) for k in (1, 2, 3)) + rows.append((wnt, a_po / a_ob, a_oc / a_ob)) +rows.sort() +wnts = [r[0] for r in rows] +fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) +axes[0].plot(wnts, [r[1] for r in rows], "o-", color="C0", lw=1.5) +axes[0].set_ylabel("pre-osteoblast : osteoblast AUC") +axes[0].set_title("A") +axes[1].plot(wnts, [r[2] for r in rows], "o-", color="C3", lw=1.5) +axes[1].set_ylabel("osteoclast : osteoblast AUC") +axes[1].set_title("B") +for ax in axes: + ax.set_xlabel("Wnt-10b (fold change)") +fig.suptitle("Cook et al. 2022 Fig. 6 - cell-population AUC ratios vs Wnt-10b") +save("fig6_auc_ratios.png") diff --git a/Benchmark-Models/Cook_AIChE2022/measurementData_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/measurementData_Cook_AIChE2022.tsv new file mode 100644 index 0000000..6d606f4 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/measurementData_Cook_AIChE2022.tsv @@ -0,0 +1,5 @@ +experimentId observableId measurement time datasetId +exp_Wnt_m1 obs_BV -29.7 400 Bennett2005_4cycles +exp_Wnt_m1 obs_BV -41.9 600 Bennett2007_6cycles_Wnt_m1 +exp_Wnt_5 obs_BV 69.2 600 Bennett2007_6cycles_Wnt_5 +exp_Wnt_50 obs_BV 339.0 1200 Bennett2005_12cycles diff --git a/Benchmark-Models/Cook_AIChE2022/measurementData_validation_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/measurementData_validation_Cook_AIChE2022.tsv new file mode 100644 index 0000000..e66f53a --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/measurementData_validation_Cook_AIChE2022.tsv @@ -0,0 +1,3 @@ +experimentId observableId measurement time noiseParameters datasetId +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 diff --git a/Benchmark-Models/Cook_AIChE2022/model_Cook_AIChE2022.xml b/Benchmark-Models/Cook_AIChE2022/model_Cook_AIChE2022.xml new file mode 100644 index 0000000..028b8de --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/model_Cook_AIChE2022.xml @@ -0,0 +1,399 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + K_S + rho + + + + + + + + + + Wnt + + + Wnt + K + + + + + + + + + + + + + + + + + + + + Bone + alpha_1 + + + Osteoblasts__B + g_31 + + + + 1 + + + Osteocytes__S + K_S + + + + + + + + + + + + + + + + + + + Bone + alpha_2 + + + Osteocytes__S + g_21 + + + + 1 + + + Osteocytes__S + K_S + + + g_22 + + + + + + + + + + + + + + + + + + + Bone + alpha_3 + + + Pre_Osteoblasts__P + g_32 + + + + 1 + + + Osteocytes__S + K_S + + + + + + alpha3adj + piwnta + + + Pre_Osteoblasts__P + f_12 + + + + + + + + + + + + + + + + + + + + + + Bone + beta_1 + + + Pre_Osteoblasts__P + f_12 + + + + Osteoclasts__C + f_14 + + + + + beta1adj + piwnta + + + Pre_Osteoblasts__P + f_12 + + + + + + + + + + + + + + + Bone + delta + Pre_Osteoblasts__P + + + + + + + + + + + + + Bone + + + beta_2 + + + beta2adj + piwnta + + + + + Osteoblasts__B + f_23 + + + + + + + + + + + + + + + + + + + + + + Bone + alpha_4 + + + Osteocytes__S + g_41 + + + + Pre_Osteoblasts__P + g_42 + + + + + + epsilon + Osteoblasts__B + + g_43 + + + + + + 1 + + + Osteocytes__S + K_S + + + g_44 + + + + + + + + + + + + + + Bone + beta_3 + + + Osteoclasts__C + f_34 + + + + + + + + + + + + + + + + + + Bone + k1 + Osteoclasts__C + + + + + + + + + + + + + + + + + Bone + k2 + Osteoblasts__B + + + + + + + diff --git a/Benchmark-Models/Cook_AIChE2022/observables_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/observables_Cook_AIChE2022.tsv new file mode 100644 index 0000000..143df31 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/observables_Cook_AIChE2022.tsv @@ -0,0 +1,2 @@ +observableId observableName observableFormula noiseFormula +obs_BV relative bone volume change Bone_volume__z - 100 1 diff --git a/Benchmark-Models/Cook_AIChE2022/parameters_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/parameters_Cook_AIChE2022.tsv new file mode 100644 index 0000000..21bb24a --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/parameters_Cook_AIChE2022.tsv @@ -0,0 +1,5 @@ +parameterId parameterName lowerBound upperBound nominalValue estimate +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 diff --git a/Benchmark-Models/Cook_AIChE2022/simulatedData_Cook_AIChE2022.tsv b/Benchmark-Models/Cook_AIChE2022/simulatedData_Cook_AIChE2022.tsv new file mode 100644 index 0000000..911c7a8 --- /dev/null +++ b/Benchmark-Models/Cook_AIChE2022/simulatedData_Cook_AIChE2022.tsv @@ -0,0 +1,5 @@ +experimentId observableId simulation time datasetId +exp_Wnt_m1 obs_BV -22.0103479953 400 Bennett2005_4cycles +exp_Wnt_m1 obs_BV -33.0155221687 600 Bennett2007_6cycles_Wnt_m1 +exp_Wnt_5 obs_BV 80.9940457566 600 Bennett2007_6cycles_Wnt_5 +exp_Wnt_50 obs_BV 337.5034729186 1200 Bennett2005_12cycles diff --git a/README.md b/README.md index 307a57a..1f60733 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,7 @@ A sortable version of the table is available | [Brannmark_JBC2010](Benchmark-Models/Brannmark_JBC2010/) | 8 | 22 | 0 | ✓ | 1 | 0 | 43 | 3 | normal | | 9 | [\[1\]](http://identifiers.org/doi/10.1074/jbc.M110.106849) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Brannmark_JBC2010/model_Brannmark_JBC2010.xml) | | [Bruno_JExpBot2016](Benchmark-Models/Bruno_JExpBot2016/) | 6 | 13 | 0 | | 0 | 0 | 77 | 5 | normal | | 7 | [\[1\]](http://identifiers.org/doi/10.1093/jxb/erw356) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Bruno_JExpBot2016/model_Bruno_JExpBot2016.xml) | | [Chen_MSB2009](Benchmark-Models/Chen_MSB2009/) | 4 | 155 | 0 | ✓ | 0 | 0 | 120 | 3 | normal | | 500 | [\[1\]](http://identifiers.org/doi/10.1038/msb.2008.74) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Chen_MSB2009/model_Chen_MSB2009.xml) | +| [Cook_AIChE2022](Benchmark-Models/Cook_AIChE2022/) | 4 | 4 | 0 | | 0 | 0 | 6 | 1 | normal | | 5 | [\[1\]](http://identifiers.org/doi/10.1002/aic.17809) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Cook_AIChE2022/model_Cook_AIChE2022.xml) | | [Crauste_CellSystems2017](Benchmark-Models/Crauste_CellSystems2017/) | 1 | 12 | 0 | | 0 | 0 | 21 | 4 | normal | | 5 | [\[1\]](http://identifiers.org/doi/10.1016/j.cels.2017.01.014) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Crauste_CellSystems2017/model_Crauste_CellSystems2017.xml) | | [Elowitz_Nature2000](Benchmark-Models/Elowitz_Nature2000/) | 1 | 21 | 0 | | 0 | 0 | 58 | 1 | log10-normal | | 8 | [\[1\]](http://identifiers.org/doi/10.1038/35002125) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Elowitz_Nature2000/model_Elowitz_Nature2000.xml) | | [Fiedler_BMCSystBiol2016](Benchmark-Models/Fiedler_BMCSystBiol2016/) | 3 | 22 | 0 | | 0 | 0 | 72 | 2 | normal | | 6 | [\[1\]](http://identifiers.org/doi/10.1186/s12918-016-0319-7) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Fiedler_BMCSystBiol2016/model_Fiedler_BMCSystBiol2016.xml) | diff --git a/src/python/benchmark_models_petab/base.py b/src/python/benchmark_models_petab/base.py index 8a81329..806ee88 100644 --- a/src/python/benchmark_models_petab/base.py +++ b/src/python/benchmark_models_petab/base.py @@ -2,12 +2,16 @@ from pathlib import Path -import petab.v1 as petab +import petab.v1 +import petab.v2 +from petab.versions import get_major_version from .C import MODELS_DIR import pandas as pd +Problem = petab.v1.Problem | petab.v2.Problem + def get_problem_yaml_path(id_: str) -> Path: """Get the path to the PEtab problem YAML file. @@ -28,7 +32,7 @@ def get_problem_yaml_path(id_: str) -> Path: return yaml_path -def get_problem(id_: str) -> petab.Problem: +def get_problem(id_: str) -> Problem: """Read PEtab problem from the benchmark collection by name. Parameters @@ -37,11 +41,13 @@ def get_problem(id_: str) -> petab.Problem: Returns ------- - The PEtab problem. + The PEtab problem (a ``petab.v1`` or ``petab.v2`` problem, depending on + the problem's ``format_version``). """ yaml_file = get_problem_yaml_path(id_) - petab_problem = petab.Problem.from_yaml(yaml_file) - return petab_problem + if get_major_version(yaml_file) >= 2: + return petab.v2.Problem.from_yaml(yaml_file) + return petab.v1.Problem.from_yaml(yaml_file) def get_simulation_df(id_: str) -> pd.DataFrame | None: @@ -58,6 +64,6 @@ def get_simulation_df(id_: str) -> pd.DataFrame | None: """ for filename in (f"simulatedData_{id_}.tsv", "simulations.tsv"): if (path := Path(MODELS_DIR, id_, filename)).is_file(): - return petab.get_simulation_df(path) + return petab.v1.get_simulation_df(path) return None diff --git a/src/python/benchmark_models_petab/check_sbml_metadata.py b/src/python/benchmark_models_petab/check_sbml_metadata.py index 7e46cf0..c667f3b 100644 --- a/src/python/benchmark_models_petab/check_sbml_metadata.py +++ b/src/python/benchmark_models_petab/check_sbml_metadata.py @@ -4,7 +4,7 @@ from pathlib import Path from . import MODELS, get_problem -from .overview import get_reference_uris +from .overview import get_reference_uris, get_sbml_model def main(): @@ -14,11 +14,10 @@ def main(): for petab_problem_id in MODELS: petab_problem = get_problem(petab_problem_id) - model_id = petab_problem.sbml_model.getId() - model_name = petab_problem.sbml_model.getName() - reference_uris = get_reference_uris( - sbml_model=petab_problem.sbml_model - ) + sbml_model = get_sbml_model(petab_problem) + model_id = sbml_model.getId() + model_name = sbml_model.getName() + reference_uris = get_reference_uris(sbml_model=sbml_model) errors = [] diff --git a/src/python/benchmark_models_petab/overview.py b/src/python/benchmark_models_petab/overview.py index 52f5e55..0d4e38b 100644 --- a/src/python/benchmark_models_petab/overview.py +++ b/src/python/benchmark_models_petab/overview.py @@ -39,21 +39,32 @@ index_column = "petab_problem_id" +def get_sbml_model(petab_problem) -> libsbml.Model: + """Get the libSBML model of a v1 or v2 PEtab problem.""" + return getattr(petab_problem, "sbml_model", None) or ( + petab_problem.model.sbml_model + ) + + +def _num_conditions(petab_problem) -> int: + """Number of simulation conditions (v1) or experiments (v2).""" + if hasattr(petab_problem, "get_simulation_conditions_from_measurement_df"): + df = petab_problem.get_simulation_conditions_from_measurement_df() + return df.shape[0] + return len(petab_problem.experiments) + + def get_summary( - petab_problem: petab.Problem, + petab_problem, petab_problem_id: str = None, ) -> Dict: """Get dictionary with stats for the given PEtab problem""" print(petab_problem_id) return { "petab_problem_id": petab_problem_id, - "conditions": petab_problem.get_simulation_conditions_from_measurement_df().shape[ - 0 - ], - "estimated_parameters": np.sum( - petab_problem.parameter_df[petab.ESTIMATE] - ), - "events": len(petab_problem.sbml_model.getListOfEvents()), + "conditions": _num_conditions(petab_problem), + "estimated_parameters": len(petab_problem.x_free_ids), + "events": len(get_sbml_model(petab_problem).getListOfEvents()), "possible_discontinuities": guess_discontinuities(petab_problem), "preequilibration": 0 if petab.PREEQUILIBRATION_CONDITION_ID @@ -85,8 +96,8 @@ def get_summary( "objective_prior_distributions": get_prior_distributions( petab_problem.parameter_df ), - "species": len(petab_problem.sbml_model.getListOfSpecies()), - "reference_uris": get_reference_uris(petab_problem.sbml_model), + "species": len(get_sbml_model(petab_problem).getListOfSpecies()), + "reference_uris": get_reference_uris(get_sbml_model(petab_problem)), "sbml4humans_urls": get_sbml4humans_urls(petab_problem_id), } @@ -110,11 +121,21 @@ def get_sbml4humans_urls(petab_problem_id: str) -> List[str]: yaml_file = get_problem_yaml_path(petab_problem_id) yaml_dict = load_yaml(yaml_file) repo_root = "https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master" + # collect model files across PEtab v1 (problems/sbml_files) and + # v2 (model_files//location) layouts + model_filenames = [] + if petab.PROBLEMS in yaml_dict: + for problem_dict in yaml_dict[petab.PROBLEMS]: + model_filenames.extend(problem_dict.get(petab.SBML_FILES, [])) + else: + for model_dict in yaml_dict.get("model_files", {}).values(): + model_filenames.append(model_dict["location"]) urls = [] - for problem_dict in yaml_dict[petab.PROBLEMS]: - for model_filename in problem_dict[petab.SBML_FILES]: - gh_raw_url = f"{repo_root}/Benchmark-Models/{petab_problem_id}/{model_filename}" - urls.append(f"https://sbml4humans.de/model_url?url={gh_raw_url}") + for model_filename in model_filenames: + gh_raw_url = ( + f"{repo_root}/Benchmark-Models/{petab_problem_id}/{model_filename}" + ) + urls.append(f"https://sbml4humans.de/model_url?url={gh_raw_url}") return urls @@ -187,7 +208,7 @@ def guess_discontinuities(petab_problem: petab.Problem) -> bool: discontinuities might play a role in the simulation; the presence of discontinuities might be parameter-dependent; ... . """ - model: libsbml.Model = petab_problem.sbml_model + model: libsbml.Model = get_sbml_model(petab_problem) for event in model.getListOfEvents(): for ea in event.getListOfEventAssignments(): @@ -195,7 +216,7 @@ def guess_discontinuities(petab_problem: petab.Problem) -> bool: return True # convert reactions to rate rules - sbml_doc = petab_problem.sbml_model.getSBMLDocument().clone() + sbml_doc = get_sbml_model(petab_problem).getSBMLDocument().clone() model = sbml_doc.getModel() conversion_config = libsbml.ConversionProperties() conversion_config.addOption("replaceReactions") diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml index bb87e6a..98933b7 100644 --- a/src/python/pyproject.toml +++ b/src/python/pyproject.toml @@ -24,11 +24,13 @@ dependencies = [ "python-libsbml", "tabulate", "sbmlmath", + "scipy", ] [project.optional-dependencies] dev = ["pre-commit", "pytest", "ruff"] site = ["bokeh>=3.7.3"] +simulate = ["amici[petab]"] [project.scripts] bmp-petablint = "benchmark_models_petab.check_petablint:main" diff --git a/src/python/simulate.py b/src/python/simulate.py index e867a20..442288b 100755 --- a/src/python/simulate.py +++ b/src/python/simulate.py @@ -2,23 +2,113 @@ """ Script to simulate a model using AMICI. +Handles both PEtab v1 and PEtab v2 problems (see :func:`simulate`). + Basic usage: simulate.py See --help for more options. """ -import os import argparse -import benchmark_models_petab -import amici.sim.sundials as ass -from amici.importers.petab.v1 import ( - import_petab_problem, - simulate_petab, - LLH, - RDATAS, -) import logging -from petab.v1.lint import measurement_table_has_timepoint_specific_mappings +import os + +import amici.sim.sundials as ass +import benchmark_models_petab +import petab.v2 +from amici.importers.petab import PetabImporter +from amici.importers.petab.v1 import import_petab_problem +from amici.sim.sundials.petab.v1 import LLH, RDATAS, simulate_petab from petab.v1.core import flatten_timepoint_specific_output_overrides +from petab.v1.lint import measurement_table_has_timepoint_specific_mappings + + +# Per-problem AMICI solver overrides for PEtab v2 problems that need looser +# tolerances than AMICI's very tight defaults. Problems not listed here keep +# the AMICI defaults, so this is not applied globally to all v2 problems. +_V2_SOLVER_OVERRIDES = { + # The cycle-reset state change makes this system stiff; the default + # atol=1e-16 triggers a too-small step size right after a reset. + "Cook_AIChE2022": dict(rtol=1e-10, atol=1e-12, max_steps=10**6), +} + + +def _apply_solver_settings(solver, solver_settings: dict | None) -> None: + """Apply optional AMICI solver tolerance overrides. + + Leaves AMICI's defaults untouched when ``solver_settings`` is ``None``. + Recognized keys: ``rtol``, ``atol``, ``max_steps``. + """ + if not solver_settings: + return + if "rtol" in solver_settings: + solver.set_relative_tolerance(solver_settings["rtol"]) + if "atol" in solver_settings: + solver.set_absolute_tolerance(solver_settings["atol"]) + if "max_steps" in solver_settings: + solver.set_max_steps(solver_settings["max_steps"]) + + +def create_v2_simulator( + problem: "petab.v2.Problem", + *, + num_threads: int = 1, + verbose: int = logging.INFO, + solver_settings: dict | None = None, +): + """Create an AMICI PEtab simulator for a PEtab v2 problem. + + ``solver_settings`` optionally overrides the AMICI solver tolerances + (keys ``rtol``, ``atol``, ``max_steps``); AMICI defaults are used when it + is ``None``. + """ + importer = PetabImporter(problem, verbose=verbose) + simulator = importer.create_simulator() + simulator.num_threads = num_threads + _apply_solver_settings(simulator.solver, solver_settings) + return simulator + + +def _simulate_v1(problem, num_threads): + """Simulate a PEtab v1 problem; returns (total_llh, rdatas).""" + if measurement_table_has_timepoint_specific_mappings( + problem.measurement_df, + ): + flatten_timepoint_specific_output_overrides(problem) + + model = import_petab_problem( + problem, generate_sensitivity_code=False, verbose=logging.INFO + ) + solver = model.create_solver() + solver.set_return_data_reporting_mode( + ass.RDataReporting.observables_likelihood + ) + res = simulate_petab( + problem, amici_model=model, solver=solver, num_threads=num_threads + ) + return res[LLH], res[RDATAS] + + +def _simulate_v2(problem, num_threads, solver_settings=None): + """Simulate a PEtab v2 problem; returns (total_llh, rdatas).""" + simulator = create_v2_simulator( + problem, num_threads=num_threads, solver_settings=solver_settings + ) + res = simulator.simulate() + return res.llh, res.rdatas + + +def simulate( + problem, num_threads: int = 1, solver_settings: dict | None = None +): + """Simulate a PEtab problem at its nominal parameters with AMICI. + + Dispatches to the PEtab v1 or v2 simulation route depending on the type of + ``problem``. ``solver_settings`` (v2 only) optionally overrides the AMICI + solver tolerances. Returns ``(total_log_likelihood, rdatas)``. + """ + if isinstance(problem, petab.v2.Problem): + return _simulate_v2(problem, num_threads, solver_settings) + return _simulate_v1(problem, num_threads) def parse_cli_args(): @@ -66,29 +156,21 @@ def main(): problem = benchmark_models_petab.get_problem(args.problem_id) - if measurement_table_has_timepoint_specific_mappings( - problem.measurement_df, - ): - flatten_timepoint_specific_output_overrides(problem) - - model = import_petab_problem( - problem, generate_sensitivity_code=False, verbose=logging.INFO - ) - solver = model.create_solver() - solver.set_return_data_reporting_mode( - ass.RDataReporting.observables_likelihood - ) - res = simulate_petab( - problem, model, solver=solver, num_threads=args.num_threads + total_llh, rdatas = simulate( + problem, + num_threads=args.num_threads, + solver_settings=_V2_SOLVER_OVERRIDES.get(args.problem_id), ) - for rdata in res[RDATAS]: + for rdata in rdatas: + chi2 = getattr(rdata, "chi2", None) print( - f"{rdata.id}: llh = {rdata.llh}, chi2 = {rdata.chi2}, status = {rdata.status}" + f"{rdata.id}: llh = {rdata.llh}, chi2 = {chi2}, " + f"status = {rdata.status}" ) print() - print("Total log-likelihood:", res[LLH]) + print("Total log-likelihood:", total_llh) print()