---
title: "Maribavir (Sun 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Maribavir (Sun 2023)}
  %\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 readModelDb("Sun_2023_maribavir")()$reference`
* Description: `r readModelDb("Sun_2023_maribavir")()$description`
* Article: <https://doi.org/10.1002/psp4.12943>
* Supplement (Tables S1-S3, Figure S1, and the NONMEM control streams for both
  candidate models): <https://doi.org/10.1002/psp4.12943> Supporting Information

Maribavir is an oral benzimidazole riboside that inhibits the cytomegalovirus
(CMV) UL97 protein kinase. It is approved for post-transplant CMV
infection/disease refractory (with or without resistance) to valganciclovir,
ganciclovir, cidofovir or foscarnet, at 400 mg twice daily (b.i.d.). Because no
paediatric clinical data existed and adolescent recruitment was not feasible,
Sun et al. used the adult population PK (PopPK) model to simulate exposures in
virtual adolescents and support extending the indication to patients aged
12 years and older weighing at least 35 kg.

# Population

The adult PopPK model was developed in NONMEM 7.4.3 / 7.5.0 from 12 clinical
studies (Table S1): nine phase I studies, two phase II dose-ranging studies
(SHP620-202 and SHP620-203; 400, 800 and 1200 mg b.i.d.) and one phase III study
(SHP620-303; 400 mg b.i.d.). The pooled analysis therefore mixes healthy
volunteers, phase I special populations (renal impairment in study 1263-101;
moderate hepatic impairment in study 1263-103), dedicated drug-drug-interaction
cohorts (ketoconazole in 1263-102, rifampin in 1263-110) and hematopoietic cell
transplant (HCT) or solid organ transplant (SOT) recipients with CMV infection.

Table S1 accounts for 182 participants across the nine phase I studies and 235
across the two phase II studies (417 combined); the phase III participant count
is not given. Table S3 reports model-estimated steady-state exposures for 253
phase III (400 mg b.i.d.) and 232 phase II (1200 mg b.i.d.) transplant
recipients with CMV. Neither the article nor the supplement reports numeric age,
sex, race or body-weight summaries -- the baseline body weight distribution
appears only as a bar chart in Figure S1.

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

```{r population}
readModelDb("Sun_2023_maribavir")()$population$disease_state
readModelDb("Sun_2023_maribavir")()$population$dose_range
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Sun_2023_maribavir.R`. The table below
collects them in one place for review.

Two on-disk sources are cited. **Table S2** is the published parameter table
(back-transformed estimates, % RSE, 95% CI and IIV CV%). **`s005.txt`** is the
supplement's NONMEM control stream for run 171 -- the "fixed weight effect
exponents" model -- whose `$THETA` block is annotated `;from 171.cnv`, i.e. the
final converged estimates. The control stream carries full precision and the
complete `$OMEGA BLOCK(6)` and `$SIGMA`, so it is the primary source for the
`ini()` values; every one of them was cross-checked against Table S2 (see the
"Cross-validation" section below).

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` | 1.32774 (= log 3.77 L/h) | s005.txt `THETA(1) [ln_CL]`; Table S2 `CL/F` 3.77 (95% CI 3.50-4.06) |
| `lvc` | 2.92052 (= log 18.6 L) | s005.txt `THETA(2) [ln_V2]`; Table S2 `Vc/F` 18.6 (95% CI 17.3-19.8) |
| `lq` | -0.096712 (= log 0.908 L/h) | s005.txt `THETA(3) [ln_Q]`; Table S2 `Q/F` 0.908 (95% CI 0.705-1.17) |
| `lvp` | 2.15842 (= log 8.66 L) | s005.txt `THETA(4) [ln_V3]`; Table S2 `Vp/F` 8.66 (95% CI 7.05-10.6) |
| `lka` | -1.092 (= log 0.336 1/h) | s005.txt `THETA(5) [ln_Ka]`; Table S2 `Ka` 0.336 (95% CI 0.271-0.415) |
| `ltlag` | -1.30567 (= log 0.271 h) | s005.txt `THETA(6) [ln_ALAG1]`; Table S2 `Lag time` 0.271 (95% CI 0.241-0.304) |
| `e_wt_cl` | 0.75 fixed | s005.txt `THETA(7) (0.75 FIX)`; Table S2 `CL/F~weight 0.75 fixed` |
| `e_wt_vc` | 1 fixed | s005.txt `THETA(8) (1 FIX)`; Table S2 `Vc/F~weight 1 fixed` |
| `e_wt_q` | 0.75 fixed | s005.txt `THETA(9) (0.75 FIX)`; Table S2 `Q/F~weight 0.75 fixed` |
| `e_wt_vp` | 1 fixed | s005.txt `THETA(10) (1 FIX)`; Table S2 `Vp/F~weight 1 fixed` |
| `e_dis_cmv_cl` | -0.280346 (= log 0.756) | s005.txt `THETA(17) [CL~CMV]`; Table S2 `CL/F~transplant patients with CMV` 0.756 (95% CI 0.690-0.827) |
| `e_dose_ka` | -1.9439 | s005.txt `THETA(16) [KA~dose]`; Table S2 `Ka~dose` -1.94 (95% CI -2.19 to -1.70) |
| `e_cyp3a4_inh_cl` | log(0.700292) | s005.txt `THETA(11) [CL~CYP3AINH]`; **not** in Table S2 |
| `e_cyp3a4_ind_cl` | log(2.24181) | s005.txt `THETA(12) [CL~CYP3AIND]`; **not** in Table S2 |
| IIV block (6x6) | 21 lower-triangle elements | s005.txt `$OMEGA BLOCK(6)` over ETA(1..6) = CL, Vc, Q, Vp, Ka, ALAG1; diagonals reproduce all six Table S2 `IIV CV%` entries |
| `propSd` | sqrt(0.136898) = 0.370 | s005.txt `$SIGMA` third element (EPS(3), studies 202/203/303) |
| Structural model (2-cmt, first-order absorption, lag time) | n/a | s005.txt `$SUBROUTINES ADVAN4 TRANS4` + `ALAG1`; Results "Adult PopPK model" paragraph |
| `cl` covariate equation | n/a | s005.txt `$PK` `MU_1`/`CL` lines; Table S2 footnote `CL/F = 3.77 x (WT/70)^0.75 x 0.756^(patients with CMV)` |
| `ka` covariate equation | n/a | s005.txt `$PK` `KADOSE`/`KA` lines; Table S2 footnote `Ka = 0.336 x (dose/800)^-1.94` |
| Reference subject (70 kg, no CMV, 800 mg) | n/a | Table S2 footnote: "The reference population is a 70-kg individual without CMV administered a 800 mg maribavir dose" |

