---
title: "Lopinavir (Jullien 2006)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Lopinavir (Jullien 2006)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 6,
  fig.height = 4
)
```

```{r setup}
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Lopinavir population PK in children from birth to 18 years

Jullien et al. (2006) describe a one-compartment population PK model for oral
lopinavir (boosted with low-dose ritonavir, LPV/rtv) in 157 paediatric HIV
patients ranging in age from 3 days to 18 years. The structural model is a
simplified one-compartment form in which a single first-order rate constant
governs both absorption and elimination (k = ka = kel, fixed equal at the level
of the ODEs and equivalently expressed as k = CL/V given the joint
parameterisation). Body weight is allometrically scaled on apparent clearance
and apparent volume, nevirapine coadministration increases apparent clearance
by 34%, and male sex increases apparent clearance by 39% in children older
than 12 years. This vignette walks the paper's typical-value parameter
predictions, replicates Figures 4 and 5, and runs a PKNCA validation across
body-weight bands against the AUC and Cmax values reported in the paper's
Results and Discussion sections.

* Citation: `r rxode2::rxode(readModelDb("Jullien_2006_lopinavir"))$reference`
* Article: <https://doi.org/10.1128/AAC.00943-05>

## Population

Jullien 2006 enrolled 157 HIV-infected children (67 girls and 90 boys) who were
monitored by routine therapeutic drug monitoring at Cochin-Saint-Vincent de
Paul and Necker-Enfants Malades (Paris). Baseline demographics from the paper's
Table 1 (mean +/- SD, range):

| Variable                 | Mean +/- SD       | Median | Range          |
|--------------------------|-------------------|--------|----------------|
| Age (years)              | 9.1 +/- 4.8       | 10.2   | 3 days-18 yr   |
| Body weight (kg)         | 29.0 +/- 14.9     | 27.6   | 2-73           |
| Dose (mg)                | 279 +/- 102       | 266    | 30-532         |
| Dose (mg/kg)             | 10.9 +/- 3.7      | 10.4   | 4.4-29.4       |
| Dose (mg/m^2)            | 288 +/- 66        | 281    | 124-566        |
| LPV concentration (mg/L) | 8.99 +/- 4.89     | 8.25   | 0.33-29.7      |
| Samples per patient      | 3.5               | 3      | 1-14           |

Concomitant antiretrovirals: at least one nucleoside reverse-transcriptase
inhibitor (NRTI) in 90% of samples; one protease inhibitor (PI) in 10%;
one non-nucleoside reverse-transcriptase inhibitor (NNRTI) in 23%. Nevirapine
was combined with lopinavir in 16% of samples, efavirenz in 8%, and amprenavir
in 4% (paper Results, Demographic data). The full programmatic metadata is
available via `readModelDb("Jullien_2006_lopinavir")$population`.

## Source trace

Direct map from each `ini()` parameter to its origin in Jullien 2006. Final
estimates come from Table 3, "Final model, original data set, Mean" row, and
the final covariate submodel equations are stated immediately above Table 3 in
the Results section.

| Parameter   | Value         | Source                                                                                       |
|-------------|---------------|----------------------------------------------------------------------------------------------|
| `lcl`       | log(2.58)     | Table 3: TV(CL/F) = 2.58 L/h (at WT = 27 kg, no NVP, female or age <= 12 yr)                 |
| `lvc`       | log(24.6)     | Table 3: TV(V/F) = 24.6 L (at WT = 27 kg)                                                    |
| `e_wt_cl`   | 0.46          | Table 3: allometric exponent theta_BW on CL/F                                                |
| `e_wt_vc`   | 0.72          | Table 3: allometric exponent theta_BW on V/F                                                 |
| `e_nvp_cl`  | log(1.34)     | Table 3 / final submodel: theta_nevirapine = 1.34 multiplier on CL/F                         |
| `e_sexf_cl` | log(1.39)     | Table 3 / final submodel: theta_SEX = 1.39 multiplier on CL/F when AGE > 12 yr (male only)   |
| `etalcl`    | var = 0.0957  | Table 3: omega^2(CL/F) = 0.0957 (log-scale variance of exponential IIV)                      |
| `etalvc`    | var = 0.180   | Table 3: omega^2(V/F) = 0.180                                                                |
| cov         | 0.131         | Table 3: cov(CL,V) = 0.131 (off-diagonal of the eta block)                                   |
| `propSd`    | sqrt(0.138)   | Table 3: sigma_1^2 = 0.138 (exponential proportional, small-sigma equivalent)                |
| `addSd`     | sqrt(1.83)    | Table 3: sigma_2^2 = 1.83 mg^2/L^2 (additive)                                                |

