Skip to content

Perelson_Science1996 data and model fixes #323

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
m-philipps wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Perelson_Science1996 data and model fixes #323

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
cd82660
Perelson_Science1996: use Wu data release instead of digitised data (…
m-philipps Sep 19, 2026
1ab59c7
Perelson_Science1996: rename observable to patient_105
m-philipps Sep 19, 2026
947445f
Perelson_Science1996: name species and parameters after the paper
m-philipps Sep 19, 2026
f28d008
Perelson_Science1996: fit on the paper's time base and estimate V0
m-philipps Sep 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 149 additions & 0 deletions Benchmark-Models/Perelson_Science1996/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
31 changes: 14 additions & 17 deletions Benchmark-Models/Perelson_Science1996/measurementData_Perelson_Science1996.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
73 changes: 58 additions & 15 deletions Benchmark-Models/Perelson_Science1996/model_Perelson_Science1996.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,33 +37,76 @@
<compartment metaid="metaid_0" sboTerm="SBO:0000410" id="default" size="1" constant="true"/>
</listOfCompartments>
<listOfSpecies>
<species metaid="metaid_1" id="Tstar" compartment="default" initialConcentration="15061.32075" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_2" id="V" compartment="default" initialConcentration="1860000" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_3" id="Vin" compartment="default" initialConcentration="1860000" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_4" id="Vni" compartment="default" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_1" id="Tstar" compartment="default" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_2" id="V" compartment="default" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_3" id="V_I" compartment="default" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species metaid="metaid_4" id="V_NI" compartment="default" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
</listOfSpecies>
<listOfParameters>
<parameter metaid="metaid_5" id="NN" value="480" constant="true"/>
<parameter metaid="metaid_5" id="N_virions" value="480" constant="true"/>
<parameter metaid="metaid_6" id="T0" value="11000" constant="true"/>
<parameter metaid="metaid_7" id="c" value="2.06" constant="true"/>
<parameter metaid="metaid_8" id="delta" value="0.53" constant="true"/>
<parameter metaid="metaid_9" id="K0" value="3.9e-07" constant="true"/>
<parameter metaid="metaid_9" id="k" constant="true"/>
<parameter metaid="metaid_25" id="V0" value="1860000" constant="true"/>
</listOfParameters>
<listOfInitialAssignments>
<initialAssignment symbol="k">
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<divide/>
<ci> c </ci>
<apply>
<times/>
<ci> N_virions </ci>
<ci> T0 </ci>
</apply>
</apply>
</math>
</initialAssignment>
<initialAssignment symbol="Tstar">
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<divide/>
<apply>
<times/>
<ci> c </ci>
<ci> V0 </ci>
</apply>
<apply>
<times/>
<ci> N_virions </ci>
<ci> delta </ci>
</apply>
</apply>
</math>
</initialAssignment>
<initialAssignment symbol="V_I">
<math xmlns="http://www.w3.org/1998/Math/MathML">
<ci> V0 </ci>
</math>
</initialAssignment>
<initialAssignment symbol="V">
<math xmlns="http://www.w3.org/1998/Math/MathML">
<ci> V0 </ci>
</math>
</initialAssignment>
</listOfInitialAssignments>
<listOfReactions>
<reaction metaid="metaid_12" id="v1" reversible="true" fast="false">
<listOfProducts>
<speciesReference metaid="metaid_10" species="Tstar" stoichiometry="1"/>
</listOfProducts>
<listOfModifiers>
<modifierSpeciesReference metaid="metaid_11" species="Vin"/>
<modifierSpeciesReference metaid="metaid_11" species="V_I"/>
</listOfModifiers>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> K0 </ci>
<ci> k </ci>
<ci> T0 </ci>
<ci> Vin </ci>
<ci> V_I </ci>
</apply>
</math>
</kineticLaw>
Expand All @@ -85,37 +128,37 @@
<reaction metaid="metaid_17" id="v3" reversible="true" fast="false">
<listOfReactants>
<speciesReference metaid="metaid_15" species="V" stoichiometry="1"/>
<speciesReference metaid="metaid_16" species="Vin" stoichiometry="1"/>
<speciesReference metaid="metaid_16" species="V_I" stoichiometry="1"/>
</listOfReactants>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> c </ci>
<ci> Vin </ci>
<ci> V_I </ci>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction metaid="metaid_20" id="v4" reversible="true" fast="false">
<listOfReactants>
<speciesReference metaid="metaid_18" species="V" stoichiometry="1"/>
<speciesReference metaid="metaid_19" species="Vni" stoichiometry="1"/>
<speciesReference metaid="metaid_19" species="V_NI" stoichiometry="1"/>
</listOfReactants>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> c </ci>
<ci> Vni </ci>
<ci> V_NI </ci>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction metaid="metaid_24" id="v5" reversible="true" fast="false">
<listOfProducts>
<speciesReference metaid="metaid_21" species="V" stoichiometry="1"/>
<speciesReference metaid="metaid_22" species="Vni" stoichiometry="1"/>
<speciesReference metaid="metaid_22" species="V_NI" stoichiometry="1"/>
</listOfProducts>
<listOfModifiers>
<modifierSpeciesReference metaid="metaid_23" species="Tstar"/>
Expand All @@ -125,7 +168,7 @@
<apply>
<times/>
<ci> delta </ci>
<ci> NN </ci>
<ci> N_virions </ci>
<ci> Tstar </ci>
</apply>
</math>
Expand Down
2 changes: 1 addition & 1 deletion Benchmark-Models/Perelson_Science1996/observables_Perelson_Science1996.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
5 changes: 2 additions & 3 deletions Benchmark-Models/Perelson_Science1996/parameters_Perelson_Science1996.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
31 changes: 14 additions & 17 deletions Benchmark-Models/Perelson_Science1996/simulatedData_Perelson_Science1996.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion Benchmark-Models/Perelson_Science1996/visualizationSpecification_Perelson_Science1996.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading