Skip to contents

Vancomycin (Jung 2024)

Source: vignettes/articles/Jung_2024_vancomycin.Rmd
Jung_2024_vancomycin.Rmd

Model and source

  • Citation: Jung D, Kishk OA, Bhutta AT, Cummings GE, El Sahly HM, Virk MK, Moffett BS, Morris Daniel JL, Watanabe A, Fishbane N, Kotloff KL, Gu K, Ghazaryan V, Gobburu JVS, Akcan-Arikan A, Campbell JD. Evaluation of Vancomycin Dose Needed to Achieve 24-Hour Area Under the Concentration-Time Curve to Minimum Inhibitory Concentration Ratio Greater Than or Equal to 400 Using Pharmacometric Approaches in Pediatric Intensive Care Patients. Crit Care Explor. 2024;6(10):e1159. doi:10.1097/CCE.0000000000001159
  • Description: Two-compartment IV population PK model for vancomycin in critically ill children (90 days to <18 years) in the PICU who are not on extracorporeal therapy (Jung 2024). Clearance and intercompartmental clearance scale allometrically with body weight (exponent 0.75, reference 20 kg) and clearance additionally scales as a power function of bedside-Schwartz estimated GFR (exponent 0.5259, reference 141 mL/min/1.73 m^2); central and peripheral volumes scale linearly with body weight (exponent 1, reference 20 kg). Residual variability is a log error model.
  • Article: https://doi.org/10.1097/CCE.0000000000001159
  • Supplement (eMethods, eTables 1-4, eFigures 1-10): http://links.lww.com/CCX/B406

Jung 2024 is a prospective, two-site population PK study of IV vancomycin in critically ill children. The paper has three quantitative products, and this vignette exercises all three:

  1. A two-compartment population PK model (Table 2) - the model packaged here as Jung_2024_vancomycin.
  2. A nonlinear regression predicting AUC24 from a single trough concentration (eTable 3, tabulated in Table 3). This is a regression, not an ODE model, so it is not packaged as a separate model file; it is reproduced below as a plain function and used as an independent cross-check.
  3. A target-attainment simulation of four dosing regimens against AUC24/MIC >= 400 (Table 4, with the underlying AUC24 distributions in eTable 4) - reproduced below from the packaged model.

Population

The model was fit to 301 critically ill children (302 enrolled between 2018-04-15 and 2020-02-06; one excluded for unavailable PK data) aged 90 days to under 18 years, admitted to the PICU at the University of Maryland Children’s Hospital or Texas Children’s Hospital and prescribed IV vancomycin for any indication (Jung 2024 Table 1). Median age was 6.0 years (IQR 1.6-13.0) and median body weight 20.1 kg (IQR 10.7-41.6); 58.3% were male; 72.8% White, 16.9% Black or African American, 10.3% other, and 40.7% Hispanic or Latino. Median bedside-Schwartz eGFR was 138.3 mL/min/1.73 m^2 (IQR 101.5-179.8), with 81.8% of subjects in the normal renal-function stratum and only 1.0% severe - this is a cohort whose renal function skews high, which is what makes the paper’s “adjust the initial dose upwards for supranormal eGFR” conclusion follow.

Patients on any form of extracorporeal therapy (ECMO, CRRT, extracorporeal liver support), those with cardiopulmonary bypass within 7 days of starting vancomycin, and those on chronic dialysis were excluded; a separate 33-subject extracorporeal-therapy cohort was recruited but not analysed in this paper. Dosing was chosen by the treating clinicians, giving actual regimens of 5-31 mg/kg q6h, most frequently 1-hour infusions of 15 mg/kg every 6-9 hours. A median of 4 samples per subject (IQR 3-6) gave 1027 serum concentration records, assayed by validated LC-MS/MS.

The same information is available programmatically via the model’s population metadata (readModelDb("Jung_2024_vancomycin")()$population).

Source trace

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