## Which of the paper's two models is extracted

The authors fitted the same structural model twice, differing only in whether
the four body-weight exponents were estimated or fixed:

* **Run 174** (supplement `s002.txt`) estimated them (CL/F 0.114, Vc/F 0.407,
  Q/F 1.95, Vp/F 0.663). The paper rejects this model: the estimates for
  `CL/F~weight` and `Vp/F~weight` were "imprecise" and the model was "not
  considered as sufficiently robust for subsequent simulation of adolescent
  exposure based on body weight."
* **Run 171** (supplement `s005.txt`) fixed them at the allometric values
  (0.75 for the clearance terms, 1 for the volume terms). This is the model
  reported in Table S2 and used for every simulation and exposure estimate in
  the paper.

Only run 171 is extracted, per the "final model only" rule for
model-development papers. Run 174 is a rejected candidate, not an independent
model, so it is deliberately not packaged.

```{r model-print}
mod <- readModelDb("Sun_2023_maribavir")
mod
```

## Cross-validation of the control stream against Table S2

Because the `ini()` values are taken from the control stream, it is worth
demonstrating explicitly that they reproduce the published table. All six
structural parameters back-transform to the Table S2 estimates, and all six IIV
variances reproduce the Table S2 `IIV CV%` column under that table's footnote-c
rule.

```{r cross-validate}
theta <- c(lcl = 1.32774, lvc = 2.92052, lq = -0.096712,
           lvp = 2.15842, lka = -1.092, ltlag = -1.30567)
table_s2_estimate <- c(3.77, 18.6, 0.908, 8.66, 0.336, 0.271)

# NONMEM $OMEGA BLOCK(6) lower triangle, row by row (s005.txt)
omega_lt <- c(
  0.243536,
  0.133832, 0.11584,
  -0.134868, -0.051648, 0.600237,
  0.0432159, 0.0911984, 0.429564, 0.722801,
  0.098576, 0.133826, -0.439322, -0.314392, 1.19927,
  0.00377063, -0.0215698, 0.150093, 0.114468, -0.32376, 0.177434
)
omega <- matrix(0, 6, 6)
k <- 1
for (i in 1:6) for (j in 1:i) { omega[i, j] <- omega_lt[k]; omega[j, i] <- omega_lt[k]; k <- k + 1 }

# Table S2 footnote c: CV = sqrt(omega^2)*100 if omega^2 <= 0.15, else sqrt(exp(omega^2)-1)*100
iiv_cv <- vapply(diag(omega), function(v) {
  if (v <= 0.15) sqrt(v) * 100 else sqrt(exp(v) - 1) * 100
}, numeric(1))

tibble(
  Parameter          = c("CL/F (L/h)", "Vc/F (L)", "Q/F (L/h)", "Vp/F (L)",
                         "Ka (1/h)", "Lag time (h)"),
  `Back-transformed` = round(exp(theta), 4),
  `Table S2`         = table_s2_estimate,
  `IIV CV% computed` = round(iiv_cv, 1),
  `Table S2 IIV CV%` = c(52.5, 34.0, 90.7, 103, 152, 44.1)
) |>
  knitr::kable(caption = "Control-stream run 171 values back-transformed against Table S2.")
```

The CMV effect and the dose effect on Ka also reproduce exactly, and the IIV
block is positive definite as published (no nudging was needed):

```{r cross-validate-2}
c(`exp(THETA(17)) vs Table S2 0.756` = round(exp(-0.280346), 4),
  `THETA(16) vs Table S2 -1.94`      = -1.9439)

c(`OMEGA min eigenvalue` = round(min(eigen(omega, only.values = TRUE)$values), 6),
  `positive definite`    = all(eigen(omega, only.values = TRUE)$values > 0))
```

# Virtual cohort

Original observed data are not publicly available, and the adult body-weight
distribution is published only as a bar chart (Figure S1). The adult cohorts
below therefore fix body weight at the model's own 70 kg allometric reference,
which is a principled choice rather than one fitted to the validation target
(see "Assumptions and deviations"). The adolescent cohorts follow the paper's
Methods exactly.

Cohort sizes are capped at 200 participants per arm.

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

n_per_arm <- 200L
tau <- 12 # dosing interval (h), b.i.d.

# Observation grids. The paper sampled adolescents at 0, 1, 2, 3, 4, 6, 8 and
# 12 h (Methods, "Adolescent PopPK model simulations"), so the adolescent
# replication uses that grid. Adult steady-state exposures in Table S3 were
# derived as model "individual estimates" rather than from a sparse grid, so
# the adult cohorts use a dense grid for an accurate AUC / Cmax / Cmin.
grid_adolescent <- c(0, 1, 2, 3, 4, 6, 8, 12)
grid_adult <- sort(unique(c(seq(0, tau, by = 0.05), grid_adolescent)))

# Build one steady-state b.i.d. cohort. Dose lands in `depot`; observations are
# taken on `central` (the ODE state), never on the algebraic observable `Cc`.
make_cohort <- function(n, dose, wt, label, obs_times, id_offset = 0L,
                        dis_cmv = 1) {
  subj <- tibble(
    id      = id_offset + seq_len(n),
    WT      = wt,
    DOSE    = dose,
    DIS_CMV = dis_cmv,
    CONMED_CYP3A4_INH = 0,
    CONMED_CYP3A4_IND = 0,
    arm     = label
  )
  doses <- subj |>
    mutate(time = 0, amt = dose, evid = 1L, cmt = "depot", ii = tau, ss = 1L)
  obs <- subj |>
    tidyr::crossing(time = obs_times) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "central", ii = 0, ss = 0L)
  bind_rows(doses, obs) |>
    arrange(id, time, desc(evid)) |>
    # Canonical event-table columns must come FIRST so rxode2's column
    # heuristic classifies `DOSE` as a model covariate taken from the event
    # table rather than as a missing parameter. Without this reordering
    # rxSolve() fails with "The following parameter(s) are required for
    # solving: DOSE". Same pattern as Zheng_2016_sifalimumab.Rmd.
    select(id, time, amt, evid, cmt, ii, ss, everything())
}

adult_events <- bind_rows(
  make_cohort(n_per_arm, dose = 400,  wt = 70, label = "400 mg b.i.d.",
              obs_times = grid_adult, id_offset = 0L),
  make_cohort(n_per_arm, dose = 1200, wt = 70, label = "1200 mg b.i.d.",
              obs_times = grid_adult, id_offset = 1000L)
)
stopifnot(!anyDuplicated(unique(adult_events[, c("id", "time", "evid")])))
```

The adolescent cohort follows the paper's Methods: virtual weights sampled from
a uniform distribution within each 2.5 kg increment from 25 to <30 kg and each
5 kg increment from 30 to <100 kg, dosed at 400 mg b.i.d. to steady state, with
interindividual variability included and **no** residual variability.

```{r cohort-adolescent}
band_edges <- c(seq(25, 30, by = 2.5), seq(35, 100, by = 5))
bands <- tibble(
  lower = head(band_edges, -1),
  upper = tail(band_edges, -1)
) |>
  mutate(band = sprintf("%g to <%g kg", lower, upper),
         band = factor(band, levels = band))

adolescent_events <- lapply(seq_len(nrow(bands)), function(i) {
  b <- bands[i, ]
  make_cohort(
    n_per_arm, dose = 400,
    wt = runif(n_per_arm, b$lower, b$upper),
    label = as.character(b$band),
    obs_times = grid_adolescent,
    id_offset = 10000L + (i - 1L) * 1000L
  )
}) |>
  bind_rows() |>
  mutate(band = factor(arm, levels = levels(bands$band)))
stopifnot(!anyDuplicated(unique(adolescent_events[, c("id", "time", "evid")])))

nrow(bands)
```

# Simulation

`Cc` is the individual prediction without residual error, which is what the
paper simulated ("No residual variability was included in the predictions").

```{r simulate}
# rxSolve() returns observation rows only (no dose rows, hence no `evid`
# column), so no post-filtering is needed.
sim_adult <- rxode2::rxSolve(
  mod, events = adult_events, keep = c("arm", "WT", "DOSE")
) |>
  as.data.frame()

sim_adolescent <- rxode2::rxSolve(
  mod, events = adolescent_events, keep = c("arm", "band", "WT")
) |>
  as.data.frame()

c(adult_rows = nrow(sim_adult), adolescent_rows = nrow(sim_adolescent))
```

Steady state is confirmed by the concentration at the start and end of the
dosing interval agreeing to within numerical tolerance:

```{r steady-state-check}
sim_adult |>
  filter(time %in% c(0, tau)) |>
  group_by(arm, time) |>
  summarise(median_Cc = median(Cc), .groups = "drop") |>
  tidyr::pivot_wider(names_from = time, values_from = median_Cc,
                     names_prefix = "t") |>
  mutate(`relative difference` = abs(t12 - t0) / t0) |>
  knitr::kable(digits = 5,
               caption = "Steady-state check: median Cc at t = 0 vs t = tau.")
```

# Structural checks

## Dose-dependent absorption rate

The most distinctive feature of this model is that Ka *decreases* as the
maribavir dose increases (`Ka = 0.336 * (dose/800)^-1.94`). This flattens the
concentration-time profile at higher doses and is why steady-state Cmax is
markedly less than dose-proportional even though AUC is dose-proportional.

```{r ka-dose}
tibble(DOSE = c(100, 200, 400, 800, 1200)) |>
  mutate(`Ka (1/h)` = 0.336 * (DOSE / 800)^-1.9439,
         `t1/2 of absorption (h)` = log(2) / `Ka (1/h)`) |>
  rename("Dose (mg)" = DOSE) |>
  knitr::kable(digits = 3,
               caption = "Dose effect on Ka (Table S2 footnote equation).")
```

```{r figure-profiles}
sim_adult |>
  group_by(arm, time) |>
  summarise(Q05 = quantile(Cc, 0.05), Q50 = quantile(Cc, 0.50),
            Q95 = quantile(Cc, 0.95), .groups = "drop") |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~arm) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Time within the dosing interval (h)", y = "Maribavir Cc (ug/mL)",
       title = "Simulated steady-state profiles, 70 kg transplant recipient with CMV",
       caption = paste("Median with 5th-95th percentile band, n = 200 per arm.",
                       "The 1200 mg profile is visibly flatter because Ka is",
                       "8.5-fold lower at 1200 mg than at 400 mg."))