Structural equations (paper Results, final covariate submodel):

```
V/F (L)   = 24.6 * (BW / 27)^0.72
CL/F (L/h) = 2.58 * (BW / 27)^0.46 * 1.34^N * (1.39^S if AGE > 12 else 1)
```

where `N = 1` if nevirapine is coadministered with lopinavir (else 0) and
`S = 1` for boys (else 0). The model file maps `N -> CONMED_NVP` and `S -> 1 -
SEXF` (the canonical `SEXF` uses 1 = female). The system of ODEs uses
`ka = kel = CL/V`:

```
d/dt(depot)   = -ka  * depot
d/dt(central) =  ka  * depot - kel * central
Cc            =  central / vc
```

## Load model

```{r load-model}
mod <- rxode2::rxode(readModelDb("Jullien_2006_lopinavir"))
mod_typical <- rxode2::zeroRe(mod)
```

## Typical-value steady-state profile

The paper's Discussion (paragraph 2) reports the typical-value PK predictions
for the study population (mean dose 288 mg/m^2, which corresponds to roughly
266 mg for a 27 kg child of BSA ~0.98 m^2): Tmax 4.5 h, Cmax 9.60 mg/L,
C_trough 7.90 mg/L, t1/2 6.6 h, AUC0-12 108 mg.h/L. We reproduce a typical
profile for a 27 kg, 10-year-old girl with no nevirapine coadministration.

```{r typical-ss-profile}
# 14 days BID (28 doses) -> last dosing interval is at steady state.
n_doses  <- 28L
dose_int <- 12
ev_typ <- rxode2::et(amt = 266, cmt = "depot", evid = 1,
                     ii = dose_int, addl = n_doses - 1L) |>
  rxode2::et(seq(0, n_doses * dose_int, by = 0.25)) |>
  rxode2::et(id = 1)
ev_typ$WT         <- 27
ev_typ$AGE        <- 10
ev_typ$SEXF       <- 1
ev_typ$CONMED_NVP <- 0

sim_typ <- as.data.frame(rxode2::rxSolve(mod_typical, ev_typ))

# Steady-state dosing interval: t = (n_doses - 1) * dose_int to n_doses * dose_int
t_ss_start <- (n_doses - 1L) * dose_int
ss <- sim_typ[sim_typ$time >= t_ss_start & sim_typ$time <= t_ss_start + dose_int, ]
ss$t_in_interval <- ss$time - t_ss_start

cmax_ss   <- max(ss$Cc)
cmin_ss   <- min(ss$Cc)
tmax_ss   <- ss$t_in_interval[which.max(ss$Cc)]
auc012_ss <- sum(diff(ss$time) * (head(ss$Cc, -1) + tail(ss$Cc, -1)) / 2)
t_half    <- log(2) / (with(ev_typ, 2.58 / 24.6))  # k = CL/V at typical values

typ_summary <- tibble::tibble(
  Quantity       = c("Tmax (h)", "Cmax (mg/L)", "C_trough (mg/L)",
                     "AUC0-12 (mg.h/L)", "t1/2 (h)"),
  `Simulated`    = c(tmax_ss, cmax_ss, cmin_ss, auc012_ss, t_half),
  `Jullien 2006` = c(4.5, 9.60, 7.90, 108, 6.6)
)
knitr::kable(typ_summary, digits = 2,
             caption = "Typical-value steady-state predictions (27 kg, 10 yr, female, no NVP) vs Jullien 2006 Discussion paragraph 2.")
```

