---
title: "Ritonavir (Kappelhoff 2005)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ritonavir (Kappelhoff 2005)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Kappelhoff_2005_ritonavir"))$reference`
* Description: `r rxode2::rxode(readModelDb("Kappelhoff_2005_ritonavir"))$description`
* Article: <https://doi.org/10.1111/j.1365-2125.2004.02241.x>

The model is the final one-compartment popPK model published by Kappelhoff
et al. (2005) for oral ritonavir in HIV-1-infected adults. The paper retains
a single covariate -- concomitant lopinavir, which multiplies the apparent
oral clearance of ritonavir by 2.72-fold via a power-form effect on `CL/F`.

# Population

Kappelhoff et al. pooled 1228 plasma ritonavir concentrations from 186
ambulatory HIV-1-infected adults followed at the Slotervaart Hospital
outpatient clinic in Amsterdam between January 1999 and June 2003 (Methods,
"Patients" and "Sampling and bioanalysis"). The cohort was predominantly
male (78.5%) and Caucasian (64.5%); demographics in Table 1 summarise
median age 39.4 years (IQR 35.0-46.0), median body weight 71.5 kg
(IQR 63.0-79.8), median baseline CD4 240 cells/mm^3 (IQR 110-380), and
median log10 HIV-1 RNA 4.86 copies/mL (IQR 3.48-5.47).

Of the 186 subjects, 115 (62%) received ritonavir as a booster of another
protease inhibitor at 100 mg QD, 100 mg BID, 133 mg BID, or 200 mg BID;
71 (38%) received ritonavir as a therapeutic antiviral at 300, 400, 500,
600, or 750 mg twice daily. Co-administration breakdown (Table 1):
indinavir/ritonavir n = 49, saquinavir/ritonavir n = 78,
lopinavir/ritonavir n = 36. Per-visit single-time-point samples (n = 505)
were supplemented with 55 full PK profiles (12-15 sampling time points
each); average 3-4 samples per patient over 7-12 months follow-up. All
samples were collected at steady state, at least 2 weeks after initiation
of a ritonavir-containing regimen.

The same demographic information is available programmatically via the
model's `population` metadata:

```{r show-population}
rxode2::rxode(readModelDb("Kappelhoff_2005_ritonavir"))$meta$population[c(
  "n_subjects", "age_range", "weight_range",
  "sex_female_pct", "race_ethnicity"
)]
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in `inst/modeldb/specificDrugs/Kappelhoff_2005_ritonavir.R`.
The table below collects them in one place for review. All values come from
Kappelhoff et al. (2005) Table 2 ("Parameter estimates from the population
pharmacokinetic model and the results of the bootstrap analysis"); the
power-form lopinavir equation is taken from the Results section paragraph
that immediately follows Table 2.

| Equation / parameter | Value | Source location |
|---|---|---|
| `lcl` (log CL/F, no LPV) | log(10.5) | Table 2 row "CL/F (L/h)" = 10.5, RSE 5.55% |
| `lvc` (log V/F) | log(96.6) | Table 2 row "V/F (L)" = 96.6, RSE 10.7% |
| `lka` (log ka) | log(0.871) | Table 2 row "ka (h-1)" = 0.871, RSE 23.1% |
| `ltlag` (log Tlag) | log(0.778) | Table 2 row "Lag-time (h)" = 0.778, RSE 4.91% |
| `e_lpv_cl` (LPV power factor) | 2.72 | Table 2 row "theta_LPV" = 2.72, RSE 11.7% |
| IIV CL/F | 38.3% CV | Table 2 row "IIV CL/F (%)" = 38.3 |
| IIV V/F | 80.0% CV | Table 2 row "IIV V/F (%)" = 80.0 |
| IIV ka | 169% CV | Table 2 row "IIV ka (%)" = 169 |
| Correlation eta_V--eta_ka | 0.868 | Table 2 row "Correlation eta_V - eta_ka" = 0.868 |
| Fraction in P1 | 64.8% | Table 2 row "Fraction in P1 (%)" = 64.8 |
| `addSd_p1` | 0.0600 mg/L | Table 2 row "Additive error P1 (mg/L)" = 0.0600 |
| `addSd_p2` | 0.199 mg/L | Table 2 row "Additive error P2 (mg/L)" = 0.199 |
| `propSd` | 0.154 (15.4%) | Table 2 row "Proportional error (%)" = 15.4 |
| ODE structure | 1-compartment, first-order absorption with lag, first-order elimination | Results paragraph 1 ("best described by a one-compartment model with first-order absorption and elimination ... addition of an absorption lag-time (0.778 h)") |
| `CL/F` covariate equation | `CL/F = 10.5 * 2.72^LPV` | Results equation after Table 2 ("The following equation describes the final model for clearance:") |
| IOV F | 59.1% CV | Table 2 row "IOV F (%)" = 59.1 -- NOT propagated; see Assumptions and deviations |

# Virtual cohort

The original observed data are not publicly available. Three regimens
mirror the dominant dosing scenarios in the Kappelhoff cohort:

* `RTV100_noLPV` -- ritonavir 100 mg BID booster, no lopinavir
  (representative of the indinavir/ritonavir 100/800 mg or
  saquinavir/ritonavir 100/1000 mg booster regimens).
* `RTV100_LPV` -- ritonavir 100 mg BID co-administered with lopinavir
  (the LPV/r 100/400 mg BID regimen tracked by 36 subjects in the
  source dataset).
* `RTV600_noLPV` -- ritonavir 600 mg BID as a therapeutic antiviral
  dose (mid-range of the 300-750 mg BID antiviral cohort).

```{r cohort}
set.seed(2024)

make_cohort <- function(n, dose_mg, lpv, label, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  # Dose every 12 hours for 7 days; observation grid spans the last day at
  # steady state for NCA, with denser sampling around the final peak.
  dose_times <- seq(0, 156, by = 12)
  obs_times <- sort(unique(c(
    seq(0, 168, by = 1),
    seq(144, 168, by = 0.25)
  )))
  dose_rows <- expand.grid(id = ids, time = dose_times) |>
    mutate(amt = dose_mg, evid = 1L, cmt = "depot")
  obs_rows <- expand.grid(id = ids, time = obs_times) |>
    mutate(amt = 0, evid = 0L, cmt = NA_character_)
  bind_rows(dose_rows, obs_rows) |>
    arrange(id, time, desc(evid)) |>
    mutate(
      CONMED_LPV    = lpv,
      MIX_LARGE_RUV = 0L,
      treatment     = label
    )
}

events <- bind_rows(
  make_cohort(50, dose_mg = 100, lpv = 0L, label = "RTV100_noLPV", id_offset =   0L),
  make_cohort(50, dose_mg = 100, lpv = 1L, label = "RTV100_LPV",   id_offset =  50L),
  make_cohort(50, dose_mg = 600, lpv = 0L, label = "RTV600_noLPV", id_offset = 100L)
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

Each cohort has 50 simulated subjects (150 total) on a 7-day BID regimen.
`MIX_LARGE_RUV = 0` selects the dominant subpopulation P1 (smaller additive
residual error) for the per-subject mixture indicator; to draw the mixture
class on a per-subject basis, sample `MIX_LARGE_RUV ~ Bernoulli(0.352)` per
ID before binding the rows.

# Simulation

```{r simulate}
mod <- readModelDb("Kappelhoff_2005_ritonavir")
sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("treatment", "CONMED_LPV", "MIX_LARGE_RUV"),
  returnType = "data.frame"
)
```

For a deterministic typical-value comparison (reproducing the structural
shape without between-subject variability), zero out the random effects:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typ <- rxode2::rxSolve(
  mod_typical, events = events,
  keep = c("treatment"), returnType = "data.frame"
) |>
  group_by(treatment, time) |>
  summarise(Cc = mean(Cc), .groups = "drop")
```

