From cd82660517963cf28c970c80ddafb23eca175cb0 Mon Sep 17 00:00:00 2001 From: Maren Philipps Date: Sat, 19 Sep 2026 13:09:08 +0000 Subject: [PATCH 1/4] Perelson_Science1996: use Wu data release instead of digitised data (#322) The measurement table held crudely digitised values read off the figure in Perelson et al. (1996). The underlying data are publicly released by Hulin Wu's group as "Early Intensive Viral Dynamic Data"; the PEtab problem uses patient ID 4 there (median ratio of the digitised to the released series 1.010; no other patient in the release is within a factor of 8). Replace all 16 values and their sampling times with the released ones. The digitised values were accurate to a few percent except at t=0.25 d, which was off by a factor of 1.78 (1.86e6 instead of 1.04e6 copies/ml). The times were read off the x-axis rather than taken from the nominal sampling schedule (hours 0/2/4/6/12/18/24/30/36/42/48, then days 3-7), drifting by up to 46 min. Refitting the problem as-is to the corrected data moves the optimum to c = 1.907, delta = 0.529, closer to the published c = 2.06, delta = 0.53 than the digitised data's optimum (c = 1.861, delta = 0.547). Data source: https://sph.uth.edu/dept/bads/faculty-home/hulinwu/datasets/early-intensive-viral-dynamic-data Co-Authored-By: Claude Opus 5 --- .../measurementData_Perelson_Science1996.tsv | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv index 48225c17..2094a65d 100644 --- a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv @@ -1,17 +1,17 @@ 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 +task0_model0_perelson1_V condition1 1022000 0 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 1079000 0.083 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 2148000 0.166 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 1043000 0.25 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 2170000 0.5 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 1538000 0.75 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 3098000 1 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 2262000 1.25 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 1315000 1.5 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 1200000 1.75 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 980000 2 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 710000 3 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 500000 4 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 280000 5 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 170000 6 sd_task0_model0_perelson1_V +task0_model0_perelson1_V condition1 89000 7 sd_task0_model0_perelson1_V From 1ab59c7dbd58ad52a9aee18fa9df8ca378c9fe4b Mon Sep 17 00:00:00 2001 From: Maren Philipps Date: Sat, 19 Sep 2026 13:45:45 +0000 Subject: [PATCH 2/4] Perelson_Science1996: rename observable to patient_105 The problem implements the plasma HIV-1 RNA time course of a single patient, patient 105 of the publication (ID 4 in the data release). Name the observable accordingly instead of task0_model0_perelson1_V, which carried no information about what is measured. PEtab requires the noise placeholder to be spelled noiseParameter${n}_${observableId}, so the noise parameter is renamed along with it: sd_task0_model0_perelson1_V -> sd_patient_105. The id does not occur in the SBML model. petablint passes. Co-Authored-By: Claude Opus 5 --- .../measurementData_Perelson_Science1996.tsv | 32 +++++++++---------- .../observables_Perelson_Science1996.tsv | 2 +- .../parameters_Perelson_Science1996.tsv | 2 +- .../simulatedData_Perelson_Science1996.tsv | 32 +++++++++---------- ...tionSpecification_Perelson_Science1996.tsv | 2 +- 5 files changed, 35 insertions(+), 35 deletions(-) diff --git a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv index 2094a65d..36d421f9 100644 --- a/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv @@ -1,17 +1,17 @@ observableId preequilibrationConditionId simulationConditionId measurement time observableParameters noiseParameters -task0_model0_perelson1_V condition1 1022000 0 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1079000 0.083 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2148000 0.166 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1043000 0.25 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2170000 0.5 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1538000 0.75 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 3098000 1 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 2262000 1.25 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1315000 1.5 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 1200000 1.75 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 980000 2 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 710000 3 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 500000 4 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 280000 5 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 170000 6 sd_task0_model0_perelson1_V -task0_model0_perelson1_V condition1 89000 7 sd_task0_model0_perelson1_V +patient_105 condition1 1022000 0 sd_patient_105 +patient_105 condition1 1079000 0.083 sd_patient_105 +patient_105 condition1 2148000 0.166 sd_patient_105 +patient_105 condition1 1043000 0.25 sd_patient_105 +patient_105 condition1 2170000 0.5 sd_patient_105 +patient_105 condition1 1538000 0.75 sd_patient_105 +patient_105 condition1 3098000 1 sd_patient_105 +patient_105 condition1 2262000 1.25 sd_patient_105 +patient_105 condition1 1315000 1.5 sd_patient_105 +patient_105 condition1 1200000 1.75 sd_patient_105 +patient_105 condition1 980000 2 sd_patient_105 +patient_105 condition1 710000 3 sd_patient_105 +patient_105 condition1 500000 4 sd_patient_105 +patient_105 condition1 280000 5 sd_patient_105 +patient_105 condition1 170000 6 sd_patient_105 +patient_105 condition1 89000 7 sd_patient_105 diff --git a/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv index cec307df..355c5c12 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 noiseParameter1_patient_105 log10 diff --git a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv index d11162a5..d99c42f0 100644 --- a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv @@ -4,4 +4,4 @@ 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 +sd_patient_105 sigma_patient_105 log10 1E-10 10000000000 100000 1 diff --git a/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv index e68cd7a8..e5255382 100644 --- a/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv @@ -1,17 +1,17 @@ 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 +patient_105 condition1 1860000 0 1 log10 normal +patient_105 condition1 1979067.19116451 0.105 1 log10 normal +patient_105 condition1 2036200.82349015 0.169 1 log10 normal +patient_105 condition1 2110605.53393048 0.282 1 log10 normal +patient_105 condition1 2170386.58687064 0.492 1 log10 normal +patient_105 condition1 2133014.95269381 0.757 1 log10 normal +patient_105 condition1 2007892.04532337 1.029 1 log10 normal +patient_105 condition1 1866529.16992371 1.253 1 log10 normal +patient_105 condition1 1667316.74545722 1.533 1 log10 normal +patient_105 condition1 1508578.7515073 1.75 1 log10 normal +patient_105 condition1 1304224.95565882 2.038 1 log10 normal +patient_105 condition1 745147.675932322 3.013 1 log10 normal +patient_105 condition1 404793.605719818 3.987 1 log10 normal +patient_105 condition1 214718.67877913 4.968 1 log10 normal +patient_105 condition1 111378.497695453 5.972 1 log10 normal +patient_105 condition1 57731.3816378903 6.973 1 log10 normal diff --git a/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv index 262254a7..6ceca44b 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 [days] lin RNA copies/ml lin HIV-mRNA copies - patient_105 patient_105 time 0 0 From 947445f3adffa33cc2755206985d19413965d34c Mon Sep 17 00:00:00 2001 From: Maren Philipps Date: Sat, 19 Sep 2026 13:56:08 +0000 Subject: [PATCH 3/4] Perelson_Science1996: name species and parameters after the paper Rename four ids to the notation of the publication: Vin -> V_I infectious virions, present before the drug took effect Vni -> V_NI non-infectious virions, produced after the drug took effect NN -> N virions produced per infected cell over its lifetime K0 -> k infection rate constant Vin/Vni were particularly easy to misread, being one transposition apart and resembling the paper's V_NI. c, delta, T0 and V already matched. Tstar keeps its name, since SBML ids cannot contain "*". Only the id, species and ci occurrences were edited, so the model file keeps its formatting. The right-hand sides of the ODE system, assembled from the stoichiometry and the kinetic laws, are bit-identical to those of the previous model under the renaming; libsbml reports no new warnings and petablint passes. Co-Authored-By: Claude Opus 5 --- .../model_Perelson_Science1996.xml | 26 +++++++++---------- .../parameters_Perelson_Science1996.tsv | 4 +-- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml b/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml index 88a16256..b9e68a8d 100644 --- a/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml +++ b/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml @@ -39,15 +39,15 @@ - - + + - + - + @@ -55,15 +55,15 @@ - + - K0 + k T0 - Vin + V_I @@ -85,14 +85,14 @@ - + c - Vin + V_I @@ -100,14 +100,14 @@ - + c - Vni + V_NI @@ -115,7 +115,7 @@ - + @@ -125,7 +125,7 @@ delta - NN + N Tstar diff --git a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv index d99c42f0..33f9a233 100644 --- a/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv +++ b/Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv @@ -1,7 +1,7 @@ parameterId parameterName parameterScale lowerBound upperBound nominalValue estimate -NN NN log10 1E-05 100000 480.000000026982 0 +N N 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 +k k log10 1E-05 100000 3.89999999976204E-07 0 sd_patient_105 sigma_patient_105 log10 1E-10 10000000000 100000 1 From f28d0083168b5dbd9f1a67d2b441650dcd626e42 Mon Sep 17 00:00:00 2001 From: Maren Philipps Date: Sat, 19 Sep 2026 15:23:13 +0000 Subject: [PATCH 4/4] Perelson_Science1996: fit on the paper's time base and estimate V0 Three changes to make the problem the one solved in the publication, plus the bookkeeping they entail. Time base. Eq. 6 of the paper is not a function of time since dosing: its t = 0 is the onset of the drug effect, which for patient 105 follows dosing by a measured pharmacologic delay of 6 h. Subtract 0.25 d from the sampling times and drop the three measurements taken before the onset, which lie outside the model's domain (16 -> 13 measurements). Initial conditions. V0 becomes an estimated parameter, as the Fig. 1 caption states it was in the publication. The other initial values follow from the pre-treatment quasi steady state and are now SBML initial assignments rather than literals, so they track V0, c and delta: V_I(0) = V0, V(0) = V0 and Tstar(0) = c*V0/(N_virions*delta). Steady-state coupling of k. The same steady state gives c = N_virions*k*T0, so k is an initial assignment c/(N_virions*T0) and no longer a free parameter. This is required, not cosmetic: with k held at a literal the model drifts from Eq. 6 as soon as the estimated c moves -- by up to 13 % at c = 1.5 -- so V0 would have been estimated against the wrong model. Noise. The observable's noiseFormula is the literal 1 and the sd_ parameter is gone, making the objective the unweighted sum of squared log10 residuals, which is the nonlinear least squares problem the paper solves. The noiseParameters column is dropped from both measurement and simulation tables. The burst size is named N_virions rather than N because N is SymPy's numeric-evaluation function, and AMICI parses SBML initial assignments through SymPy; a parameter called N makes the model fail to import. Validation: AMICI at the nominal parameters agrees with Eq. 6 to 2e-8, and re-estimating recovers the values published for patient 105 -- c = 2.062, delta = 0.525, V0 = 1.832e6 against 2.06, 0.53, 1.86e6 -- so the nominal parameters now sit essentially at the optimum. simulatedData is regenerated from the model at the nominal parameters, which also fixes the measurement/ simulation row mismatch that made petab.v1.calculate.calculate_llh raise IndexError for this problem (see #278); calculate_llh and AMICI now agree on -199.60959155685103. petablint passes, README updated with bmp-create-overview --update. Co-Authored-By: Claude Opus 5 --- .../Perelson_Science1996/README.md | 149 ++++++++++++++++++ .../measurementData_Perelson_Science1996.tsv | 31 ++-- .../model_Perelson_Science1996.xml | 55 ++++++- .../observables_Perelson_Science1996.tsv | 2 +- .../parameters_Perelson_Science1996.tsv | 5 +- .../simulatedData_Perelson_Science1996.tsv | 31 ++-- ...tionSpecification_Perelson_Science1996.tsv | 2 +- README.md | 2 +- 8 files changed, 231 insertions(+), 46 deletions(-) create mode 100644 Benchmark-Models/Perelson_Science1996/README.md 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 36d421f9..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 -patient_105 condition1 1022000 0 sd_patient_105 -patient_105 condition1 1079000 0.083 sd_patient_105 -patient_105 condition1 2148000 0.166 sd_patient_105 -patient_105 condition1 1043000 0.25 sd_patient_105 -patient_105 condition1 2170000 0.5 sd_patient_105 -patient_105 condition1 1538000 0.75 sd_patient_105 -patient_105 condition1 3098000 1 sd_patient_105 -patient_105 condition1 2262000 1.25 sd_patient_105 -patient_105 condition1 1315000 1.5 sd_patient_105 -patient_105 condition1 1200000 1.75 sd_patient_105 -patient_105 condition1 980000 2 sd_patient_105 -patient_105 condition1 710000 3 sd_patient_105 -patient_105 condition1 500000 4 sd_patient_105 -patient_105 condition1 280000 5 sd_patient_105 -patient_105 condition1 170000 6 sd_patient_105 -patient_105 condition1 89000 7 sd_patient_105 +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 b9e68a8d..f50551b0 100644 --- a/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml +++ b/Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml @@ -37,18 +37,61 @@ - - - + + + - + - + + + + + + + + c + + + N_virions + T0 + + + + + + + + + + + c + V0 + + + + N_virions + delta + + + + + + + V0 + + + + + V0 + + + @@ -125,7 +168,7 @@ delta - N + N_virions Tstar diff --git a/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv b/Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv index 355c5c12..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 -patient_105 V normal noiseParameter1_patient_105 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 33f9a233..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 -N N 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 -k k log10 1E-05 100000 3.89999999976204E-07 0 -sd_patient_105 sigma_patient_105 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 e5255382..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 -patient_105 condition1 1860000 0 1 log10 normal -patient_105 condition1 1979067.19116451 0.105 1 log10 normal -patient_105 condition1 2036200.82349015 0.169 1 log10 normal -patient_105 condition1 2110605.53393048 0.282 1 log10 normal -patient_105 condition1 2170386.58687064 0.492 1 log10 normal -patient_105 condition1 2133014.95269381 0.757 1 log10 normal -patient_105 condition1 2007892.04532337 1.029 1 log10 normal -patient_105 condition1 1866529.16992371 1.253 1 log10 normal -patient_105 condition1 1667316.74545722 1.533 1 log10 normal -patient_105 condition1 1508578.7515073 1.75 1 log10 normal -patient_105 condition1 1304224.95565882 2.038 1 log10 normal -patient_105 condition1 745147.675932322 3.013 1 log10 normal -patient_105 condition1 404793.605719818 3.987 1 log10 normal -patient_105 condition1 214718.67877913 4.968 1 log10 normal -patient_105 condition1 111378.497695453 5.972 1 log10 normal -patient_105 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 6ceca44b..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 - patient_105 patient_105 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) |