```

## Allometry and the CMV effect

```{r allometry}
tibble(WT = c(25, 35, 50, 70, 90)) |>
  mutate(`CL/F, healthy (L/h)`  = 3.77 * (WT / 70)^0.75,
         `CL/F, CMV (L/h)`      = 3.77 * (WT / 70)^0.75 * 0.756,
         `Vc/F (L)`             = 18.6 * (WT / 70)^1) |>
  rename("Body weight (kg)" = WT) |>
  knitr::kable(digits = 3,
               caption = "Allometric scaling and the 0.756x CMV effect on CL/F.")
```

# PKNCA validation

Steady-state NCA is computed over the final dosing interval with PKNCA
(`cmax`, `tmax`, `cmin`, `auclast`, `cav`), stratified by treatment arm.

```{r pknca}
sim_nca <- sim_adult |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, arm)

# Guarantee a time = 0 row per (id, arm); pre-dose Cc for an extravascular
# steady-state interval is the trough carried over from the previous dose,
# which the simulation already provides, so existing rows win.
sim_nca <- bind_rows(
  sim_nca,
  sim_nca |> distinct(id, arm) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, arm, time, .keep_all = TRUE) |>
  arrange(id, arm, time)

dose_df <- adult_events |>
  filter(evid == 1) |>
  select(id, time, amt, arm)

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

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

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

## Comparison against published NCA (Table S3)

Table S3 reports geometric means of steady-state AUC0-tau, Cmax,ss and Cmin,ss
for adult transplant recipients with CMV at 400 mg b.i.d. (phase III, n = 253)
and 1200 mg b.i.d. (phase II, n = 232). Those are the reference values below.

```{r nca_compare}
published <- tibble::tribble(
  ~arm,              ~auclast, ~cmax, ~cmin,
  "400 mg b.i.d.",   129,      16.9,  5.13,
  "1200 mg b.i.d.",  379,      36.7,  23.7
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "arm",
  units         = c(auclast = "ug*h/mL", cmax = "ug/mL", cmin = "ug/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Simulated (70 kg reference subject) vs Table S3 geometric",
                  "means. * differs from reference by >20%."),
  align = c("l", "l", "r", "r", "r")
)
attr(cmp, "footnote")
```

Five of the six comparisons fall inside the 20% tolerance; steady-state Cmax at
400 mg b.i.d. is flagged marginally above it (roughly +21%). Critically, **all
six deviations are in the same direction** -- every simulated value is *above*
its reference -- which points to a single systematic cause rather than a
transcription error in any individual parameter.

That cause is body weight. The simulated cohort is fixed at the model's 70 kg
allometric reference, whereas the adult transplant populations summarised in
Table S3 were almost certainly heavier than 70 kg on average; the distribution
is published only as the Figure S1 bar chart, so it cannot be reproduced. Because
CL/F scales as WT^0.75, a heavier cohort has higher clearance and therefore lower
exposure, which is exactly the offset seen. The sensitivity table below shows
that a cohort median weight near 78 kg reproduces all six Table S3 values
closely.

This is offered as an *explanation* of the residual bias, **not** as a fitted
adjustment: no parameter was tuned, and the packaged model is unchanged. The
independent evidence that the transcription itself is correct is the
cross-validation section above, where all six structural parameters, all six IIV
CV% values, the CMV effect and the dose-on-Ka exponent reproduce Table S2
exactly.

```{r weight-sensitivity}
wt_sens <- lapply(c(70, 75, 78, 80, 85), function(w) {
  ev <- bind_rows(
    make_cohort(1L, dose = 400,  wt = w, label = "400 mg b.i.d.",
                obs_times = grid_adult, id_offset = 0L),
    make_cohort(1L, dose = 1200, wt = w, label = "1200 mg b.i.d.",
                obs_times = grid_adult, id_offset = 10L)
  )
  s <- rxode2::rxSolve(mod, ev, omega = NA, keep = "arm") |>
    as.data.frame() |>
    filter(!is.na(Cc))
  s |>
    group_by(arm) |>
    summarise(
      Cmax = max(Cc),
      Cmin = min(Cc),
      AUCtau = PKNCA::pk.calc.auc(conc = Cc, time = time,
                                  interval = c(0, tau),
                                  method = "lin up/log down"),
      .groups = "drop"
    ) |>
    mutate(WT = w)
}) |>
  bind_rows()

wt_sens |>
  select(WT, arm, AUCtau, Cmax, Cmin) |>
  arrange(arm, WT) |>
  rename("Body weight (kg)" = WT, "Arm" = arm,
         "AUC0-tau (ug*h/mL)" = AUCtau, "Cmax,ss (ug/mL)" = Cmax,
         "Cmin,ss (ug/mL)" = Cmin) |>
  knitr::kable(digits = 2,
               caption = paste("Typical-value steady-state exposure vs body",
                               "weight. Table S3 geometric means are 129 /",
                               "16.9 / 5.13 (400 mg) and 379 / 36.7 / 23.7",
                               "(1200 mg)."))
```