# Replicate published figures

Kappelhoff Figure 1 displays the raw concentration-time data overlay
(single-time-point samples plus full PK profiles) across all 186 subjects.
Below is the corresponding panel built from the virtual cohort: the typical
trajectories per regimen are overlaid on percentile ribbons that capture
the simulated between-subject variability.

```{r figure-1, fig.cap = "Simulated ritonavir steady-state concentration-time profiles by regimen. Solid line = typical-value (zeroRe) trajectory; ribbon = 5th-95th percentile across the virtual cohort. Replicates the spread shown in Figure 1 of Kappelhoff (2005)."}
sim_summary <- sim |>
  filter(time >= 144, time <= 168) |>
  mutate(t_in_interval = time - 144) |>
  group_by(treatment, t_in_interval) |>
  summarise(
    q05 = quantile(Cc, 0.05, na.rm = TRUE),
    q50 = quantile(Cc, 0.50, na.rm = TRUE),
    q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

sim_typ_int <- sim_typ |>
  filter(time >= 144, time <= 168) |>
  mutate(t_in_interval = time - 144)

ggplot(sim_summary, aes(t_in_interval, q50, fill = treatment, colour = treatment)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), alpha = 0.20, colour = NA) +
  geom_line(data = sim_typ_int, aes(t_in_interval, Cc, colour = treatment),
            inherit.aes = FALSE, linewidth = 0.9) +
  facet_wrap(~treatment, scales = "free_y") +
  labs(
    x = "Time within steady-state dosing interval (h)",
    y = "Ritonavir plasma concentration (mg/L)",
    title = "Simulated steady-state PK by regimen",
    caption = "Cohorts: 50 virtual subjects each, MIX_LARGE_RUV = 0. Replicates the visual envelope of Kappelhoff (2005) Figure 1."
  ) +
  theme_minimal() +
  theme(legend.position = "none")
```