Equation / parameter Value Source location
lcl (CL at WT 20 kg, eGFR 141) 2.763 L/h (RSE 2.746%) Jung 2024 Table 2, clearance row
lvc (Vc at WT 20 kg) 8.63 L (RSE 14.39%) Jung 2024 Table 2, central volume row
lq (Q at WT 20 kg) 3.808 L/h (RSE 20.45%) Jung 2024 Table 2, intercompartmental clearance row
lvp (Vp at WT 20 kg) 9.566 L (RSE 14.46%) Jung 2024 Table 2, peripheral volume row
e_wt_cl_q 0.75, fixed Jung 2024 Table 2: (WT/20)^0.75 printed on both CL and Q, no RSE; eTable 2 names the reference model “2-compartment model with allometric scaling”
e_wt_vc_vp 1, fixed Jung 2024 Table 2: (WT/20) printed on both Vc and Vp with no exponent, i.e. linear
e_crcl_cl 0.5259 (RSE 14.92%) Jung 2024 Table 2: (eGFR/141)^0.5259
Reference weight 20 kg Jung 2024 Table 2 and Results (“typical values for a subject weighing 20 kg”)
Reference eGFR 141 mL/min/1.73 m^2 Jung 2024 Table 2 and Results (“CL 2.76 L/h for eGFR of 141 mL/min/1.73 m^2”)
eGFR definition 0.413 * height (cm) / Scr (mg/dL) Jung 2024 Table 1 footnote a (bedside Schwartz)
etalcl omega^2 = 0.135678 Jung 2024 Table 2: 38.12% CV (RSE 15.94%); supplement eMethods 4.1 gives %CV = sqrt(exp(omega^2) - 1) * 100, so omega^2 = log(CV^2 + 1)
etalvc omega^2 = 0.719144 Jung 2024 Table 2: 102.6% CV (RSE 28.64%); same conversion
BSV on Q, Vp absent Jung 2024 Table 2: reported as “NE” (not estimated)
expSd (lnorm) 0.3310 Jung 2024 Table 2: residual variability (log scale) variance 0.1096 (RSE 12.69%), sd 0.3310; supplement eMethods 4.2.4 specifies the log error model log(Y) = log(Yhat) + eps
Two-compartment linear-elimination ODEs n/a Jung 2024 Results and Conclusions: “A two-compartment model allometrically scaled with body weight with linear elimination”
IV infusion, no depot / bioavailability n/a Jung 2024 Methods: IV vancomycin, 1-hour infusions (Table 4 footnote)
AUC24 identity Dose * (24/tau) / CL n/a Jung 2024 Methods and supplement eMethods “Target Attainment Analysis”
AUC24-from-trough regression coefficients 366, 54.9, 0.0894, 0.495 Jung 2024 eTable 3

Two source details deserve emphasis because they are easy to get wrong and both are resolved only by the supplement, not the main text:

  • The BSV column is an exact log-normal CV, not 100 * omega. eMethods 4.1 states %CV = sqrt(exp(omega^2) - 1) * 100. Treating 102.6% as 100 * omega would have set omega^2 = 1.0527 instead of the correct 0.7191 - a 46% overstatement of the Vc variance.
  • The residual error is a log error model, not a proportional one. eMethods 4.2.4 gives log(Y) = log(Yhat) + eps, which is nlmixr2’s lnorm(), and the “Estimate” column of Table 2 is the variance (sqrt(0.1096) = 0.3311, matching the separately-printed 0.3310 sd).

The packaged model reproduces Table 2 exactly 

mod <- readModelDb("Jung_2024_vancomycin")

typical <- rxode2::rxSolve(
  mod,
  events = rxode2::et(amt = 0, cmt = "central") |> rxode2::et(0, cmt = "central"),
  params = c(WT = 20, CRCL = 141),
  omega  = NA,
  returnType = "data.frame"
)
#> ℹ parameter labels from comments will be replaced by 'label()'

table2 <- tibble::tibble(
  Parameter = c("Clearance (L/h)", "Central volume (L)",
                "Intercompartmental clearance (L/h)", "Peripheral volume (L)"),
  Published = c(2.763, 8.63, 3.808, 9.566),
  Model     = round(unlist(typical[1, c("cl", "vc", "q", "vp")]), 4)
)
stopifnot(max(abs(table2$Model - table2$Published)) < 1e-8)

knitr::kable(
  table2,
  caption = paste("Typical values at the Jung 2024 reference subject (WT 20 kg,",
                  "eGFR 141 mL/min/1.73 m^2) versus Jung 2024 Table 2."),
  align = c("l", "r", "r")
)
Typical values at the Jung 2024 reference subject (WT 20 kg, eGFR 141 mL/min/1.73 m^2) versus Jung 2024 Table 2.
Parameter Published Model
Clearance (L/h) 2.763 2.763
Central volume (L) 8.630 8.630
Intercompartmental clearance (L/h) 3.808 3.808
Peripheral volume (L) 9.566 9.566

Virtual cohort

Original subject-level data are not publicly available. Jung 2024 built its virtual population of 2250 subjects “by randomly sampling with replacement the baseline covariate vector of subjects from the study”; only the marginal median and IQR of each covariate are published (Table 1). The cohort below therefore uses log-normal marginals matched exactly to the published median and IQR of body weight and eGFR, coupled by a seeded Latin-hypercube permutation, with the two etas placed on their own exact quantiles. This is deterministic - no random draw - so the numbers below are reproducible and cannot be moved by reseeding (see Assumptions and deviations for the independence caveat).

n_per_arm <- 200L # cohort cap; 200/arm is ample and keeps the render fast
set.seed(20240908)

quantile_grid <- (seq_len(n_per_arm) - 0.5) / n_per_arm

# Log-normal quantile function pinned to a published median and IQR.
lognormal_at <- function(p, median, q25, q75) {
  sdlog <- (log(q75) - log(q25)) / (2 * stats::qnorm(0.75))
  median * exp(stats::qnorm(p) * sdlog)
}

base_cohort <- tibble::tibble(
  # Jung 2024 Table 1: WT median 20.1 kg (IQR 10.7-41.6)
  WT     = lognormal_at(quantile_grid[sample(n_per_arm)], 20.1, 10.7, 41.6),
  # Jung 2024 Table 1: eGFR median 138.3 (IQR 101.5-179.8) mL/min/1.73 m^2
  CRCL   = lognormal_at(quantile_grid[sample(n_per_arm)], 138.3, 101.5, 179.8),
  # etas on their own exact normal quantiles (omega^2 from the model file)
  etalcl = stats::qnorm(quantile_grid[sample(n_per_arm)], 0, sqrt(0.135678)),
  etalvc = stats::qnorm(quantile_grid[sample(n_per_arm)], 0, sqrt(0.719144))
)

tibble::tibble(
  Covariate = c("Body weight (kg)", "eGFR (mL/min/1.73 m^2)"),
  `Published median (IQR)` = c("20.1 (10.7-41.6)", "138.3 (101.5-179.8)"),
  `Cohort median (IQR)` = c(
    sprintf("%.1f (%.1f-%.1f)", median(base_cohort$WT),
            quantile(base_cohort$WT, 0.25), quantile(base_cohort$WT, 0.75)),
    sprintf("%.1f (%.1f-%.1f)", median(base_cohort$CRCL),
            quantile(base_cohort$CRCL, 0.25), quantile(base_cohort$CRCL, 0.75))
  )
) |>
  knitr::kable(caption = "Virtual cohort covariates versus Jung 2024 Table 1.")
Virtual cohort covariates versus Jung 2024 Table 1.
Covariate Published median (IQR) Cohort median (IQR)
Body weight (kg) 20.1 (10.7-41.6) 20.1 (10.2-39.5)
eGFR (mL/min/1.73 m^2) 138.3 (101.5-179.8) 138.3 (104.1-183.8)

Simulation of the four dosing regimens

The four regimens are the ones Jung 2024 simulated (Table 4): 15 mg/kg q6h, 15 mg/kg q8h, 20 mg/kg q6h, and 30 mg/kg q8h, each as a 1-hour infusion (Table 4 footnote). Steady state is imposed exactly with ss = 1 on the first dose rather than by simulating a burn-in, so the AUC over the 0-24 h window is the steady-state AUC24 the paper’s target-attainment analysis used. Each arm gets a disjoint block of subject IDs.

regimens <- tibble::tibble(
  treatment = c("15 mg/kg q6h", "15 mg/kg q8h", "20 mg/kg q6h", "30 mg/kg q8h"),
  mg_per_kg = c(15, 15, 20, 30),
  tau       = c(6, 8, 6, 8),
  id_offset = c(0L, 200L, 400L, 600L)
)

make_arm <- function(mg_per_kg, tau, id_offset, treatment) {
  subj <- base_cohort |>
    dplyr::mutate(
      id        = id_offset + dplyr::row_number(),
      treatment = treatment,
      dose_mg   = mg_per_kg * WT
    )

  doses <- subj |>
    tidyr::expand_grid(dose_time = seq(0, 24 - tau, by = tau)) |>
    dplyr::transmute(
      id, treatment, WT, CRCL, etalcl, etalvc,
      time = dose_time,
      amt  = dose_mg,
      rate = dose_mg / 1, # 1-hour infusion (Jung 2024 Table 4 footnote)
      evid = 1L,
      cmt  = "central",
      # ss = 1 on the time-zero dose imposes exact steady state for the
      # q-tau regimen; later doses inside the window are ordinary doses.
      ss   = ifelse(dose_time == 0, 1L, 0L),
      ii   = ifelse(dose_time == 0, tau, 0)
    )

  # Observation grid: 0.25 h throughout, refined to 0.05 h over the infusion
  # and the immediate post-infusion peak so that trapezoidal / lin-up-log-down
  # AUC reproduces the analytic integral. Records are also placed immediately
  # BEFORE each dose: an observation sharing a dose's timestamp is evaluated
  # after that dose, so without these the interval minimum would be the
  # concentration a quarter-hour early rather than the true trough.
  peak_times <- as.numeric(outer(seq(0, 24 - tau, by = tau),
                                 seq(0, 1.5, by = 0.05), "+"))
  pre_dose_times <- seq(tau, 24, by = tau) - 1e-4
  grid <- sort(unique(c(seq(0, 24, by = 0.25), pre_dose_times,
                        peak_times[peak_times <= 24])))

  obs <- subj |>
    tidyr::expand_grid(time = grid) |>
    dplyr::transmute(
      id, treatment, WT, CRCL, etalcl, etalvc, time,
      amt = NA_real_, rate = 0, evid = 0L, cmt = "central", ss = 0L, ii = 0
    )

  dplyr::bind_rows(doses, obs) |> dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(lapply(
  seq_len(nrow(regimens)),
  function(i) {
    with(regimens[i, ], make_arm(mg_per_kg, tau, id_offset, treatment))
  }
))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
stopifnot(dplyr::n_distinct(events$id) == 4L * n_per_arm)
# etalcl / etalvc are supplied as columns, so omega = NA gives a fully
# deterministic solve (no resampling of the random effects).
sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("treatment", "WT", "CRCL"),
  omega  = NA,
  returnType = "data.frame"
)
#> Warning: multi-subject simulation without without 'omega'

stopifnot(!anyNA(sim$Cc), all(sim$Cc >= 0))

Replicate published figures and findings

Steady-state concentration-time profiles

Jung 2024 does not publish a typical-value concentration-time figure (its eFigures 4-7 are diagnostics on observed data), so the panel below is illustrative rather than a figure replication. It shows the biphasic decline the paper describes, “with the second phase beginning at approximately 8 hours after the end of infusion”.

# rxSolve returns observation records only (no evid column), so no filter is
# needed here.
sim |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05), Q50 = median(Cc), Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~treatment) +
  scale_y_log10() +
  labs(
    x = "Time within the steady-state dosing interval window (h)",
    y = "Serum vancomycin (ug/mL)",
    title = "Steady-state profiles, median and 5th-95th percentile",
    caption = "Illustrative; Jung 2024 publishes no typical-value profile figure."
  )

Clearance is a less-than-proportional function of eGFR

This is the paper’s headline structural finding (Key Points: “As the eGFR increased, the clearance of vancomycin showed a less than proportional increase”). The exponent 0.5259 is what encodes it.

egfr_grid <- seq(20, 250, by = 1)
tibble::tibble(
  CRCL = egfr_grid,
  `Jung 2024 (exponent 0.5259)` = 2.763 * (egfr_grid / 141)^0.5259,
  `Proportional (exponent 1)`   = 2.763 * (egfr_grid / 141)
) |>
  tidyr::pivot_longer(-CRCL, names_to = "Scaling", values_to = "cl") |>
  ggplot(aes(CRCL, cl, linetype = Scaling)) +
  geom_line() +
  geom_vline(xintercept = 141, colour = "grey50") +
  labs(
    x = "Bedside-Schwartz eGFR (mL/min/1.73 m^2)",
    y = "Clearance at WT 20 kg (L/h)",
    title = "Less-than-proportional eGFR effect on vancomycin clearance",
    caption = "Vertical line: the 141 mL/min/1.73 m^2 reference. Jung 2024 Table 2."
  ) +
  theme(legend.position = "bottom")

PKNCA validation

AUC24 is computed with PKNCA over the steady-state 0-24 h window. The dose object declares duration = 1 so PKNCA treats the records as 1-hour infusions rather than bolus doses.

sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# ss = 1 means the time-zero record already carries the steady-state trough,
# so a time-zero concentration exists for every subject; assert it rather than
# patching one in.
stopifnot(sum(sim_nca$time == 0) == 4L * n_per_arm)

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

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

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

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

auc24 <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "auclast", start == 0, end == 24) |>
  dplyr::select(id, treatment, AUC24 = PPORRES)
stopifnot(nrow(auc24) == 4L * n_per_arm)

Gate 1 - PKNCA AUC24 equals the paper’s own analytic AUC24

Jung 2024 never integrates a simulated profile: both its target-attainment analysis and its AUC24-versus-trough regression define AUC24 = Dose / CL * (24 / tau) (Methods; supplement eMethods “Target Attainment Analysis”). For a linear model at steady state that identity is exact, so reproducing it from the ODE solution end-to-end validates the structural model, the steady-state setup, and the NCA all at once.

analytic <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::group_by(id, treatment) |>
  dplyr::summarise(dose_24h = sum(amt), .groups = "drop") |>
  dplyr::left_join(
    sim |> dplyr::group_by(id) |> dplyr::summarise(cl = dplyr::first(cl), .groups = "drop"),
    by = "id"
  ) |>
  dplyr::mutate(AUC24_analytic = dose_24h / cl)

gate1 <- auc24 |>
  dplyr::left_join(analytic, by = c("id", "treatment")) |>
  dplyr::mutate(rel_err = AUC24 / AUC24_analytic - 1)

max_rel_err <- max(abs(gate1$rel_err))
max_rel_err
#> [1] 0.0007839027
stopifnot(max_rel_err < 0.005)

