---
title: "Abacavir (Zhao 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Abacavir (Zhao 2013)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Zhao_2013_abacavir"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Zhao_2013_abacavir"))$description`
* Article: <https://doi.org/10.1111/bcp.12024>

This vignette accompanies the **abacavir** popPK model from Zhao 2013
(pooled paediatric meta-analysis across three trials). A separate
nlmixr2lib entry, `Archary_2019_abacavir`, holds the same drug in a
different paediatric population (severely malnourished, South African
inpatients) and is appropriate for that population -- both models are
abacavir popPK; choose the one closer to the population you are
simulating.

# Population

Zhao 2013 reports a pooled paediatric meta-analysis of three clinical
trials in HIV type-1-infected children: **PENTA 13** (n = 14,
2.14-12.84 years), **PENTA 15** (n = 18, 0.42-2.81 years), and the
**ARROW** pharmacokinetic substudy (n = 37, 3.62-12.54 years). The
overall cohort is 69 children spanning 0.42-12.84 years and 7.6-60.9 kg
(median 5.66 years, 17.6 kg; Table 1 of the source). Children received
either a weight-normalised 16 mg/kg/day regimen or WHO weight-band
fixed dosing (300 / 450 / 600 mg per day) split as twice daily
(8 mg/kg BID or 150-300 mg BID) or once daily; both tablet and oral
solution formulations were used. Steady-state pharmacokinetic samples
were collected at T0, T1, T2, T3, T4, T6, T8, and T12 h after dosing
(plus T24 h for the once-daily regimen). The pooled dataset contains
1065 abacavir concentrations across 138 pharmacokinetic profiles.

The same information is available programmatically:
`readModelDb("Zhao_2013_abacavir")$population` after the model is loaded.

# Source trace

Per-parameter origin (also recorded as in-file comments next to each
`ini()` entry of `inst/modeldb/specificDrugs/Zhao_2013_abacavir.R`):

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka`     | log(0.913) | Zhao 2013 Table 2 (`Ka` = 0.913 /h) |
| `lcl`     | log(20.1)  | Zhao 2013 Table 2 (`CL/F_ref` = 20.1 L/h at 17.6 kg) |
| `lvc`     | log(13.0)  | Zhao 2013 Table 2 (`V1/F_ref` = 13.0 L at 17.6 kg) |
| `lvp`     | log(13.5)  | Zhao 2013 Table 2 (`V2/F` = 13.5 L) |
| `lq`      | log(2.0)   | Zhao 2013 Table 2 (`Q/F` = 2.0 L/h) |
| `e_wt_cl` | 0.802      | Zhao 2013 Table 2 (theta_1 allometric exponent on CL/F) |
| `e_wt_vc` | 0.810      | Zhao 2013 Table 2 (theta_2 allometric exponent on V1/F) |
| `etalcl`  | 0.08763    | Zhao 2013 Table 2 (combined 21.9 % BSV + 20.4 % IOV on CL/F; see Assumptions and deviations) |
| `etalvc`  | 0.20509    | Zhao 2013 Table 2 (47.7 % CV IIV on V1/F; omega^2 = log(1 + 0.477^2)) |
| `etalvp`  | 0.28555    | Zhao 2013 Table 2 (57.5 % CV IIV on V2/F; omega^2 = log(1 + 0.575^2)) |
| `etalq`   | 0.16608    | Zhao 2013 Table 2 (42.5 % CV IIV on Q/F; omega^2 = log(1 + 0.425^2)) |
| `propSd`  | 0.382      | Zhao 2013 Table 2 (38.2 % proportional residual error) |
| Reference weight | 17.6 kg | Zhao 2013 Results ("the reference weight was the median value of our population, 17.6 kg") |
| `d/dt(depot)`, `d/dt(central)`, `d/dt(peripheral1)` | n/a | Zhao 2013 Results (two-compartment model with first-order absorption and first-order elimination) |
| `Cc <- central / vc` | n/a | Standard apparent-clearance parameterisation; dose mg, volume L -> mg/L |
| `Cc ~ prop(propSd)` | n/a | Zhao 2013 Results ("Residual variability was best described by a proportional model") |

# Virtual cohort

Original observed data are not publicly available. The cohort below is
a virtual reproduction of the paper's two subpopulation simulations
(infants + toddlers, and children) used to compare against the
published steady-state Cmax and AUC0-12 values:

* **Infants and toddlers** (n = 21 in the source; age 0.42-2.81 years;
  WT range ~7.6-15.8 kg, median 11.6 kg -- taken from PENTA 15
  demographics in Table 1).
* **Children** (n = 48 in the source; age 3.58-12.84 years; pooled from
  ARROW (n = 37, mean 20.3 kg, SD 4.0, range 14-29.8 kg) and PENTA 13
  (n = 14, mean 23.9 kg, SD 13.2, range 14-60.9 kg) -- the PENTA 13
  long right tail (up to 60.9 kg) makes the children subgroup
  considerably wider than ARROW alone).

Both subgroups are simulated at the standard regimen of 8 mg/kg twice
daily until steady state.

```{r cohort}
set.seed(20130428L)

