Dolutegravir (Waalewijn 2024) • nlmixr2lib

Model and source

Citation: Waalewijn H, Wasmann RE, Bamford A, Gibb DM, McIlleron HM, Colbers A, Burger DM, Denti P, and the CHAPAS-4 trial team. Population Pharmacokinetics of Dolutegravir in African Children: Results From the CHAPAS-4 Trial. J Pediatric Infect Dis Soc. 2024;13(10):533-536. doi:10.1093/jpids/piae076.
Article: https://doi.org/10.1093/jpids/piae076
Supplement (Tables S1-S5, Figure S1 and the complete NONMEM control stream): https://www.ebi.ac.uk/europepmc/webservices/rest/PMC11424990/supplementaryFiles

CHAPAS-4 (ISRCTN22964075) randomised African children living with HIV to second-line antiretroviral therapy. This secondary analysis characterises the population pharmacokinetics of the integrase inhibitor dolutegravir and asks whether the nucleos(t)ide backbone changes dolutegravir exposure.

The published supplement is unusually complete: Supplementary Table S3 carries the final parameter estimates with sampling-importance-resampling confidence intervals, and the supplement additionally reproduces the entire NONMEM control stream. Every structural decision below is therefore traced to the control stream, and every value to Table S3.

The control stream’s $THETA block holds INITIAL estimates, not final ones. $THETA(1) CL = 1.21 L/h and $THETA(2) V = 13.2 L both fall far outside the Table S3 sampling-importance-resampling intervals (0.645-0.797 and 5.84-7.46). Table S3 is the authority for every value in the model file; the control stream is the authority for the model’s structure. The one structural $THETA that does equal its final value is $THETA(14) = 0.804, the FTC/TAF bioavailability.

Population

Forty-two children from Uganda, Zambia and Zimbabwe contributed 358 dolutegravir concentrations (2 below the limit of quantification) between January 2019 and March 2021. Median age was 10.9 years (range 5.46-15.5), median weight 27.8 kg (range 15.9-53.0), and 47.6% were female (Supplementary Table S2). Children were randomised to an emtricitabine / tenofovir alafenamide (FTC/TAF) backbone (n = 21) or a standard-of-care (SOC) backbone – lamivudine/zidovudine (n = 12) or lamivudine/abacavir (n = 9).

Dolutegravir was dosed once daily per WHO weight bands: 25 mg as five 5 mg dispersible tablets for children 14 to <20 kg (n = 10), and one 50 mg film-coated tablet above 20 kg (n = 32). On the pharmacokinetic assessment day the dose was taken with a standardised low-fat breakfast (250 kcal, 5% fat) and observed by study staff; samples were drawn pre-dose and at 1, 2, 4, 6, 8, 12 and 24 h, with an extra 0.5 h sample for children taking TAF.

The same information is available programmatically via readModelDb("Waalewijn_2024_dolutegravir")()$population.

Source trace

Equation / parameter	Value	Source location
`lcl` (CL/F at 30 kg)	0.722 L/h	Table S3 “Clearance (L/h)”, 95% CI 0.645-0.797
`lvc` (Vc/F at 30 kg)	6.66 L	Table S3 “Central distribution volume (L)”, 5.84-7.46
`lq` (Q/F at 30 kg)	0.278 L/h	Table S3 “Intercompartmental clearance (L/h)”, 0.157-0.457
`lvp` (Vp/F at 30 kg)	2.00 L	Table S3 “Peripheral distribution volume (L)”, 1.45-2.71
`lka`	0.950 1/h	Table S3 “Absorption rate constant (h-1)”, 0.704-1.31
`lmtt`	1.08 h	Table S3 “Mean transit time (h)”, 0.837-1.33
`lnn`	10.1	Table S3 “Number of transit compartments”, 7.66-14.5
`lfdepot`	1, fixed	Table S3 “Relative bioavailability [F] = 1 fixed”; `$THETA(4) (1) FIX`
`e_wt_cl`, `e_wt_vc`	0.75, 1.0, fixed	Supplementary Methods para. 2; `$PK ALLMCL_WT`, `ALLMV_WT`
Reference weight	30 kg	Table S3 footnote b; `$PK TVWT = 30`
`e_taf_fdepot`	-0.196	Table S3 “Effect of FTC/TAF on bioavailability” -19.6% (-30.8%, -8.13%); `$THETA(14) = 0.804`
`etalcl` (BSV CL)	21.1% CV -> 0.044521	Table S3 “Clearance” BSV, 16.0-27.0; footnote c defines %CV = sqrt(omega^2) x 100
BOV on F (observed doses)	22.1% CV -> 0.048841	Table S3 “Bioavailability of observed doses”, 17.2-28.6
BOV on F (pre-dose occasion)	22.1% x 3.35 -> 0.548118	Table S3 “Scaling factor for increased variability of pre-dose bioavailability” 3.35 (2.25, 4.86); `$PK BOVBIO = ETA(13)*PDBOVBIO`
BOV on ka	98.3% CV -> 0.966289	Table S3 “Absorption rate constant”, 77.8-121
BOV on MTT	75.3% CV -> 0.567009	Table S3 “Mean transit time”, 59.7-95.0
`propSd`	0.091	Table S3 “Proportional error (%)” 9.1 (7.9-10.7)
`addSd`	0.07 mg/L	Table S3 “Additive error (mg/L)” 0.07 (0.04-0.10)
Transit absorption equation	n/a	Control stream `$PK KTR = (NN+1)/MTT`, `PIZZA = LOG(BIOPDKTR + 1E-12) - GAMLN(NN+1)`; `$DES TRANSIT = EXP(PIZZA + NN*LOG(KTT) - KTT)`
`d/dt(depot)`, `d/dt(central)`, `d/dt(peripheral1)`	n/a	Control stream `$DES DADT(1..3)`; `$MODEL COMP=(ABS DEFDOSE)`, `(CENTRAL DEFOBSERVATION)`, `(PERI1)`
`f(depot) <- 0`	n/a	Control stream `$PK F1 = 0` (“I need to set bioavailability in compartment 1 to 0 for this implementation of the transit compartment absorption”)
Two-compartment (not three)	n/a	`$THETA(10) V4` and `$THETA(11) Q2` are `0 FIX`; `COMP=(PERI2)` is commented out
Observation `Cc <- central / vc`	n/a	Control stream `$ERROR IPRED = A(2)/V`
Published simulated GM AUC / Ctrough	see table below	Table S5

mod <- readModelDb("Waalewijn_2024_dolutegravir")

Deterministic replication of the published simulation

Waalewijn 2024 simulated (n = 3000) steady-state AUC_0-24 and C_trough for each weight band and backbone (Table S5). Before adding variability, two identities are checked against the typical-value model.

# One typical subject per Table S5 subgroup. The paper does not report the
# within-band weight distribution it simulated, so the typical weights here are
# the geometric means of a uniform draw across each published band -- the same
# assumption used for the stochastic cohort below (see Assumptions).
tau <- 24
n_lead <- 9                      # lead-in daily doses before the assessment dose
t_dose_ss <- n_lead * tau        # 216 h: the observed assessment dose

subgroups <- tibble::tribble(
  ~subgroup,              ~band,        ~backbone,  ~WT,   ~amt, ~CONMED_TAF,
  "14-<20 kg | SOC",      "14-<20 kg",  "SOC",      16.91, 25,   0,
  "14-<20 kg | FTC/TAF",  "14-<20 kg",  "FTC/TAF",  16.91, 25,   1,
  "20-<35 kg | SOC",      "20-<35 kg",  "SOC",      27.16, 50,   0,
  "20-<35 kg | FTC/TAF",  "20-<35 kg",  "FTC/TAF",  27.16, 50,   1
)