The worst per-subject discrepancy across all 800 subjects is 0.078%, which is the residual lin-up/log-down interpolation error of the observation grid, not a model discrepancy.

Gate 2 - AUC24 is exactly proportional to the total daily dose

The four regimens deliver 60, 45, 80, and 90 mg/kg/day. Because the model is linear and every arm shares one cohort, the AUC24 distributions must be exact rescalings of one another - and Jung 2024’s own eTable 4 medians are (452, 339, 603, 678), whose ratios are (1, 0.750, 1.334, 1.500) against a daily dose ratio of (1, 0.750, 1.333, 1.500). Making that a hard assertion turns a published-table property into a regression guard.

daily_dose <- regimens$mg_per_kg * 24 / regimens$tau
names(daily_dose) <- regimens$treatment

sim_median <- vapply(
  split(auc24$AUC24, auc24$treatment), median, numeric(1)
)[regimens$treatment]

gate2 <- max(abs(
  sim_median / sim_median[[1]] - daily_dose / daily_dose[[1]]
))
gate2
#> [1] 5.379163e-06
stopifnot(gate2 < 1e-4)

Comparison against the published AUC24 distribution

Jung 2024 eTable 4 reports the median [IQR] of AUC24/MIC by regimen and sampling strategy. At MIC = 1 ug/mL the ratio equals AUC24 itself, so the infinite-sampling rows are directly comparable. ncaComparisonTable() pools the simulated side by median, which is the statistic eTable 4 reports.

published <- tibble::tribble(
  ~treatment,      ~auclast,
  "15 mg/kg q6h",  452,
  "15 mg/kg q8h",  339,
  "20 mg/kg q6h",  603,
  "30 mg/kg q8h",  678
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res$result |> dplyr::filter(PPTESTCD == "auclast"),
  reference = published,
  by        = "treatment",
  params    = "auclast",
  units     = c(auclast = "ug*h/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Steady-state AUC24 versus Jung 2024 eTable 4",
                  "(MIC = 1 ug/mL, infinite-sampling rows).",
                  "* differs from reference by >20%."),
  align = c("l", "l", "r", "r", "r")
)
Steady-state AUC24 versus Jung 2024 eTable 4 (MIC = 1 ug/mL, infinite-sampling rows). * differs from reference by >20%.
NCA parameter treatment Reference Simulated % diff
AUClast (ug*h/mL) 15 mg/kg q6h 452 446 -1.3%
AUClast (ug*h/mL) 15 mg/kg q8h 339 335 -1.3%
AUClast (ug*h/mL) 20 mg/kg q6h 603 595 -1.3%
AUClast (ug*h/mL) 30 mg/kg q8h 678 669 -1.3%

The full distribution, not just the median:

auc24 |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    Simulated = sprintf("%.0f [%.0f, %.0f]", median(AUC24),
                        quantile(AUC24, 0.25), quantile(AUC24, 0.75)),
    .groups = "drop"
  ) |>
  dplyr::left_join(
    tibble::tibble(
      treatment = regimens$treatment,
      Published = c("452 [331, 630]", "339 [248, 473]",
                    "603 [442, 840]", "678 [497, 945]")
    ),
    by = "treatment"
  ) |>
  dplyr::rename(
    "Dosing regimen"                = treatment,
    "Published AUC24, median [IQR]" = Published,
    "Simulated AUC24, median [IQR]" = Simulated
  ) |>
  dplyr::relocate(`Published AUC24, median [IQR]`,
                  .before = `Simulated AUC24, median [IQR]`) |>
  knitr::kable(
    caption = paste("AUC24 (ug*h/mL) distribution versus Jung 2024 eTable 4,",
                    "MIC = 1, infinite sampling."),
    align = c("l", "r", "r")
  )
AUC24 (ug*h/mL) distribution versus Jung 2024 eTable 4, MIC = 1, infinite sampling.
Dosing regimen Published AUC24, median [IQR] Simulated AUC24, median [IQR]
15 mg/kg q6h 452 [331, 630] 446 [303, 603]
15 mg/kg q8h 339 [248, 473] 335 [228, 452]
20 mg/kg q6h 603 [442, 840] 595 [405, 803]
30 mg/kg q8h 678 [497, 945] 669 [455, 904]

Target attainment of AUC24/MIC >= 400 (Table 4) 

attainment <- auc24 |>
  tidyr::expand_grid(MIC = c(0.5, 1, 2)) |>
  dplyr::group_by(treatment, MIC) |>
  dplyr::summarise(Simulated = 100 * mean(AUC24 / MIC >= 400), .groups = "drop")