n_per_group  <- 100L    # subjects per subpopulation
ndoses       <- 14L     # 7 days of BID dosing -> steady state
tau          <- 12      # h, BID interval
sample_hours <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 10, 12)
dose_window  <- (ndoses - 1L) * tau   # 156 h: last dose time

draw_wt_mix <- function(n) {
  # Children subgroup is a mixture of ARROW (n = 37) and PENTA 13 (n = 14)
  # demographics; 37 / (37 + 14) = 72.5 % ARROW, 27.5 % PENTA 13.
  is_penta13 <- stats::runif(n) < (14 / (14 + 37))
  arrow      <- pmin(pmax(stats::rnorm(n, 20.3,  4.0), 14, 29.8))
  penta13    <- pmin(pmax(stats::rnorm(n, 23.9, 13.2), 14, 60.9))
  ifelse(is_penta13, penta13, arrow)
}

# Helper: build one cohort. WT is sampled per subject from the
# subgroup-specific weight distribution and held constant over the
# simulation. Dose is computed per subject as 8 mg/kg * WT and given
# BID for 7 days.
make_cohort <- function(n, wt_draw, cohort_label, id_offset) {
  wt <- wt_draw(n)
  ids <- id_offset + seq_len(n)
  dose_amt <- 8 * wt   # 8 mg/kg per dose

  dose_rows <- tibble::tibble(
    id   = rep(ids, each = ndoses),
    time = rep(seq(0, by = tau, length.out = ndoses), times = n),
    amt  = rep(dose_amt, each = ndoses),
    evid = 1L,
    cmt  = 1L,
    WT   = rep(wt, each = ndoses)
  )
  obs_rows <- tibble::tibble(
    id   = rep(ids, each = length(sample_hours)),
    time = rep(dose_window + sample_hours, times = n),
    amt  = 0,
    evid = 0L,
    cmt  = NA_integer_,
    WT   = rep(wt, each = length(sample_hours))
  )
  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::mutate(cohort = cohort_label) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(
  make_cohort(
    n_per_group,
    wt_draw = function(n) pmin(pmax(stats::rnorm(n, 11.5, 2.3), 7.6, 15.8)),
    cohort_label = "Infants and toddlers",
    id_offset = 0L
  ),
  make_cohort(
    n_per_group,
    wt_draw = draw_wt_mix,
    cohort_label = "Children",
    id_offset = 1000L
  )
)

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode2(readModelDb("Zhao_2013_abacavir"))

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("WT", "cohort")
) |>
  as.data.frame()
```

For deterministic typical-value lines (reproducing the paper's
geometric-mean profile without IIV / residual scatter):

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = events,
  keep   = c("WT", "cohort")
) |>
  as.data.frame()
```

# Steady-state concentration-time profile by subpopulation

The plot below shows the simulated last-dose-interval VPC (median + 5-95
percentile band) for both subpopulations at the standard 8 mg/kg BID
regimen.

```{r figure-ss, fig.width = 7, fig.height = 4}
sim_ss <- sim |>
  dplyr::filter(time >= dose_window, !is.na(Cc)) |>
  dplyr::mutate(tau_time = time - dose_window)

sim_quantiles <- sim_ss |>
  dplyr::group_by(cohort, tau_time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(sim_quantiles, aes(tau_time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(linewidth = 0.6) +
  facet_wrap(~ cohort) +
  labs(
    x = "Time after last dose (h)",
    y = "Cc (mg/L)",
    title = "Steady-state abacavir profile by subpopulation (8 mg/kg BID, day 7)",
    caption = "Median + 5-95% interval, n = 100 per cohort. Companion to the AUC0-12 / Cmax comparison below."
  )
```

# PKNCA validation

NCA on the last dosing interval (steady-state) for each subpopulation:

```{r pknca}
pkn_in <- sim |>
  dplyr::filter(time >= dose_window, !is.na(Cc)) |>
  dplyr::mutate(tau_time = time - dose_window,
                treatment = cohort) |>
  dplyr::select(id, time = tau_time, Cc, treatment)

# Defensive guarantee: time = 0 row per (id, treatment); pre-dose at
# steady state the concentration is Ctrough, not 0, but we already have
# tau_time = 0 in sample_hours so the row is present. The bind-rows
# pattern below is a no-op safety net (distinct collapses duplicates).
pkn_in <- dplyr::bind_rows(
  pkn_in,
  pkn_in |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0,
                  Cc = pkn_in |>
                    dplyr::filter(time == 0) |>
                    dplyr::pull(Cc) |>
                    stats::median(na.rm = TRUE))
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

# One dose per subject at tau_time = 0 representing the steady-state dose.
dose_pkn <- events |>
  dplyr::filter(evid == 1L, time == dose_window) |>
  dplyr::mutate(treatment = cohort, time = 0) |>
  dplyr::select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(pkn_in, Cc ~ time | treatment + id)
dose_obj <- PKNCA::PKNCAdose(dose_pkn, amt ~ time | treatment + id,
                             route = "extravascular")

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

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

## Comparison against published values

Zhao 2013 reports the model-predicted steady-state geometric-mean Cmax
and AUC0-12 for the standard 8 mg/kg BID regimen for each
subpopulation:

* **Infants and toddlers**: Cmax = 2.5 mg/L, AUC0-12 = 6.1 mg\*h/L
* **Children**: Cmax = 3.6 mg/L, AUC0-12 = 8.7 mg\*h/L

The simulated values are summarised as the geometric mean across the
virtual cohort (n = 100 per subpopulation) to match the paper's
geometric-mean summary statistic.

```{r nca_compare}
# Compute geometric mean of simulated Cmax / AUClast per cohort, in the
# same row shape ncaComparisonTable() expects (treatment + PPTESTCD +
# PPORRES; one row per subject-parameter pair).
sim_long <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "auclast")) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(
    PPORRES = exp(mean(log(PPORRES), na.rm = TRUE)),
    .groups = "drop"
  )