# Build one subject's records: `n_lead + 1` daily doses onto the depot, then an
# observation grid across the final (assessment) dosing interval.
#
# `occ` sets the occasion for every dose. Table S5 reports exposure "for
# children taking dolutegravir with food", i.e. the observed-dose condition, so
# the replication below uses occasion 2 (observed-dose BOV, 22.1% CV on F)
# throughout. Occasion 1 carries the 3.35-fold inflated bioavailability BOV that
# the authors estimated for the unobserved at-home dose the pre-dose sample came
# from; it is a data-quality feature of the fit rather than a property of
# clinical dosing, and it is exercised separately below. See Assumptions.
make_records <- function(id, WT, amt, CONMED_TAF, obs_by, occ = 2) {
  doses <- tibble(
    id = id, time = seq(0, t_dose_ss, by = tau), amt = amt,
    evid = 1L, cmt = "depot"
  )
  obs <- tibble(
    id = id, time = seq(t_dose_ss, t_dose_ss + tau, by = obs_by), amt = NA_real_,
    evid = 0L, cmt = "central"
  )
  bind_rows(doses, obs) |>
    arrange(time, desc(evid)) |>
    mutate(WT = WT, CONMED_TAF = CONMED_TAF, OCC = occ)
}

ev_typ <- do.call(bind_rows, lapply(seq_len(nrow(subgroups)), function(i) {
  make_records(i, subgroups$WT[i], subgroups$amt[i], subgroups$CONMED_TAF[i],
               obs_by = 0.02) |>
    mutate(subgroup = subgroups$subgroup[i])
}))
stopifnot(!anyDuplicated(unique(ev_typ[, c("id", "time", "evid")])))

sim_typ <- rxode2::rxSolve(rxode2::zeroRe(mod), events = ev_typ,
                           keep = c("subgroup", "WT", "CONMED_TAF")) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3'
#> Warning: multi-subject simulation without without 'omega'

trap <- function(t, y) sum(diff(t) * (head(y, -1) + tail(y, -1)) / 2)

typ <- sim_typ |>
  filter(!is.na(Cc), time >= t_dose_ss) |>
  group_by(subgroup) |>
  summarise(
    WT      = unique(WT),
    taf     = unique(CONMED_TAF),
    auc     = trap(time, Cc),
    ctrough = Cc[which.min(abs(time - (t_dose_ss + tau)))],
    .groups = "drop"
  ) |>
  left_join(subgroups |> select(subgroup, amt), by = "subgroup") |>
  mutate(
    # Closed form for a linear model at steady state: AUC0-tau = F * Dose / CL.
    cl_i        = 0.722 * (WT / 30)^0.75,
    f_i         = 1 - 0.196 * taf,
    auc_closed  = f_i * amt / cl_i,
    ratio       = auc / auc_closed
  )

knitr::kable(
  typ |>
    select(subgroup, auc, auc_closed, ratio, ctrough) |>
    rename(
      "Subgroup"                = subgroup,
      "AUC0-24 solved"          = auc,
      "AUC0-24 = F*Dose/CL"     = auc_closed,
      "Ratio"                   = ratio,
      "Ctrough (mg/L)"          = ctrough
    ),
  digits = c(0, 2, 2, 5, 3),
  caption = "Typical-value steady-state exposure against the closed form."
)

Typical-value steady-state exposure against the closed form.
Subgroup	AUC0-24 solved	AUC0-24 = F*Dose/CL	Ratio	Ctrough (mg/L)
14-<20 kg \| FTC/TAF	42.79	42.79	1	0.536
14-<20 kg \| SOC	53.23	53.23	1	0.667
20-<35 kg \| FTC/TAF	59.99	59.99	1	0.853
20-<35 kg \| SOC	74.62	74.62	1	1.060

Two checks. First, the solved AUC over the assessment interval must equal F * Dose / CL exactly – both sides use the same parameters, so the only difference is trapezoidal error, and the bound can be tight. Second, the FTC/TAF arm’s AUC must be exactly 1 - 0.196 times its SOC counterpart at the same weight, which pins the covariate’s sign, magnitude and reference category.