published_ta <- tibble::tribble(
  ~treatment,     ~MIC, ~Published,
  "15 mg/kg q6h", 0.5,  96.2,
  "15 mg/kg q6h", 1.0,  60.4,
  "15 mg/kg q6h", 2.0,  13.2,
  "15 mg/kg q8h", 0.5,  86.8,
  "15 mg/kg q8h", 1.0,  36.5,
  "15 mg/kg q8h", 2.0,   5.1,
  "20 mg/kg q6h", 0.5,  99.4,
  "20 mg/kg q6h", 1.0,  80.8,
  "20 mg/kg q6h", 2.0,  27.8,
  "30 mg/kg q8h", 0.5,  99.8,
  "30 mg/kg q8h", 1.0,  86.8,
  "30 mg/kg q8h", 2.0,  36.5
)

ta_cmp <- published_ta |>
  dplyr::left_join(attainment, by = c("treatment", "MIC")) |>
  dplyr::mutate(
    Difference = round(Simulated - Published, 1),
    Simulated  = round(Simulated, 1)
  )

ta_cmp |>
  dplyr::rename(
    "Dosing regimen"          = treatment,
    "MIC (ug/mL)"             = MIC,
    "Published % >= 400"      = Published,
    "Simulated % >= 400"      = Simulated,
    "Difference (percentage points)" = Difference
  ) |>
  knitr::kable(
    caption = paste("Target attainment of AUC24/MIC >= 400 versus Jung 2024",
                    "Table 4 (infinite-sampling rows)."),
    align = c("l", "r", "r", "r", "r")
  )
Target attainment of AUC24/MIC >= 400 versus Jung 2024 Table 4 (infinite-sampling rows).
Dosing regimen MIC (ug/mL) Published % >= 400 Simulated % >= 400 Difference (percentage points)
15 mg/kg q6h 0.5 96.2 94.5 -1.7
15 mg/kg q6h 1.0 60.4 56.0 -4.4
15 mg/kg q6h 2.0 13.2 12.5 -0.7
15 mg/kg q8h 0.5 86.8 80.5 -6.3
15 mg/kg q8h 1.0 36.5 31.5 -5.0
15 mg/kg q8h 2.0 5.1 6.0 0.9
20 mg/kg q6h 0.5 99.4 99.5 0.1
20 mg/kg q6h 1.0 80.8 75.0 -5.8
20 mg/kg q6h 2.0 27.8 25.5 -2.3
30 mg/kg q8h 0.5 99.8 100.0 0.2
30 mg/kg q8h 1.0 86.8 80.5 -6.3
30 mg/kg q8h 2.0 36.5 31.5 -5.0 

max_ta_gap <- max(abs(ta_cmp$Difference))
max_ta_gap
#> [1] 6.3

Every simulated attainment is below the published value by 6.3 percentage points or less, and the sign is consistently negative. That is the expected consequence of the covariate approximation: the virtual cohort’s AUC24 IQR is slightly wider than the paper’s (see the distribution table above), so slightly more subjects fall below the 400 threshold. The rank ordering of the regimens, the shape of the MIC dependence, and the paper’s clinical conclusion - that 15 mg/kg q6h and q8h leave a large minority of critically ill children under target, while 20 mg/kg q6h and 30 mg/kg q8h exceed 80% at MIC = 1 - all reproduce.

auc24 |>
  ggplot(aes(x = treatment, y = AUC24)) +
  geom_boxplot(outlier.size = 0.4) +
  geom_hline(yintercept = 400, linetype = "dashed") +
  scale_y_log10() +
  labs(
    x = NULL, y = "Steady-state AUC24 (ug*h/mL)",
    title = "AUC24 by regimen against the AUC24/MIC >= 400 target at MIC = 1",
    caption = "Dashed line: AUC24 = 400. Compare Jung 2024 Table 4 and eTable 4."
  ) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Independent cross-check: the AUC24-from-trough regression

Jung 2024’s second product is a regression predicting AUC24 from a single observed trough. The main text prints the equation as an undecoded image; its coefficients are in eTable 3:

eTable 3 parameter Estimate %RSE
Coefficient of trough concentration = 7 ug/mL 366 1.98
Additive shift for dose frequency < 8 hours 54.9 24.2
Slope of (DOSE [mg] - 500) 0.0894 17.9
Exponent of trough concentration 0.495 3.92
Residual variability (CV%) 20.2 18.8

Those four coefficients reconstruct the equation as

AUC24 = (366 + 54.9 * DFRQ8 + 0.0894 * (Dose_mg - 500)) * (Ctrough / 7)^0.495

where DFRQ8 is 1 when the dosing interval is under 8 hours. The reconstruction is confirmed against the paper’s own Table 3, which tabulates, for 12 regimen-by-weight combinations, the trough concentration that yields AUC24/MIC = 400. Feeding those published troughs back through the equation must return 400 for every row.

auc24_from_trough <- function(dose_mg, interval_h, ctrough) {
  (366 + 54.9 * as.integer(interval_h < 8) + 0.0894 * (dose_mg - 500)) *
    (ctrough / 7)^0.495
}