reference <- tibble::tribble(
  ~treatment,             ~cmax, ~auclast,
  "Infants and toddlers",  2.5,   6.1,
  "Children",              3.6,   8.7
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = sim_long,
  reference     = reference,
  by            = "treatment",
  units         = c(cmax = "mg/L", auclast = "mg*h/L"),
  tolerance_pct = 20
)
knitr::kable(
  cmp,
  caption = paste(
    "Simulated steady-state geometric-mean Cmax and AUC0-12",
    "(8 mg/kg BID, n = 100 per cohort) vs published values in",
    "Zhao 2013 Results section. * differs from reference by >20%."
  )
)
```

# Assumptions and deviations

- **BSV + IOV on CL/F collapsed onto a single etalcl.** Zhao 2013
  Table 2 reports two separate variance components on apparent
  clearance: between-subject variability (21.9 % CV) and
  inter-occasion variability (20.4 % CV), coupled additively in
  NONMEM (`CL_ij = TVCL * exp(eta_i + kappa_ij)`). nlmixr2's library
  models typically simulate a single-occasion forward profile, so the
  two log-scale variances are added (omega^2 = 0.04686 + 0.04077 =
  0.08763) and applied through `etalcl`. This matches the Birgersson
  2016 artemisinin convention. Consumers who need the verbatim
  BSV-only profile can override `etalcl` to log(1 + 0.219^2) = 0.04686.
- **AGE, gender, and formulation tested but not retained.** Zhao 2013
  Results section reports that age, body weight, and formulation each
  produced a significant OFV drop in forward selection on CL/F, but
  only the effect of body weight on CL/F and V1/F survived backward
  elimination. These three covariates are recorded under
  `covariatesDataExcluded` in the model file so the provenance of the
  paper's covariate screen is preserved without producing convention
  warnings.
- **Allometric exponents estimated, not fixed at 0.75 / 1.** Zhao 2013
  estimated both allometric exponents from the paediatric cohort
  (0.802 for CL/F, RSE 11.6 %; 0.810 for V1/F, RSE 23.3 %). The model
  file uses the estimated point estimates without wrapping them in
  `fixed()`, matching the paper's reported uncertainty.
- **No bioavailability anchor.** All structural parameters in Zhao
  2013 are apparent (CL/F, V1/F, V2/F, Q/F); bioavailability is folded
  into them and is not separately identifiable. The model omits
  `f(depot)` rather than fixing it at 1.
- **Steady-state PKNCA validation simulates seven days of BID
  dosing.** The terminal half-life at the 17.6 kg reference
  (alpha 1.715 /h, beta 0.134 /h -> t_half_beta = 5.2 h) is well
  below the 12 h dosing interval, so seven days (14 doses) gives a
  comfortable margin to steady state before sampling.
- **No race / ethnicity covariate.** The pooled dataset spans European
  paediatric cohorts (PENTA 13 / PENTA 15) and a Ugandan paediatric
  cohort (ARROW), but the paper's Limitations section notes that
  "Information on ethnicity and other potential demographic factors
  was not available" and ethnicity is not entered into the model.
- **Cohort weight distributions are virtual.** The infants-and-toddlers
  cohort uses mean 11.5 kg, SD 2.3 kg, clipped to 7.6-15.8 kg (PENTA
  15 demographics from Table 1). The children cohort is a mixture
  weighted by trial size: 72.5 % from ARROW (mean 20.3 kg, SD 4.0,
  clipped to 14-29.8 kg) and 27.5 % from PENTA 13 (mean 23.9 kg,
  SD 13.2, clipped to 14-60.9 kg). The PENTA 13 right tail
  (up to 60.9 kg) is important to reproduce the published children
  Cmax / AUC0-12 because abacavir AUC scales as ~WT^0.198 at fixed
  8 mg/kg dosing (CL/F exponent 0.802; AUC = Dose/CL proportional to
  WT^(1 - 0.802) = WT^0.198), so heavier children contribute
  disproportionately to the geometric mean. The published Cmax /
  AUC0-12 values in the
  Zhao 2013 Results section are the model-predicted geometric means
  using the paper's own virtual-cohort simulations; the reproduction
  here uses an independent virtual cohort drawn from the same
  demographic envelope.