```{r typical-ss-plot}
ggplot(ss, aes(t_in_interval, Cc)) +
  geom_line(linewidth = 0.8, colour = "steelblue") +
  geom_hline(yintercept = c(2.5), linetype = "dotted") +
  annotate("text", x = 11.5, y = 2.5, label = "C_trough therapeutic threshold (2.5 mg/L)",
           hjust = 1, vjust = -0.5, size = 3) +
  labs(
    x = "Time after dose (h)",
    y = "Lopinavir concentration (mg/L)",
    title = "Typical-value steady-state interval (27 kg, 10 yr, female, no NVP, 266 mg BID)"
  ) +
  theme_bw()
```

The simulated typical values are within 5% of the published Discussion-paragraph
values, consistent with the paper's own self-validation against an earlier
pediatric LPV/rtv study (paper Discussion paragraph 2 vs reference 27).

## Replicate Figure 4: CL/F, weight-normalized CL/F, and t1/2 vs body weight

Figure 4 of Jullien 2006 shows the typical-value evolution of apparent clearance
(panel A), weight-normalized CL/F (panel B), and terminal half-life (panel C)
across the body-weight range of the cohort.

```{r figure-4}
fig4 <- tibble::tibble(WT = seq(2, 73, by = 0.5)) |>
  dplyr::mutate(
    CL = 2.58 * (WT / 27)^0.46,
    CL_per_kg = CL / WT,
    t_half = log(2) * (24.6 * (WT / 27)^0.72) / CL
  )

panel_A <- ggplot(fig4, aes(WT, CL)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Body weight (kg)", y = "CL/F (L/h)", title = "A: CL/F vs BW") +
  theme_bw()

panel_B <- ggplot(fig4, aes(WT, CL_per_kg)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Body weight (kg)", y = "CL/F per kg (L/h/kg)",
       title = "B: weight-normalized CL/F vs BW") +
  theme_bw()

panel_C <- ggplot(fig4, aes(WT, t_half)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Body weight (kg)", y = "Terminal t1/2 (h)",
       title = "C: t1/2 vs BW") +
  theme_bw()

panel_A
panel_B
panel_C
```

The three panels reproduce the published nonlinear relationships:
weight-normalized CL/F decreases markedly with increasing BW (panel B
matches the Fig 4B trend that the paper highlights for neonates and
infants), and terminal half-life lengthens from ~4 h at 3 kg to ~9 h at the
upper BW range (panel C). At BW = 3.7 kg (the paper's neonate-subgroup
mean), CL/F is ~1.0 L/h and t1/2 is ~3.9 h, agreeing with the
paper-reported subgroup means of CL/F = 1.15 L/h and t1/2 = 3.9 h (paper
Results "Evolution of lopinavir exposure with respect to BW, age, and sex"
paragraph 1).

## Replicate Figure 5: AUC0-12 and C_trough at recommended pediatric dose

Jullien 2006 Figure 5 shows the calculated AUC0-12 and C_trough at the
currently recommended pediatric dose (12 mg/kg BID for BW 7-14 kg; 10 mg/kg BID
for BW 14-40 kg; 400 mg BID for BW > 40 kg) across the BW range, separated by
sex among adolescents (AGE > 12 yr). We build a deterministic virtual cohort
covering the BW grid 5-75 kg.