# PKNCA validation

```{r pknca, warning = FALSE}
sim_nca <- sim |>
  filter(time >= 144, time <= 156) |>
  mutate(time_in_interval = time - 144) |>
  select(id, time = time_in_interval, Cc, treatment) |>
  filter(!is.na(Cc))

dose_df <- events |>
  filter(evid == 1, time == 144) |>
  select(id, time, amt, treatment) |>
  mutate(time = time - 144)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | treatment + id,
  concu = "mg/L", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | treatment + id,
  doseu = "mg"
)

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

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)

nca_tbl <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "cmin", "auclast", "cav")) |>
  group_by(treatment, PPTESTCD) |>
  summarise(
    median = median(PPORRES),
    q05    = quantile(PPORRES, 0.05),
    q95    = quantile(PPORRES, 0.95),
    .groups = "drop"
  ) |>
  arrange(treatment, factor(PPTESTCD, levels = c("tmax", "cmax", "cmin", "auclast", "cav")))

knitr::kable(
  nca_tbl,
  caption = "Steady-state NCA parameters over the 12-hour dosing interval, computed by PKNCA from the simulated cohort. Cmax / Cmin in mg/L, tmax in h, AUC0-12 in mg*h/L, Cavg in mg/L. Each row = (treatment, parameter); median and 5th-95th percentile across 50 simulated subjects."
)
```

## Comparison against published quantities

Kappelhoff (2005) does not report a steady-state NCA summary table, but the
paper does report an analytic terminal half-life of 6.4 h derived from the
structural parameter estimates `t1/2 = ln(2) * (V/F) / (CL/F)`. The sanity
check below recomputes the half-life from the typical-value `cl` and `vc`
in the simulated cohort and confirms it matches the paper.

```{r halflife-check}
typ_par <- sim |>
  filter(treatment == "RTV100_noLPV", time == 144) |>
  slice_head(n = 1) |>
  mutate(t_half_h = log(2) * vc / cl) |>
  select(treatment, cl, vc, t_half_h)

knitr::kable(
  typ_par,
  caption = "Simulated typical-value clearance (L/h), volume (L), and derived terminal half-life (h) for the no-LPV regimen. Compare with Kappelhoff (2005) Results: 'The calculated value for half-life from these estimates was 6.4 h.'"
)
```