# Replicate Figure 1

Figure 1 of Sun 2023 shows simulated steady-state exposure in virtual
adolescents at 400 mg b.i.d. by body-weight band, with adult benchmarks from
Table S3 overlaid as horizontal reference lines. Panel (a) is AUC0-tau, (b) is
Cmax,ss and (c) is Cmin,ss.

```{r figure-1-metrics}
# Per-subject steady-state metrics on the paper's sparse sampling grid, using
# the same linear-up/log-down trapezoidal rule the paper specifies.
adolescent_nca <- sim_adolescent |>
  filter(!is.na(Cc)) |>
  group_by(id, band, WT) |>
  summarise(
    AUCtau = PKNCA::pk.calc.auc(conc = Cc, time = time,
                                interval = c(0, tau),
                                method = "lin up/log down"),
    Cmax = max(Cc),
    Cmin = min(Cc),
    .groups = "drop"
  )

# Adult benchmarks, all from Table S3
adult_bm <- tibble::tribble(
  ~metric,  ~geomean_400, ~pct5_400, ~pct25_400, ~geomean_1200, ~pct75_1200, ~pct95_1200,
  "AUCtau", 129,          57.0,      93.3,       379,           515,         776,
  "Cmax",   16.9,         9.16,      13.0,       36.7,          48.6,        69.9,
  "Cmin",   5.13,         1.18,      3.17,       23.7,          35.0,        55.4
)

metric_labels <- c(AUCtau = "(a) AUC0-tau (ug*h/mL)",
                   Cmax   = "(b) Cmax,ss (ug/mL)",
                   Cmin   = "(c) Cmin,ss (ug/mL)")
```

```{r figure-1, fig.height = 9, fig.width = 7.5}
long <- adolescent_nca |>
  tidyr::pivot_longer(c(AUCtau, Cmax, Cmin), names_to = "metric",
                      values_to = "value") |>
  mutate(panel = metric_labels[metric])

ref_lines <- adult_bm |>
  tidyr::pivot_longer(-metric, names_to = "benchmark", values_to = "value") |>
  mutate(
    panel = metric_labels[metric],
    colour = case_when(
      grepl("^geomean", benchmark) ~ "Adult geometric mean (400 / 1200 mg)",
      grepl("^pct5_|^pct95", benchmark) ~ "Adult 5th pct (400) / 95th pct (1200)",
      TRUE ~ "Adult 25th pct (400) / 75th pct (1200)"
    )
  )

ggplot(long, aes(band, value)) +
  geom_hline(data = ref_lines, aes(yintercept = value, colour = colour),
             linewidth = 0.4) +
  geom_boxplot(outlier.shape = NA, fill = "grey92", linewidth = 0.3,
               coef = 0) +
  stat_summary(fun.min = function(x) quantile(x, 0.05),
               fun.max = function(x) quantile(x, 0.95),
               geom = "errorbar", width = 0.35, linewidth = 0.3) +
  facet_wrap(~panel, ncol = 1, scales = "free_y") +
  scale_y_log10() +
  scale_colour_manual(values = c(
    "Adult geometric mean (400 / 1200 mg)" = "red",
    "Adult 25th pct (400) / 75th pct (1200)" = "blue",
    "Adult 5th pct (400) / 95th pct (1200)" = "orange"
  ), name = NULL) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1),
        legend.position = "bottom", legend.direction = "vertical") +
  labs(x = "Weight category", y = NULL,
       title = "Figure 1 - maribavir 400 mg b.i.d. exposure in virtual adolescents",
       caption = paste("Replicates Figure 1 of Sun 2023. Boxes are median and",
                       "interquartile range; whiskers are the 5th and 95th",
                       "percentiles. n = 200 per weight band (the paper used",
                       "1000). Reference lines are adult benchmarks from Table S3."))
```

