diff --git a/Benchmark-Models/Perelson_Science1996/README.md b/Benchmark-Models/Perelson_Science1996/README.md new file mode 100644 index 00000000..1b5cc72e --- /dev/null +++ b/Benchmark-Models/Perelson_Science1996/README.md @@ -0,0 +1,149 @@ +# Implementation details of the Perelson_Science1996 benchmark problem + +## Model + +Linear three-state model of HIV-1 dynamics in vivo under a protease inhibitor +[[1]]. Ritonavir renders newly produced virions non-infectious but does not stop +production by already infected cells, so the plasma virion population splits +into an infectious pool present when the drug took effect and a non-infectious +pool produced afterwards: + +``` +d(Tstar)/dt = k*T0*V_I - delta*Tstar (Eq. 3) +d(V_I)/dt = -c*V_I (Eq. 4) +d(V_NI)/dt = N_virions*delta*Tstar - c*V_NI (Eq. 5) +V = V_I + V_NI +``` + +The concentration of uninfected target cells is held at its pre-treatment value +`T0`, and the drug is assumed to be 100 % effective. + +| id | paper | meaning | +|----|-------|---------| +| `Tstar` | *T\** | productively infected cells | +| `V_I` | *V_I* | infectious virions, present when the drug took effect | +| `V_NI` | *V_NI* | non-infectious virions, produced after the drug took effect | +| `V` | *V* | total plasma viral load; the measured quantity | +| `c` | *c* | virion clearance rate [1/day] | +| `delta` | *δ* | loss rate of virus-producing cells [1/day] | +| `V0` | *V_0* | viral load when the drug took effect [RNA copies/ml] | +| `N_virions` | *N* | virions produced per infected cell over its lifetime | +| `T0` | *T_0* | pre-treatment concentration of uninfected target cells | +| `k` | *k* | infection rate constant | + +`Tstar` is spelled without an asterisk because SBML identifiers cannot contain +one. The burst size is `N_virions` rather than `N` because `N` is the name of +SymPy's numeric-evaluation function, and AMICI parses SBML initial assignments +through SymPy — a parameter called `N` makes the model fail to import. + +## The estimation problem + +A single condition, `condition1`, no pre-equilibration, and one observable, +`patient_105`, which is `V` with a `log10` transformation and normally +distributed noise — that is, residuals are normal in `log10` space. The +measurement table holds 13 plasma HIV-1 RNA concentrations over 6.75 days. + +Three parameters are estimated, all on a `log10` scale: `c`, `delta` and `V0`. +These are exactly the three that [[1]] states were estimated simultaneously, and +their nominal values are the ones published for this patient. `N_virions = 480` +and `T0 = 11000` are fixed. + +There is no noise parameter: the observable table gives `noiseFormula` as the +literal `1`, so the objective is the unweighted sum of squared `log10` +residuals — the nonlinear least squares problem solved in [[1]]. That 1 is a +weighting convention rather than an estimate of the measurement noise; since +the observable is `log10`-transformed it would be a residual width in `log10` +units, and the root mean square `log10` residual at the optimum is 0.117. +Weighting the residuals equally does not bias the estimates: with a single +observable and one shared noise level, estimating that level instead would +profile out to a strictly monotone function of the same sum of squares, giving +the same optimum in `c`, `delta` and `V0`. + +## Data + +The measurements are the raw plasma HIV-1 RNA concentrations (copies/ml) of +**patient 105** of [[1]], taken from Hulin Wu's *Early intensive viral dynamic +data* release [[2]], where this patient carries `ID = 4`. + +Table 1 of [[1]] reports base-line plasma **virions**/ml and states that each +virion carries two RNA copies, whereas the release reports RNA copies/ml. +Halving the day-0 release values reproduces the Table 1 base lines for all five +patients, in ascending order of ID, which is how the release IDs map onto the +patient numbers used in the publication: + +| release ID | RNA/ml at day 0 | → virions/ml | patient | Table 1 base line | ratio | +|-----------:|----------------:|-------------:|--------:|------------------:|------:| +| 1 | 610,000 | 305,000 | 102 | 294,000 | 1.04 | +| 2 | 19,000 | 9,500 | 103 | 12,000 | 0.79 | +| 3 | 100,000 | 50,000 | 104 | 52,000 | 0.96 | +| 4 | 1,022,000 | 511,000 | 105 | 643,000 | 0.80 | +| 5 | 160,000 | 80,000 | 107 | 77,000 | 1.04 | + +The residual spread is expected: the Table 1 base lines are averages over days +−7, −4, −1 and 0, while the release starts at day 0. Three of the five +assignments are confirmed independently by Fig. 1 of [[1]], which plots the +released series point for point: Fig. 1A upper is patient 104, Fig. 1A lower is +patient 107, and Fig. 1B is patient 105. + +## Time base + +Time is measured **from the onset of the drug effect, not from dosing**. The +two differ by a pharmacologic delay, 6 h = 0.25 d for patient 105, measured for +this patient from the first drop in plasma infectivity (Table 1 of [[1]] and its +legend). The model has no validity before that onset, so the three measurements +taken earlier — released days 0, 0.083 and 0.166 — are not part of the problem, +and 0.25 d has been subtracted from the sampling times of the remaining 13. + +## Initial conditions + +Only `V0` is free; the other initial values follow from the pre-treatment quasi +steady state assumed in [[1]] and are encoded as SBML initial assignments, so +they track `V0`, `c` and `delta` as those are estimated: + +| state | initial value | origin | +|-------|---------------|--------| +| `V_I` | `V0` | Eq. 4, `V_I(0) = V0` | +| `V_NI` | `0` | Eq. 5, no virions produced under the drug yet | +| `V` | `V0` | `V = V_I + V_NI` | +| `Tstar` | `c*V0/(N_virions*delta)` | `dV/dt = 0` at steady state (Eqs. 1–2) | + +The same steady state fixes the infection rate constant, `c = N_virions*k*T0`, +so `k` is likewise an initial assignment, `c/(N_virions*T0)`, and not a free +parameter. This matters because `c` is estimated: holding `k` at a literal value +would let the model drift away from Eq. 6 of [[1]] as `c` moves — by up to 13 % +at `c = 1.5` — whereas with `k` coupled the model reproduces Eq. 6 at any +parameter value. + +`T0` cancels out of the dynamics under that coupling (`k*T0 = c/N_virions`). It +is kept as a fixed parameter because [[1]] reports it and because the coupling is +clearer written with it. + +## Validation + +`simulate_Perelson_Science1996.py` in this directory simulates the problem with +AMICI at the nominal parameter values and writes +`simulatedData_Perelson_Science1996.tsv`, so the simulated-data table can be +regenerated rather than taken on trust. Those simulations agree with the +closed-form solution, Eq. 6 of [[1]], to a relative 2e-8 at AMICI's default +tolerances. + +Re-estimating `c`, `delta` and `V0` from the measurement table recovers the +values published for patient 105: + +| | `c` | `delta` | `V0` | +|---|---|---|---| +| published (Table 1) | 2.06 | 0.53 | 1.86e6 | +| re-estimated here | 2.062 | 0.525 | 1.832e6 | + +The nominal parameters therefore sit essentially at the optimum: the sum of +squared `log10` residuals is 0.1775 at the nominal values and 0.1771 at the +optimum. + +## Deviations from the publication + +- [[1]] eliminated one outlying data point per patient by a jackknife procedure + before fitting. That step is not reproduced here; it is not needed to recover + the published parameter values. + +[1]: https://doi.org/10.1126/science.271.5255.1582 +[2]: https://sph.uth.edu/dept/bads/faculty-home/hulinwu/datasets/early-intensive-viral-dynamic-data diff --git a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv index 48225c17..f96cfecc 100644 --- a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv @@ -1,17 +1,14 @@ -observableId preequilibrationConditionId simulationConditionId measurement time observableParameters noiseParameters -task0_model0_perelson1_V condition1 1029000 0 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1087000 0.105 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2158000 0.169 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1860000 0.282 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2195000 0.492 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1567000 0.757 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 3208000 1.029 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2293000 1.253 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1342000 1.533 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1197000 1.75 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 987100 2.038 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 697300 3.013 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 503800 3.987 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 282200 4.968 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 172700 5.972 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 91618.242 6.973 sd_task0_model0_perelson1_V +observableId preequilibrationConditionId simulationConditionId measurement time observableParameters +patient_105 condition1 1043000 0 +patient_105 condition1 2170000 0.25 +patient_105 condition1 1538000 0.5 +patient_105 condition1 3098000 0.75 +patient_105 condition1 2262000 1 +patient_105 condition1 1315000 1.25 +patient_105 condition1 1200000 1.5 +patient_105 condition1 980000 1.75 +patient_105 condition1 710000 2.75 +patient_105 condition1 500000 3.75 +patient_105 condition1 280000 4.75 +patient_105 condition1 170000 5.75 +patient_105 condition1 89000 6.75 diff --git a/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml b/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml index 88a16256..f50551b0 100644 --- a/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml +++ b/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml @@ -37,33 +37,76 @@ - - - - + + + + - + - + + + + + + + + c + + + N_virions + T0 + + + + + + + + + + + c + V0 + + + + N_virions + delta + + + + + + + V0 + + + + + V0 + + + - + - K0 + k T0 - Vin + V_I @@ -85,14 +128,14 @@ - + c - Vin + V_I @@ -100,14 +143,14 @@ - + c - Vni + V_NI @@ -115,7 +158,7 @@ - + @@ -125,7 +168,7 @@ delta - NN + N_virions Tstar diff --git a/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv index cec307df..40d9b677 100644 --- a/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv @@ -1,2 +1,2 @@ observableId observableFormula observableName noiseDistribution noiseFormula observableTransformation -task0_model0_perelson1_V V normal noiseParameter1_task0_model0_perelson1_V log10 +patient_105 V normal 1 log10 diff --git a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv index d11162a5..c163363b 100644 --- a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv @@ -1,7 +1,6 @@ parameterId parameterName parameterScale lowerBound upperBound nominalValue estimate -NN NN log10 1E-05 100000 480.000000026982 0 +N_virions N_virions log10 1E-05 100000 480.000000026982 0 T0 T0 log10 1E-05 100000 11000.0000010581 0 c c log10 1E-05 100000 2.06000000014632 1 delta delta log10 1E-05 100000 0.529999999999037 1 -K0 K0 log10 1E-05 100000 3.89999999976204E-07 0 -sd_task0_model0_perelson1_V sigma_task0_model0_perelson1_V log10 1E-10 10000000000 100000 1 +V0 V0 log10 1E+03 1E+09 1860000 1 diff --git a/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv index e68cd7a8..7cee41f3 100644 --- a/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv @@ -1,17 +1,14 @@ -observableId preequilibrationConditionId simulationConditionId simulation time observableParameters noiseParameters observableTransformation noiseDistribution -task0_model0_perelson1_V condition1 1860000 0 1 log10 normal -task0_model0_perelson1_V condition1 1979067.19116451 0.105 1 log10 normal -task0_model0_perelson1_V condition1 2036200.82349015 0.169 1 log10 normal -task0_model0_perelson1_V condition1 2110605.53393048 0.282 1 log10 normal -task0_model0_perelson1_V condition1 2170386.58687064 0.492 1 log10 normal -task0_model0_perelson1_V condition1 2133014.95269381 0.757 1 log10 normal -task0_model0_perelson1_V condition1 2007892.04532337 1.029 1 log10 normal -task0_model0_perelson1_V condition1 1866529.16992371 1.253 1 log10 normal -task0_model0_perelson1_V condition1 1667316.74545722 1.533 1 log10 normal -task0_model0_perelson1_V condition1 1508578.7515073 1.75 1 log10 normal -task0_model0_perelson1_V condition1 1304224.95565882 2.038 1 log10 normal -task0_model0_perelson1_V condition1 745147.675932322 3.013 1 log10 normal -task0_model0_perelson1_V condition1 404793.605719818 3.987 1 log10 normal -task0_model0_perelson1_V condition1 214718.67877913 4.968 1 log10 normal -task0_model0_perelson1_V condition1 111378.497695453 5.972 1 log10 normal -task0_model0_perelson1_V condition1 57731.3816378903 6.973 1 log10 normal +observableId preequilibrationConditionId simulationConditionId simulation time observableParameters observableTransformation noiseDistribution +patient_105 condition1 1860000 0 log10 normal +patient_105 condition1 1851810.36689883 0.25 log10 normal +patient_105 condition1 1810225.3405842 0.5 log10 normal +patient_105 condition1 1731023.03438251 0.75 log10 normal +patient_105 condition1 1622815.16935044 1 log10 normal +patient_105 condition1 1496900.50271573 1.25 log10 normal +patient_105 condition1 1363269.14615056 1.5 log10 normal +patient_105 condition1 1229442.38119697 1.75 log10 normal +patient_105 condition1 767132.663106335 2.75 log10 normal +patient_105 condition1 459200.33082893 3.75 log10 normal +patient_105 condition1 271534.501982976 4.75 log10 normal +patient_105 condition1 160019.934111288 5.75 log10 normal +patient_105 condition1 94217.5663873072 6.75 log10 normal diff --git a/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv index 262254a7..98a4216e 100644 --- a/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv @@ -1,2 +1,2 @@ plotId plotTypeData plotTypeSimulation xLabel xScale yLabel yScale legendEntry yValues xValues xOffset yOffset -plot1 MeanAndSD LinePlot time [days] lin RNA copies/ml lin HIV-mRNA copies - task0_model0_perelson1_V task0_model0_perelson1_V time 0 0 +plot1 MeanAndSD LinePlot time since drug effect [days] lin RNA copies/ml lin HIV-mRNA copies - patient_105 patient_105 time 0 0 diff --git a/README.md b/README.md index 307a57a4..3a9af9f0 100644 --- a/README.md +++ b/README.md @@ -38,7 +38,7 @@ A sortable version of the table is available | [Lucarelli_CellSystems2018](Benchmark-Models/Lucarelli_CellSystems2018/) | 16 | 84 | 0 | | 0 | 0 | 1755 | 65 | normal; log10-normal | | 33 | [\[1\]](http://identifiers.org/doi/10.1016/j.cels.2017.11.010) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Lucarelli_CellSystems2018/model_Lucarelli_CellSystems2018.xml) | | [Okuonghae_ChaosSolitonsFractals2020](Benchmark-Models/Okuonghae_ChaosSolitonsFractals2020/) | 1 | 16 | 0 | | 0 | 0 | 92 | 2 | normal | | 9 | [\[1\]](http://identifiers.org/doi/10.1016/j.chaos.2020.110032) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Okuonghae_ChaosSolitonsFractals2020/model_Okuonghae_ChaosSolitonsFractals2020.xml) | | [Oliveira_NatCommun2021](Benchmark-Models/Oliveira_NatCommun2021/) | 1 | 12 | 0 | ✓ | 0 | 0 | 120 | 2 | normal | | 9 | [\[1\]](http://identifiers.org/doi/10.1038/s41467-020-19798-3) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Oliveira_NatCommun2021/Oliveira_NatCommun2021_model.xml) | -| [Perelson_Science1996](Benchmark-Models/Perelson_Science1996/) | 1 | 3 | 0 | | 0 | 0 | 16 | 1 | log10-normal | | 4 | [\[1\]](http://identifiers.org/doi/10.1126/science.271.5255.1582) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml) | +| [Perelson_Science1996](Benchmark-Models/Perelson_Science1996/) | 1 | 3 | 0 | | 0 | 0 | 13 | 1 | log10-normal | | 4 | [\[1\]](http://identifiers.org/doi/10.1126/science.271.5255.1582) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml) | | [Rahman_MBS2016](Benchmark-Models/Rahman_MBS2016/) | 1 | 9 | 0 | | 0 | 0 | 23 | 1 | normal | | 7 | [\[1\]](http://identifiers.org/doi/10.1016/j.mbs.2016.07.009) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Rahman_MBS2016/model_Rahman_MBS2016.xml) | | [Raia_CancerResearch2011](Benchmark-Models/Raia_CancerResearch2011/) | 4 | 39 | 0 | | 0 | 0 | 205 | 8 | normal | | 14 | [\[1\]](http://identifiers.org/doi/10.1158/0008-5472.CAN-10-2987) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Raia_CancerResearch2011/model_Raia_CancerResearch2011.xml) | | [Raimundez_PCB2020](Benchmark-Models/Raimundez_PCB2020/) | 170 | 136 | 0 | ✓ | 4 | 0 | 627 | 79 | normal | parameterScaleNormal; laplace | 22 | [\[1\]](http://identifiers.org/doi/10.1371/journal.pcbi.1007147) | [\[1\]](https://sbml4humans.de/model_url?url=https://raw.githubusercontent.com/Benchmarking-Initiative/Benchmark-Models-PEtab/master/Benchmark-Models/Raimundez_PCB2020/model_Raimundez_PCB2020.xml) |