table3 <- tibble::tribble(
  ~treatment,     ~tau, ~WT, ~ctrough,
  "15 mg/kg q6h",  6,   10,   7.4,
  "15 mg/kg q6h",  6,   20,   6.9,
  "15 mg/kg q6h",  6,   40,   6.1,
  "15 mg/kg q8h",  8,   10,  10.0,
  "15 mg/kg q8h",  8,   20,   9.3,
  "15 mg/kg q8h",  8,   40,   8.0,
  "20 mg/kg q6h",  6,   10,   7.2,
  "20 mg/kg q6h",  6,   20,   6.6,
  "20 mg/kg q6h",  6,   40,   5.6,
  "30 mg/kg q8h",  8,   10,   9.3,
  "30 mg/kg q8h",  8,   20,   8.0,
  "30 mg/kg q8h",  8,   40,   6.1
) |>
  dplyr::mutate(
    mg_per_kg = as.numeric(sub(" .*", "", treatment)),
    dose_mg   = mg_per_kg * WT,
    AUC24     = round(auc24_from_trough(dose_mg, tau, ctrough), 1)
  )

# Table 3 troughs are printed to one decimal, so the round trip lands near 400
# rather than exactly on it; 2 ug*h/mL is the rounding budget.
stopifnot(max(abs(table3$AUC24 - 400)) < 2)

table3 |>
  dplyr::select(treatment, WT, ctrough, AUC24) |>
  dplyr::rename(
    "Dosing regimen"               = treatment,
    "Body weight (kg)"             = WT,
    "Published trough (ug/mL)"     = ctrough,
    "Regression AUC24 (ug*h/mL)"   = AUC24
  ) |>
  knitr::kable(
    caption = paste("Jung 2024 Table 3 round trip: the published trough",
                    "achieving AUC24/MIC = 400 at MIC = 1, fed back through the",
                    "eTable 3 regression, returns 400."),
    align = c("l", "r", "r", "r")
  )
Jung 2024 Table 3 round trip: the published trough achieving AUC24/MIC = 400 at MIC = 1, fed back through the eTable 3 regression, returns 400.
Dosing regimen Body weight (kg) Published trough (ug/mL) Regression AUC24 (ug*h/mL)
15 mg/kg q6h 10 7.4 400.5
15 mg/kg q6h 20 6.9 400.2
15 mg/kg q6h 40 6.1 401.5
15 mg/kg q8h 10 10.0 399.3
15 mg/kg q8h 20 9.3 400.7
15 mg/kg q8h 40 8.0 400.6
20 mg/kg q6h 10 7.2 399.6
20 mg/kg q6h 20 6.6 400.1
20 mg/kg q6h 40 5.6 400.9
30 mg/kg q8h 10 9.3 400.7
30 mg/kg q8h 20 8.0 400.6
30 mg/kg q8h 40 6.1 400.4

The two products can now be checked against each other over the whole virtual cohort, which is how the regression was actually built: each subject’s steady-state trough (PKNCA cmin over the dosing window) is fed through the regression, and the result compared with that subject’s own model AUC24. No extra simulation is needed - both quantities come from the PKNCA run above.

cross <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "cmin", start == 0, end == 24) |>
  dplyr::select(id, treatment, ctrough = PPORRES) |>
  dplyr::left_join(auc24, by = c("id", "treatment")) |>
  dplyr::left_join(
    events |> dplyr::filter(evid == 1) |>
      dplyr::group_by(id) |>
      dplyr::summarise(dose_mg = dplyr::first(amt),
                       tau = 24 / dplyr::n(), .groups = "drop"),
    by = "id"
  ) |>
  dplyr::mutate(
    AUC24_regression = auc24_from_trough(dose_mg, tau, ctrough),
    pct_diff = 100 * (AUC24_regression / AUC24 - 1)
  )

pooled_bias <- median(cross$pct_diff)
pooled_cv <- 100 * sd(cross$AUC24_regression / cross$AUC24) /
  mean(cross$AUC24_regression / cross$AUC24)

cross |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    `Median % difference` = round(median(pct_diff), 1),
    `IQR of % difference` = sprintf("%.1f to %.1f", quantile(pct_diff, 0.25),
                                    quantile(pct_diff, 0.75)),
    .groups = "drop"
  ) |>
  dplyr::rename("Dosing regimen" = treatment) |>
  knitr::kable(
    caption = paste("Regression AUC24 (from each subject's simulated",
                    "steady-state trough) versus that subject's model AUC24,",
                    "by regimen."),
    align = c("l", "r", "r")
  )