stopifnot(
  # Internal identity: pure numerical error on a dense grid.
  all(abs(typ$ratio - 1) < 0.005),
  # The covariate effect, isolated from every other model component.
  abs(with(typ, auc[subgroup == "14-<20 kg | FTC/TAF"] /
                 auc[subgroup == "14-<20 kg | SOC"]) - 0.804) < 1e-4,
  abs(with(typ, auc[subgroup == "20-<35 kg | FTC/TAF"] /
                 auc[subgroup == "20-<35 kg | SOC"]) - 0.804) < 1e-4
)

Virtual cohort

# set.seed() seeds R's RNG (the weight draws below). It does NOT seed rxode2's
# simulation RNG, whose streams are partitioned per solver thread -- so the eta
# draws differ between a 2-core CI runner and a 16-thread workstation. Every
# assertion downstream is written to hold for any cohort this model can produce.
set.seed(20240924)

n_per_arm <- 200L

# The paper reports the weight BANDS it simulated but not the weight
# DISTRIBUTION within each band, so weights are drawn uniformly across the band.
# See Assumptions and deviations.
band_weight <- function(n, band) {
  if (band == "14-<20 kg") runif(n, 14, 20) else runif(n, 20, 35)
}

make_cohort <- function(n, sg, id_offset) {
  wt <- band_weight(n, sg$band)
  do.call(bind_rows, lapply(seq_len(n), function(k) {
    make_records(id_offset + k, wt[k], sg$amt, sg$CONMED_TAF, obs_by = 0.25) |>
      mutate(subgroup = sg$subgroup, band = sg$band, backbone = sg$backbone)
  }))
}

events <- do.call(bind_rows, lapply(seq_len(nrow(subgroups)), function(i) {
  make_cohort(n_per_arm, subgroups[i, ], id_offset = (i - 1L) * n_per_arm)
}))
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

Simulation

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("subgroup", "band", "backbone", "WT", "CONMED_TAF")
) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3
#> as a work-around try putting the mu-referenced expression on a simple line

# `Cc` is the individual prediction and carries no residual error; `sim` is the
# same prediction with the combined additive + proportional residual applied.
stopifnot(all(is.finite(sim$Cc)), all(sim$Cc >= 0))

# Replicates Figure 1 of Waalewijn 2024: simulated steady-state AUC0-24 (left)
# and C24 (right) by weight band and backbone.
exposure <- sim |>
  filter(!is.na(Cc), time >= t_dose_ss) |>
  group_by(id, subgroup, band, backbone) |>
  summarise(
    auc     = trap(time, Cc),
    ctrough = Cc[which.min(abs(time - (t_dose_ss + tau)))],
    .groups = "drop"
  )

exposure |>
  tidyr::pivot_longer(c(auc, ctrough), names_to = "metric", values_to = "value") |>
  mutate(metric = factor(metric, c("auc", "ctrough"),
                         c("AUC0-24 (mg*h/L)", "C24 (mg/L)"))) |>
  ggplot(aes(band, value, fill = backbone)) +
  geom_boxplot(outlier.shape = NA, coef = 0) +
  stat_summary(fun.data = function(x) {
    data.frame(ymin = quantile(x, 0.05), ymax = quantile(x, 0.95), y = median(x))
  }, geom = "errorbar", position = position_dodge(0.75), width = 0.25) +
  geom_hline(data = data.frame(metric = factor("AUC0-24 (mg*h/L)",
                                               c("AUC0-24 (mg*h/L)", "C24 (mg/L)")),
                               yint = 53.6),
             aes(yintercept = yint), linetype = "dashed") +
  geom_hline(data = data.frame(metric = factor(c("C24 (mg/L)", "C24 (mg/L)"),
                                               c("AUC0-24 (mg*h/L)", "C24 (mg/L)")),
                               yint = c(0.32, 0.064),
                               lt = c("dashed", "dotted")),
             aes(yintercept = yint, linetype = lt), show.legend = FALSE) +
  scale_fill_manual(values = c(SOC = "white", "FTC/TAF" = "grey65")) +
  facet_wrap(~metric, scales = "free_y") +
  labs(x = NULL, y = NULL, fill = "Backbone",
       title = "Figure 1 - simulated steady-state dolutegravir exposure",
       caption = paste("Replicates Figure 1 of Waalewijn 2024. Dashed line (AUC):",
                       "adult reference GM AUC0-24 53.6 mg*h/L. Dashed and dotted",
                       "lines (C24): EC90 0.32 mg/L and protein-adjusted IC90",
                       "0.064 mg/L.")) +
  theme_bw()

PKNCA validation

# Times are re-referenced to the assessment dose so the steady-state interval is
# [0, tau]. This is what makes the time-zero anchor row PKNCA needs for AUC (the
# t = 216 record, i.e. the pre-dose trough of the previous interval), and it
# makes `end` land exactly on an observation, which `pk.calc.ctrough()` requires
# -- it returns NA unless a concentration record sits precisely at the interval
# end. The only row filter is `!is.na(Cc)`.
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::filter(time >= t_dose_ss) |>
  dplyr::mutate(time = time - t_dose_ss) |>
  dplyr::select(id, time, Cc, subgroup) |>
  as.data.frame()

stopifnot(all(table(sim_nca$id, sim_nca$time == 0)[, "TRUE"] == 1))

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | subgroup + id,
  concu = "mg/L", timeu = "h"
)

dose_df <- events |>
  dplyr::filter(evid == 1, time == t_dose_ss) |>
  dplyr::mutate(time = 0) |>
  dplyr::select(id, time, amt, subgroup) |>
  as.data.frame()

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | subgroup + id, doseu = "mg")

# Steady-state interval: the final (observed) dosing interval.
intervals <- data.frame(
  start   = 0,
  end     = tau,
  cmax    = TRUE,
  tmax    = TRUE,
  auclast = TRUE,
  ctrough = TRUE
)

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