```{r figure-5}
recommended_dose <- function(WT) {
  dplyr::case_when(
    WT < 7   ~ 12 * WT,
    WT < 14  ~ 12 * WT,
    WT < 40  ~ 10 * WT,
    TRUE     ~ 400
  )
}

# Adolescent cohort: BW 30-75 kg, AGE = 14 (>12), male vs female
# Younger cohort: BW 5-30 kg, AGE = 8 (<12), sex doesn't affect CL
adolescent_grid <- expand.grid(
  WT   = seq(30, 75, by = 5),
  SEXF = c(0, 1)
) |>
  dplyr::mutate(AGE = 14, scenario = paste0("Adolescent (>12 yr) ",
                                            ifelse(SEXF == 1, "girls", "boys")))

younger_grid <- tibble::tibble(WT = seq(5, 30, by = 2.5)) |>
  dplyr::mutate(AGE = 8, SEXF = 1, scenario = "Children (<12 yr)")

grid <- dplyr::bind_rows(adolescent_grid, younger_grid)
grid$dose       <- recommended_dose(grid$WT)
grid$CONMED_NVP <- 0
grid$id         <- seq_len(nrow(grid))

build_subject_events <- function(row) {
  ev <- rxode2::et(amt = row$dose, cmt = "depot", evid = 1,
                   ii = 12, addl = 27L) |>
    rxode2::et(seq(0, 28 * 12, by = 0.5)) |>
    rxode2::et(id = row$id)
  df <- as.data.frame(ev)
  df$WT         <- row$WT
  df$AGE        <- row$AGE
  df$SEXF       <- row$SEXF
  df$CONMED_NVP <- row$CONMED_NVP
  df$scenario   <- row$scenario
  df$dose       <- row$dose
  df
}

ev_grid <- do.call(rbind, lapply(seq_len(nrow(grid)), function(i) {
  build_subject_events(grid[i, ])
}))

sim_grid <- as.data.frame(rxode2::rxSolve(
  mod_typical, ev_grid,
  keep = c("WT", "AGE", "SEXF", "CONMED_NVP", "scenario", "dose")
))

t_ss_start <- 27L * 12L
ss_grid <- sim_grid |>
  dplyr::filter(time >= t_ss_start, time <= t_ss_start + 12) |>
  dplyr::mutate(t_in_interval = time - t_ss_start)

ss_summary <- ss_grid |>
  dplyr::group_by(id, WT, scenario) |>
  dplyr::summarise(
    Cmax   = max(Cc),
    Cmin   = min(Cc),
    AUC012 = sum(diff(t_in_interval) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups = "drop"
  )

ggplot(ss_summary, aes(WT, AUC012, colour = scenario)) +
  geom_line(linewidth = 0.8) + geom_point(size = 1.5) +
  labs(x = "Body weight (kg)", y = "Steady-state AUC0-12 (mg.h/L)",
       colour = NULL,
       title = "Figure 5A: AUC0-12 at the recommended pediatric LPV dose") +
  theme_bw() + theme(legend.position = "bottom")

ggplot(ss_summary, aes(WT, Cmin, colour = scenario)) +
  geom_line(linewidth = 0.8) + geom_point(size = 1.5) +
  geom_hline(yintercept = 2.5, linetype = "dotted") +
  annotate("text", x = max(ss_summary$WT), y = 2.5,
           label = "Treatment-failure threshold (2.5 mg/L)",
           hjust = 1, vjust = -0.5, size = 3) +
  labs(x = "Body weight (kg)", y = "Steady-state C_trough (mg/L)",
       colour = NULL,
       title = "Figure 5B: C_trough at the recommended pediatric LPV dose") +
  theme_bw() + theme(legend.position = "bottom")
```