The lopinavir effect is directly observable in the AUC ratio between
`RTV100_noLPV` and `RTV100_LPV`. The published equation
`CL/F = 10.5 * 2.72^LPV` predicts an AUC ratio of 2.72 (since AUC
inversely scales with CL/F at fixed dose and fixed F). The simulation
recovers that ratio:

```{r lpv-effect}
auc_by_arm <- nca_tbl |>
  filter(PPTESTCD == "auclast") |>
  select(treatment, median_auc = median)

ratio <- auc_by_arm |>
  pivot_wider(names_from = treatment, values_from = median_auc) |>
  mutate(
    ratio_noLPV_over_LPV = RTV100_noLPV / RTV100_LPV,
    paper_predicted      = 2.72
  )

knitr::kable(
  ratio,
  caption = "Steady-state AUC0-12 (median, mg*h/L) for 100 mg BID with and without lopinavir co-administration, and the AUC ratio. Paper-predicted value 2.72 from the CL/F covariate equation."
)
```

# Assumptions and deviations

* **Interoccasion variability omitted.** Kappelhoff (2005) Table 2 reports
  IOV on the apparent bioavailability F of 59.1% (CV). Apparent CL/F and
  V/F both depend on F, so per-visit variability in F propagates to both
  the typical CL/F and V/F at each occasion. nlmixr2 does not have a
  natural single-statement representation for IOV in a simulation-only
  context; this model omits the IOV term and carries only the IIV
  structure on `lcl`, `lvc`, and `lka`. Re-introducing IOV requires
  generating per-visit random effects on F outside the model file.
* **Race / ethnicity for the simulated cohort.** The virtual cohort is
  unweighted -- it does not reproduce the published race / ethnicity
  distribution (Caucasian 64.5%, Black 10.2%, Asian 4.3%, Hispanic 5.4%,
  missing 15.6%). Race was tested as a covariate in the Kappelhoff
  screening and was not retained in the final model, so this omission
  does not affect the simulated PK; it is documented here only so future
  users who want a demographic-faithful cohort know to add it.
* **Mixture residual error gated by `MIX_LARGE_RUV`.** The published
  `$MIX` block assigns each subject to one of two residual-error classes
  (P1: smaller addSd = 0.0600 mg/L, 64.8% of subjects; P2: larger addSd =
  0.199 mg/L, 35.2%); both classes share the same 15.4% proportional
  residual. Kappelhoff Discussion notes the mechanism behind the two
  classes could not be identified. The vignette sets
  `MIX_LARGE_RUV = 0` for all subjects (dominant class). To draw the
  class per subject, sample `MIX_LARGE_RUV ~ Bernoulli(0.352)` before
  binding cohort rows. The pattern matches the
  `Allegaert_2015_paracetamol.R` OCC = 5 piecewise-RUV implementation
  (`Cc ~ add(add_sd_eff) + prop(propSd)` where `add_sd_eff` is computed
  inside `model()`).
* **Power-form encoding of the lopinavir effect.** The published
  equation `CL/F = 10.5 * 2.72^LPV` is implemented as
  `cl <- exp(lcl + etalcl) * e_lpv_cl^CONMED_LPV` with `e_lpv_cl = 2.72`,
  so the LPV = 0 typical value is `exp(lcl) = 10.5 L/h` and the LPV = 1
  typical value is `exp(lcl) * 2.72 = 28.56 L/h`.
* **Dosing scenarios.** The virtual cohort uses three representative
  regimens (100 mg BID no-LPV, 100 mg BID with LPV, 600 mg BID no-LPV).
  Kappelhoff Methods documents at least eight distinct dose levels (100
  QD, 100 BID, 133 BID, 200 BID, 300, 400, 500, 600, 750 BID); other
  doses can be simulated by changing the `dose_mg` argument to
  `make_cohort()`.
* **Indinavir and saquinavir co-administration.** These were tested
  as covariates on ritonavir CL/F in the source paper and were not
  retained ("Only the introduction of lopinavir resulted in a
  statistically significant increase in goodness-of-fit"). The model
  carries `CONMED_LPV` only; subjects on indinavir or saquinavir use
  the LPV = 0 typical-value clearance.