nca_wide <- as.data.frame(nca_res) |>
  dplyr::select(subgroup, id, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

ctrough is the concentration at the end of the dosing interval – the paper’s C_trough. cmin is deliberately not used: under a transit-absorption model the minimum over an interval can fall in the absorption phase rather than at its end.

Comparison against published NCA

Table S5 reports geometric-mean AUC_0-24 and C_trough (with CV%) from the paper’s own n = 3000 simulation. Those are the reference values below.

gm <- function(x) exp(mean(log(x)))

simulated_gm <- nca_wide |>
  group_by(subgroup) |>
  summarise(auclast = gm(auclast), ctrough = gm(ctrough), .groups = "drop")

published <- tibble::tribble(
  ~subgroup,             ~auclast, ~ctrough,
  "14-<20 kg | SOC",     54.0,     0.737,
  "14-<20 kg | FTC/TAF", 42.5,     0.572,
  "20-<35 kg | SOC",     74.4,     1.11,
  "20-<35 kg | FTC/TAF", 60.8,     0.934
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = simulated_gm,
  reference     = published,
  by            = "subgroup",
  units         = c(auclast = "mg*h/L", ctrough = "mg/L"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  digits  = 3,
  caption = paste("Simulated geometric means vs Waalewijn 2024 Supplementary",
                  "Table S5. * differs from reference by >20%.")
)

Simulated geometric means vs Waalewijn 2024 Supplementary Table S5. * differs from reference by >20%.
NCA parameter	subgroup	Reference	Simulated	% diff
AUClast (mg*h/L)	14-<20 kg \| SOC	54	53	-1.9%
AUClast (mg*h/L)	14-<20 kg \| FTC/TAF	42.5	43.9	+3.2%
AUClast (mg*h/L)	20-<35 kg \| SOC	74.4	78	+4.9%
AUClast (mg*h/L)	20-<35 kg \| FTC/TAF	60.8	60.9	+0.2%
Ctrough (mg/L)	14-<20 kg \| SOC	0.737	0.713	-3.3%
Ctrough (mg/L)	14-<20 kg \| FTC/TAF	0.572	0.633	+10.7%
Ctrough (mg/L)	20-<35 kg \| SOC	1.11	1.21	+8.8%
Ctrough (mg/L)	20-<35 kg \| FTC/TAF	0.934	0.935	+0.1%

chk <- simulated_gm |>
  rename(auc_sim = auclast, ctau_sim = ctrough) |>
  left_join(published |> rename(auc_pub = auclast, ctau_pub = ctrough), by = "subgroup") |>
  mutate(auc_pct = 100 * (auc_sim / auc_pub - 1),
         ctau_pct = 100 * (ctau_sim / ctau_pub - 1))

stopifnot(
  # AUC is the structural check: a mis-transcribed clearance, reference weight,
  # dose or bioavailability effect moves every subgroup by tens of percent, so
  # these bounds still go red on a real regression. They are set well outside
  # the run-to-run spread rather than at it: with 200 subjects per arm the
  # geometric mean carries roughly 2.5% standard error from the clearance BSV
  # and the bioavailability BOV alone, and the within-band weight draw adds
  # more. Observed across renders: median |deviation| ~2.5%, max ~5%.
  median(abs(chk$auc_pct)) < 10,
  max(abs(chk$auc_pct)) < 15,
  # Ctrough additionally depends on the absorption parameters -- whose
  # between-occasion variability is 98% CV on ka and 75% on mean transit time --
  # and on the unpublished within-band weight distribution, so it gets a looser
  # bound. Observed across renders: median ~6%, max ~11%.
  median(abs(chk$ctau_pct)) < 20,
  max(abs(chk$ctau_pct)) < 30,
  # Direction of the backbone effect, per subgroup pair. Not a near-zero
  # comparison: the model applies a fixed 19.6% reduction, far larger than the
  # Monte-Carlo spread of either arm's geometric mean.
  chk$auc_sim[chk$subgroup == "14-<20 kg | FTC/TAF"] <
    chk$auc_sim[chk$subgroup == "14-<20 kg | SOC"],
  chk$auc_sim[chk$subgroup == "20-<35 kg | FTC/TAF"] <
    chk$auc_sim[chk$subgroup == "20-<35 kg | SOC"]
)

Efficacy targets

Waalewijn 2024 reports that geometric-mean C_trough exceeded the 0.32 mg/L EC₉₀ in every subgroup, and that over 99% of individual C_trough values exceeded the 0.064 mg/L protein-adjusted IC₉₀.

targets <- nca_wide |>
  left_join(exposure |> distinct(id, subgroup), by = c("id", "subgroup")) |>
  group_by(subgroup) |>
  summarise(
    gm_ctau      = gm(ctrough),
    pct_over_ec90 = 100 * mean(ctrough > 0.32),
    pct_over_ic90 = 100 * mean(ctrough > 0.064),
    .groups = "drop"
  )

knitr::kable(
  targets |>
    rename(
      "Subgroup"                 = subgroup,
      "GM Ctrough (mg/L)"        = gm_ctau,
      "% > EC90 (0.32 mg/L)"     = pct_over_ec90,
      "% > IC90 (0.064 mg/L)"    = pct_over_ic90
    ),
  digits = c(0, 3, 1, 1),
  caption = "Simulated attainment of the dolutegravir efficacy targets."
)

Simulated attainment of the dolutegravir efficacy targets.
Subgroup	GM Ctrough (mg/L)	% > EC90 (0.32 mg/L)	% > IC90 (0.064 mg/L)
14-<20 kg \| FTC/TAF	0.633	91.5	100
14-<20 kg \| SOC	0.713	94.0	100
20-<35 kg \| FTC/TAF	0.935	99.0	100
20-<35 kg \| SOC	1.208	99.5	100


stopifnot(
  # Paper: GM Ctrough "consistently exceeded" 0.32 mg/L in all subgroups.
  all(targets$gm_ctau > 0.32),
  # Paper: "over 99% of individuals in each subgroup" above 0.064 mg/L.
  # Bounded at 97% rather than 99% so the gate survives cohort-to-cohort
  # variation while still going red on a real regression (observed: 100%).
  all(targets$pct_over_ic90 >= 97)
)

The unobserved pre-dose occasion

The paper’s distinguishing structural finding is that bioavailability varied far more on the occasion whose dose was taken at home and not observed by study staff: 22.1% CV for observed doses versus about 75% for the pre-dose occasion. Re-running the same cohort with every dose assigned to occasion 1 rather than occasion 2 must widen the exposure distribution substantially while leaving its centre alone, because both occasions share a geometric mean of 1.

events_occ1 <- do.call(bind_rows, lapply(seq_len(nrow(subgroups)), function(i) {
  sg <- subgroups[i, ]
  wt <- band_weight(n_per_arm, sg$band)
  do.call(bind_rows, lapply(seq_len(n_per_arm), function(k) {
    make_records((i - 1L) * n_per_arm + k, wt[k], sg$amt, sg$CONMED_TAF,
                 obs_by = 1, occ = 1) |>
      mutate(subgroup = sg$subgroup)
  }))
}))

occ_cmp <- bind_rows(
  rxode2::rxSolve(mod, events = events_occ1, keep = "subgroup") |>
    as.data.frame() |> mutate(occasion = "1 (unobserved home dose)"),
  sim |> mutate(occasion = "2 (observed dose)")
) |>
  filter(!is.na(Cc), time == t_dose_ss + tau) |>
  group_by(occasion) |>
  summarise(gm_c24 = exp(mean(log(Cc))), sd_log_c24 = sd(log(Cc)), .groups = "drop")

knitr::kable(
  occ_cmp |>
    rename("Occasion" = occasion, "GM C24 (mg/L)" = gm_c24,
           "SD of log C24" = sd_log_c24),
  digits = 3,
  caption = "Spread of the 24 h concentration by dosing occasion."
)

Spread of the 24 h concentration by dosing occasion.
Occasion	GM C24 (mg/L)	SD of log C24
1 (unobserved home dose)	0.785	0.922
2 (observed dose)	0.845	0.568


spread_ratio <- occ_cmp$sd_log_c24[occ_cmp$occasion == "1 (unobserved home dose)"] /
  occ_cmp$sd_log_c24[occ_cmp$occasion == "2 (observed dose)"]

stopifnot(
  # The bioavailability BOV inflates 3.35-fold on occasion 1, but ka, MTT, CL
  # variability and the weight spread are unchanged and dilute the ratio, so the
  # bound is a floor rather than a target. Observed across renders: ~1.5-1.8.
  spread_ratio > 1.25,
  # Both occasions have a geometric mean of 1 on F, so the centres must agree.
  abs(log(occ_cmp$gm_c24[1] / occ_cmp$gm_c24[2])) < 0.30
)

The paper additionally notes that children of 14 to <20 kg taking FTC/TAF had a 21% lower GM AUC than the 53.6 mg*h/L adult reference.

adult_ref <- 53.6
pct_vs_adult <- simulated_gm |>
  mutate(pct_of_adult = 100 * (auclast / adult_ref - 1))

knitr::kable(
  pct_vs_adult |>
    select(subgroup, auclast, pct_of_adult) |>
    rename(
      "Subgroup"                        = subgroup,
      "GM AUC0-24 (mg*h/L)"             = auclast,
      "% difference vs adult reference" = pct_of_adult
    ),
  digits = c(0, 1, 1),
  caption = paste("Simulated GM AUC0-24 against the 53.6 mg*h/L adult reference",
                  "(50 mg film-coated tablet with food).")
)

Simulated GM AUC0-24 against the 53.6 mg*h/L adult reference (50 mg film-coated tablet with food).
Subgroup	GM AUC0-24 (mg*h/L)	% difference vs adult reference
14-<20 kg \| FTC/TAF	43.9	-18.2
14-<20 kg \| SOC	53.0	-1.1
20-<35 kg \| FTC/TAF	60.9	13.6
20-<35 kg \| SOC	78.0	45.6


# Paper: -21% for the 14-<20 kg FTC/TAF subgroup. Bound set well outside the
# observed spread across renders (-19% to -22%).
stopifnot(
  abs(pct_vs_adult$pct_of_adult[pct_vs_adult$subgroup == "14-<20 kg | FTC/TAF"] -
        (-21)) < 8
)

Assumptions and deviations

Within-band weight distribution. Waalewijn 2024 states the weight bands it simulated (14 to <20 kg and 20 to <35 kg) but not the weight distribution inside each band, which is what sets the geometric-mean AUC because CL scales as WT^0.75. Weights here are drawn uniformly across each band. No parameter was adjusted to improve agreement with Table S5.
transit() is expanded by hand. The model file writes the Savic transit input out explicitly with podo(depot) / tad(depot) rather than calling rxode2’s transit(n, mtt, bio) macro. Under rxode2 5.1.7 that macro evaluates to zero inside an nlmixr2 model function (it works in a plain rxode2() model), so a model using it together with the required f(depot) <- 0 simulates flat-zero concentrations with no error or warning. The hand-expanded form is numerically identical to the working transit() path and is a more literal transcription of the paper’s $DES block. This affects other nlmixr2lib models that call transit(); it is a library-wide issue, not a property of this paper.
The pre-dose bioavailability inflation is folded into an occasion variance. The paper estimates a scaling factor of 3.35 (95% CI 2.25-4.86) that multiplies the bioavailability BOV eta on the unobserved pre-dose occasion. nlmixr2lib has no canonical name for a per-occasion multiplier on an eta’s standard deviation, so occasion 1’s variance is set to (0.221 * 3.35)^2 = 0.548118 while the observed occasions keep 0.221^2. The two forms are numerically identical – this is a re-parameterisation, not an approximation – but the 3.35 folds cannot be re-estimated as a single parameter in a refit.
Occasion assignment in the simulation. The model carries three occasion slots, matching the control stream’s IF(OCC==1/2/3) blocks. The simulation above uses occasion 1 for every at-home (unobserved) lead-in dose and occasion 2 for the observed assessment dose. All lead-in doses therefore share one eta draw rather than one per day; because mean transit time is 1.08 h, absorption from any earlier dose is complete long before the assessment interval, so this affects only the small amount carried into the interval. Occasion 3 is unused here and is present for CHAPAS-4’s second pharmacokinetic visit.
Reported LLOQ conflict. The main text states an LC-MS/MS lower limit of quantification of 0.05 mg/L, while the supplement’s $ERROR block sets LLOQ = 0.01 and LOD = 0.00432 mg/L. This affects only the below-limit handling (Beal’s M6 method), which is a fitting-time concern and is not part of the packaged model; the additive residual error is taken directly from Table S3 as the assembled value (0.07 mg/L), not rebuilt from $THETA(6) + 0.2 * LOD.
Screened but unretained covariates. Fat-free mass, Schwartz creatinine clearance, ALT, AST, dolutegravir formulation (dispersible vs film-coated tablet), and TAF / tenofovir AUC as continuous replacements for the categorical backbone indicator were all tested and none was retained. They are recorded in the model file’s covariatesDataExcluded list rather than being dropped silently. The null formulation effect is the paper’s most notable negative result and is specific to dosing with food.
No BSV outside clearance. Every between-subject $OMEGA other than clearance is 0 FIX in the control stream, so the packaged model carries a single BSV eta. This is the authors’ structure, not a simplification.