Regression AUC24 (from each subject’s simulated steady-state trough) versus that subject’s model AUC24, by regimen.
Dosing regimen Median % difference IQR of % difference
15 mg/kg q6h 11.9 1.4 to 26.1
15 mg/kg q8h 4.5 -5.6 to 18.1
20 mg/kg q6h -0.1 -10.1 to 11.6
30 mg/kg q8h -20.6 -28.9 to -8.4

Three things fall out, and all three are properties of the paper rather than of this extraction:

  • The reconstruction is unbiased in aggregate. Pooled across the cohort the median difference is +0.6%, matching Jung 2024’s own claim that the equation “provided an unbiased prediction of AUC24 across the observed range of trough vancomycin concentrations.”
  • Its scatter matches the published residual error. The ratio of the two AUC24 routes has a CV of 22.1%, against the 20.2% residual variability reported for the regression in eTable 3. Recovering that number independently is a third confirmation - alongside the Table 3 round trip and the unbiasedness - that the eTable 3 coefficients were reassembled into the right functional form.
  • But the bias is regimen-dependent, and worst where the paper extrapolates. The per-regimen medians run from about +12% for 15 mg/kg q6h to about -21% for 30 mg/kg q8h. The regression carries dose as a linear additive term centred at 500 mg and was fit to the observed regimens only (5-31 mg/kg q6h, most commonly 15 mg/kg every 6-9 hours). 30 mg/kg q8h is a purely simulated regimen, reaching 1200 mg per dose in the heaviest children - well outside that calibration range - so the linear dose term is being extrapolated. A user who wants AUC24 for one of the higher simulated regimens should take it from the PK model (Dose * (24/tau) / CL), not from the trough regression.

Assumptions and deviations

  • Covariate joint distribution. Only marginal medians and IQRs are published (Table 1), so the virtual cohort uses log-normal marginals matched to those and treats body weight and eGFR as independent. The real cohort almost certainly has some weight-eGFR association. Because AUC24 scales as WT^0.25 * eGFR^-0.5259, a positive association would partly cancel and narrow the AUC24 distribution - which is the direction of the residual disagreement observed above (simulated IQR slightly wider, target attainment 4-6 percentage points lower than published). No parameter was adjusted to close that gap.
  • Cohort size. 200 subjects per arm rather than the paper’s 2250, per the repository cap. The comparison statistics are medians, IQRs, and proportions, all of which are stable at n = 200; the cohort is deterministic (exact marginal quantiles, seeded coupling) so the reported numbers do not move under reseeding.
  • Time-varying covariates held constant. Jung 2024 treats both body weight and eGFR as time-varying (eTable 1) and fit them that way. The target- attainment simulation resampled baseline covariate vectors, so this vignette holds both constant over the 24-hour steady-state window, matching the paper’s own simulation rather than its estimation dataset.
  • Steady state imposed, not simulated. ss = 1 replaces a multi-day burn-in. This matches the paper’s analytic Dose * (24/tau) / CL definition of AUC24 exactly (Gate 1) and avoids the long-half-life tail that the 102.6% CV on Vc would otherwise produce in a finite burn-in.
  • Body weight is not the dosing basis in the model file. The model takes amt in mg; the mg/kg regimens are converted to mg in the vignette’s event table using each subject’s WT.
  • eGFR equation. The packaged model uses the bedside Schwartz eGFR that Jung 2024 Table 2 reports. The paper also refit the model with the CKiD cystatin-C equation in the 248-subject subset having both serum creatinine and cystatin C, and found equivalent predictive performance (eFigures 8-9), but no separate parameter table is published for that fit, so it is not packaged as a second model.
  • The AUC24-from-trough regression is not packaged as a model file. It is a nonlinear regression on model-predicted AUC24, not an ODE model, so it does not belong in inst/modeldb/. Its coefficients are reproduced in this vignette from eTable 3 and validated three independent ways (Table 3 round trip, aggregate unbiasedness, and recovery of the published 20.2% residual CV).
  • The trough regression should not be used for the higher simulated regimens. Its dose term is linear and centred at 500 mg and it was fit only to the observed regimens, so it disagrees with the PK model by a median of about -21% at 30 mg/kg q8h (see the cross-check section). This is a property of Jung 2024, not of the packaged model, and no parameter was altered because of it; it is recorded here so a downstream user does not mistake the regression for a general-purpose AUC24 estimator.
  • The extracorporeal-therapy cohort is out of scope. Jung 2024 recruited 33 additional subjects on ECMO / CRRT / extracorporeal liver support and explicitly excluded them from this analysis (eMethods); no model was published for them.
  • No non-paper-derived parameter values. Every value in ini() comes from Jung 2024 Table 2 or the supplement’s eMethods; nothing was digitised from a figure, obtained by correspondence, or carried from another model. The supplement was retrieved from the publisher’s own permalink (http://links.lww.com/CCX/B406) and is archived alongside the lead PDF.