As in the published figure, exposure rises monotonically as body weight falls,
and every metric in every band sits between the adult 400 mg and 1200 mg
geometric means for weights at or above roughly 35-40 kg.

# Replicate Table 1

Table 1 of Sun 2023 reports, for the four weight bands below 40 kg, the
percentage of virtual adolescents whose exposure clears each adult benchmark.
All eight rows are reproduced below.

```{r table-1}
bm <- function(metric, col) adult_bm[[col]][adult_bm$metric == metric]

pct <- adolescent_nca |>
  filter(band %in% c("25 to <27.5 kg", "27.5 to <30 kg",
                     "30 to <35 kg", "35 to <40 kg")) |>
  group_by(band) |>
  summarise(
    `AUC0-tau > adult geomean AUC, 400 mg`   = 100 * mean(AUCtau > bm("AUCtau", "geomean_400")),
    `Cmin,ss > adult geomean Cmin, 400 mg`   = 100 * mean(Cmin   > bm("Cmin",   "geomean_400")),
    `AUC0-tau <= adult geomean AUC, 1200 mg` = 100 * mean(AUCtau <= bm("AUCtau", "geomean_1200")),
    `AUC0-tau <= adult 75th pct AUC, 1200 mg`  = 100 * mean(AUCtau <= bm("AUCtau", "pct75_1200")),
    `AUC0-tau <= adult 95th pct AUC, 1200 mg`  = 100 * mean(AUCtau <= bm("AUCtau", "pct95_1200")),
    `Cmax,ss <= adult geomean Cmax, 1200 mg`   = 100 * mean(Cmax   <= bm("Cmax",   "geomean_1200")),
    `Cmax,ss <= adult 75th pct Cmax, 1200 mg`  = 100 * mean(Cmax   <= bm("Cmax",   "pct75_1200")),
    `Cmax,ss <= adult 95th pct Cmax, 1200 mg`  = 100 * mean(Cmax   <= bm("Cmax",   "pct95_1200")),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(-band, names_to = "Comparison", values_to = "simulated") |>
  mutate(simulated = round(simulated))

published_t1 <- tibble::tribble(
  ~Comparison,                                  ~`25 to <27.5 kg`, ~`27.5 to <30 kg`, ~`30 to <35 kg`, ~`35 to <40 kg`,
  "AUC0-tau > adult geomean AUC, 400 mg",        94, 92, 92, 85,
  "Cmin,ss > adult geomean Cmin, 400 mg",        74, 74, 74, 71,
  "AUC0-tau <= adult geomean AUC, 1200 mg",      69, 74, 80, 83,
  "AUC0-tau <= adult 75th pct AUC, 1200 mg",     87, 91, 93, 94,
  "AUC0-tau <= adult 95th pct AUC, 1200 mg",     97, 98, 99, 99,
  "Cmax,ss <= adult geomean Cmax, 1200 mg",      34, 40, 48, 60,
  "Cmax,ss <= adult 75th pct Cmax, 1200 mg",     60, 64, 75, 82,
  "Cmax,ss <= adult 95th pct Cmax, 1200 mg",     87, 91, 95, 96
) |>
  tidyr::pivot_longer(-Comparison, names_to = "band", values_to = "published")

pct |>
  left_join(published_t1, by = c("Comparison", "band")) |>
  mutate(diff = simulated - published) |>
  tidyr::pivot_wider(names_from = band,
                     values_from = c(simulated, published, diff),
                     names_glue = "{band}: {.value}") |>
  knitr::kable(caption = paste("Replicates Table 1 of Sun 2023: percentage of",
                               "virtual adolescents under 40 kg clearing each",
                               "adult exposure benchmark. 'diff' is simulated",
                               "minus published (percentage points)."))
```

```{r table-1-agreement}
agreement <- pct |>
  left_join(published_t1, by = c("Comparison", "band")) |>
  mutate(abs_diff = abs(simulated - published))

c(n_comparisons        = nrow(agreement),
  median_abs_diff_pp   = median(agreement$abs_diff),
  max_abs_diff_pp      = max(agreement$abs_diff),
  within_5_pp          = sum(agreement$abs_diff <= 5),
  within_10_pp         = sum(agreement$abs_diff <= 10))
```