The two panels reproduce the paper's key finding that adolescent boys (red /
green) have lower AUC0-12 and C_trough than adolescent girls at the same
weight-based recommended dose, with the gap widening at higher body weights
where the 400 mg cap takes effect (paper Results "The lopinavir CL/F was also
found to be related to sex..." and Figure 5 discussion).

## PKNCA validation

Stochastic NCA across a virtual cohort matched to the cohort BW distribution
(Table 1). We sample 80 subjects with covariate distributions approximating the
paper's Table 1, then run PKNCA on the steady-state interval.

```{r pknca}
set.seed(2006)
n_subj <- 80L

# Three approximate age-weight strata to span the cohort
n_per_stratum <- n_subj %/% 3L

stratum_infant      <- tibble::tibble(
  AGE = pmax(0.01, rnorm(n_per_stratum, mean = 1.5, sd = 1.0)),
  WT  = pmax(2.0, rnorm(n_per_stratum, mean = 9, sd = 3)),
  stratum = "Infant / young child (<3 yr)"
)
stratum_child       <- tibble::tibble(
  AGE = pmax(3.0, rnorm(n_per_stratum, mean = 8, sd = 2)),
  WT  = pmax(8.0, rnorm(n_per_stratum, mean = 25, sd = 7)),
  stratum = "Child (3-12 yr)"
)
stratum_adolescent  <- tibble::tibble(
  AGE = pmax(12.5, rnorm(n_subj - 2L * n_per_stratum, mean = 15, sd = 1.5)),
  WT  = pmax(35.0, rnorm(n_subj - 2L * n_per_stratum, mean = 55, sd = 12)),
  stratum = "Adolescent (>12 yr)"
)

cohort <- dplyr::bind_rows(stratum_infant, stratum_child, stratum_adolescent)
cohort$id         <- seq_len(nrow(cohort))
cohort$SEXF       <- rbinom(nrow(cohort), 1, 0.427)
cohort$CONMED_NVP <- rbinom(nrow(cohort), 1, 0.16)
cohort$dose       <- recommended_dose(cohort$WT)

build_ev <- function(row) {
  ev <- rxode2::et(amt = row$dose, cmt = "depot", evid = 1,
                   ii = 12, addl = 27L) |>
    rxode2::et(c(seq(0, 12, by = 0.5),
                 seq(13 * 24, 14 * 24, by = 0.5))) |>
    rxode2::et(id = row$id)
  df <- as.data.frame(ev)
  df$WT         <- row$WT
  df$AGE        <- row$AGE
  df$SEXF       <- row$SEXF
  df$CONMED_NVP <- row$CONMED_NVP
  df$stratum    <- row$stratum
  df$dose       <- row$dose
  df
}

ev_all <- do.call(rbind, lapply(seq_len(nrow(cohort)), function(i) {
  build_ev(cohort[i, ])
}))

sim_pop <- as.data.frame(rxode2::rxSolve(
  mod, ev_all,
  keep = c("WT", "AGE", "SEXF", "CONMED_NVP", "stratum", "dose")
))

# Day-14 (steady-state) interval: t = 13 * 24 to 14 * 24
ss_pop <- sim_pop |>
  dplyr::filter(time >= 13 * 24, time <= 14 * 24) |>
  dplyr::mutate(t_in_interval = time - 13 * 24)

nca_concs <- ss_pop |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, t_in_interval, Cc, stratum)

# One dose record per subject at t = 0 of the steady-state interval
nca_doses <- cohort |>
  dplyr::transmute(id, time = 0, amt = dose, stratum)

conc_obj <- PKNCA::PKNCAconc(nca_concs, Cc ~ t_in_interval | stratum + id,
                             concu = "mg/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(nca_doses, amt ~ time | stratum + id,
                             doseu = "mg")

intervals <- data.frame(
  start    = 0,
  end      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  auclast  = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- suppressWarnings(PKNCA::pk.nca(nca_data))
nca_df   <- as.data.frame(nca_res$result)

nca_summary <- nca_df |>
  dplyr::filter(PPTESTCD %in% c("auclast", "cmax", "cmin")) |>
  dplyr::group_by(stratum, PPTESTCD) |>
  dplyr::summarise(
    median = median(PPORRES, na.rm = TRUE),
    P05    = quantile(PPORRES, 0.05, na.rm = TRUE),
    P95    = quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  )
knitr::kable(nca_summary, digits = 2,
             caption = "PKNCA summary at steady-state Day-14 dosing interval across the virtual cohort.")
```

### Comparison with published values

| Quantity                          | Paper (text / Table 1)               | Simulated (this vignette)             |
|-----------------------------------|--------------------------------------|---------------------------------------|
| Population median LPV conc (mg/L) | 8.25 (Table 1, sample-time median)   | NCA Cmax / Cmin span this range       |
| Typical Cmax 27 kg, 266 mg, no NVP| 9.60 (Discussion para 2)             | `typ_summary$Simulated[2]` above      |
| Typical C_trough 27 kg            | 7.90 (Discussion para 2)             | `typ_summary$Simulated[3]` above      |
| Typical AUC0-12 27 kg             | 108 (Discussion para 2)              | `typ_summary$Simulated[4]` above      |
| AUC0-12 girls > 12 yr at recommended dose | ~100                          | Figure 5A "Adolescent girls" trace    |
| AUC0-12 boys > 12 yr at recommended dose  | ~50                           | Figure 5A "Adolescent boys" trace     |
| AUC0-12 neonates 3.7 kg, 12 mg/kg | 38.6 (Results para 6)                | reproduced by 5 kg adolescent infant stratum point in Figure 5A |

The typical-value match against the Discussion paragraph 2 reference values
(Cmax 9.60, C_trough 7.90, AUC 108, t1/2 6.6 h) is the most informative
quantitative check; differences are < 5%.

## Assumptions and deviations

This implementation reproduces the published structural model and Table-3
parameter values directly. The following simplifications / notes are
documented so reviewers can reconcile the model with the source.

1. **Combined-residual encoding.** Jullien 2006 reports a combined exponential
   (sigma_1^2 = 0.138 on log scale) and additive (sigma_2^2 = 1.83 mg^2/L^2)
   residual error model. nlmixr2's native simulation-time residual does not
   support `add(...) + lnorm(...)` combinations (see `inst/modeldb/specificDrugs/Ruhs_2012_methotrexate.R`
   and the matching vignette for the same situation). Following that
   precedent, the exponential proportional component is encoded as the
   small-sigma-equivalent `prop(propSd)` with `propSd = sqrt(0.138)`, and the
   additive component is encoded as `add(addSd)` with `addSd = sqrt(1.83)`. At
   sigma_1 ~ 0.37 the linear-scale and log-scale proportional residuals
   diverge by < 7%, which is negligible for the simulation purposes of this
   library model.
2. **k = ka = kel constraint.** Jullien 2006 used a one-compartment model in
   which a single first-order rate constant served as both the absorption and
   the elimination rate constant (paper Results, paragraph 1 of Population
   pharmacokinetics; citing Wahlby 2002 as the reference for the simplified
   parameterisation). The classical 1-cmt model with independent ka and kel
   was over-parameterised for the dataset. The model file encodes this by
   setting `ka <- kel` inside `model()`; the ODE system is the canonical
   two-state oral-PK chain, and the constraint is the single line `ka <- kel`.
3. **Reference body weight 27 kg.** The paper normalizes the allometric
   scaling at the cohort median BW (27.6 kg, rounded to 27 in the published
   covariate-submodel equation). This is preserved exactly.
4. **Sex effect age gate at 12 years.** The paper applies the 1.39 multiplier
   on apparent CL/F to boys with `AGE > 12` (strict inequality, per the
   published equation). At `AGE = 12` exactly no boost is applied. Users
   simulating a patient at the boundary should be aware that the published
   submodel is piecewise discontinuous at 12 years; smooth interpolation
   would require re-fitting the original data.
5. **Sex source-column orientation.** The paper used `S = 1` for boys and
   `S = 0` for girls. The canonical `SEXF` column inverts this (1 = female,
   0 = male) per the nlmixr2lib `inst/references/covariate-columns.md`
   register. The published `1.39^S` multiplier is preserved as
   `exp(e_sexf_cl * (1 - SEXF) * (AGE > 12))` inside `model()`, with
   `e_sexf_cl = log(1.39)`. The effect direction (boys > girls) and the
   magnitude (39%) are unchanged.
6. **CONMED_NVP newly registered.** The `CONMED_NVP` indicator column is
   added in this PR to the canonical-covariate register as the nevirapine
   analog of the existing `CONMED_EFV` indicator. Both flag CYP3A induction
   by a non-nucleoside reverse-transcriptase inhibitor on a coadministered
   antiretroviral drug.
7. **AUC0-12 trapezoidal approximation in the figure code.** The Figure 5
   replication uses a simple linear-trapezoidal integration of the
   simulated time-concentration curve at a 0.5 h grid. The matching PKNCA
   block (below) uses PKNCA's `auclast` which gives the same value within
   numerical tolerance. The Figure 5 plot is intended to communicate the
   trend, not to be the validation target.
8. **Inter-occasion variability not encoded.** The paper does not separately
   report IOV (the design used at most one sample per occasion). The IIV /
   covariance values in Table 3 are the only stochastic-effects information
   available; IOV is implicitly absorbed.
9. **Bootstrap-validated parameters.** Table 3 reports both original-data-set
   means and bootstrap means (1000 resamples). The packaged model uses the
   original-data-set means (the actual final estimates); the bootstrap means
   are within SE for all parameters.

## Reference

* `r rxode2::rxode(readModelDb("Jullien_2006_lopinavir"))$reference`