Across all 32 published percentages the reproduction is close. The residual
differences are dominated by (i) Monte Carlo noise -- 200 subjects per band here
versus 1000 in the paper, giving a standard error of roughly 3 percentage points
for a percentage near 90% -- and (ii) the sparse eight-point sampling grid, which
slightly under-reads Cmax and therefore biases the Cmax rows.

# Assumptions and deviations

- **Adult body-weight distribution is not published.** Figure S1 shows the
  baseline body weight distribution as a bar chart only, and no numeric weight
  summary appears anywhere in the article or supplement. The adult validation
  cohorts therefore fix WT at 70 kg -- the model's own allometric reference
  weight -- rather than at a value chosen to match Table S3. The consequence is
  a uniform upward bias of roughly 10-21% against the Table S3 geometric means:
  five of the six comparisons sit inside the 20% tolerance and steady-state Cmax
  at 400 mg b.i.d. is flagged just outside it. Every deviation has the same sign,
  which is the signature of a size effect rather than a parameter error. The
  weight-sensitivity table quantifies this and shows that a cohort median near
  78 kg would reconcile the two; that value is *not* used in the packaged model
  and is offered only as an explanation. No parameter was tuned.
- **Age, sex and race distributions are not reported** and so are not simulated.
  Both sex effects in the final control stream are fixed to zero, so sex would
  not change the predictions in any case.
- **Only the fixed-weight-exponent model (run 171) is packaged.** The
  estimated-exponent model (run 174) is a candidate the authors explicitly
  rejected as insufficiently robust; per the "final model only" rule it is not
  extracted. Its parameter values remain available in supplement `s002.txt`.
- **The two CYP3A4 perpetrator coefficients are not in Table S2.** They appear
  only in the NONMEM control stream (`THETA(11)` = 0.700292 for an inhibitor,
  `THETA(12)` = 2.24181 for an inducer), where they are estimated. They are
  included in the model file for full fidelity, but no % RSE or confidence
  interval is available for them, and the published simulations hardcode both
  indicators to zero. Users who set `CONMED_CYP3A4_INH` or
  `CONMED_CYP3A4_IND` to 1 are extrapolating beyond what Table S2 documents.
- **Residual error: the study-specific proportional term is used.** The source
  `$ERROR` block applies `EPS(1)` (variance 0.0672821, SD 0.259) to the phase I
  studies and `EPS(3)` (variance 0.136898, SD 0.370) to studies 202/203/303,
  sharing an additive term `EPS(2)` that is fixed at exactly zero. This model
  file describes the transplant-recipient-with-CMV population, so it carries
  `EPS(3)`; the paper's own simulation control stream likewise hardcodes
  `STUDY = 202` and therefore selects `EPS(3)`. The additive term is omitted
  rather than written as `add(fixed(0))` because an additive SD of exactly zero
  contributes nothing. nlmixr2's error-model syntax cannot switch a residual SD
  on a covariate within a single endpoint, so the phase I SD is documented here
  rather than encoded.
- **Covariates screened but not retained.** Female sex (on CL/F and Vc/F) and
  Child-Pugh class B hepatic impairment (on Vc/F) are carried in the final
  control stream with coefficients `(0 FIX)`. They are recorded in the model
  file's `covariatesDataExcluded` metadata for provenance and are deliberately
  not referenced in `model()`, since a coefficient fixed at zero has no effect.
- **The adolescent replication uses 200 subjects per weight band, not 1000.**
  This respects the vignette cohort cap; the cost is roughly 3 percentage points
  of Monte Carlo noise on the Table 1 percentages.
- **No residual variability in the simulations**, matching the paper's Methods.
  All figures and tables use `Cc` (the individual prediction), not `sim`.
- **`ini()` values are taken from the run 171 control stream rather than from
  Table S2.** The control stream reports full precision and is the only on-disk
  source for the complete `$OMEGA BLOCK(6)` and `$SIGMA`. Every value was
  cross-checked against Table S2; the agreement is shown explicitly in the
  "Cross-validation" section above. No parameter value came from outside the
  article and its supplement.
