Skip to contents

Vancomycin (Yu 2023)

Source: vignettes/articles/Yu_2023_vancomycin.Rmd
Yu_2023_vancomycin.Rmd

Model and source

  • Citation: Yu Z, Liu J, Yu H, Zhou L, Zhu J, Liang G, Yang Y, Zheng Y, Han Y, Xu J, Han G, Yu L, Zhao Y. Population pharmacokinetics and individualized dosing of vancomycin for critically ill patients receiving continuous renal replacement therapy: the role of residual diuresis. Front Pharmacol. 2023;14:1298397. doi:10.3389/fphar.2023.1298397
  • Description: One-compartment IV population PK model for vancomycin in 71 critically ill adults receiving intermittent intravenous vancomycin during continuous renal replacement therapy (Yu 2023). Total clearance is a power function of 24-hour residual diuresis, CL = 1.05 * 1.90^(log10(URINE_VOL_24H + 10) / 2.3) L/h, so an anuric subject retains only the CRRT-mediated clearance (1.39 L/h) while a subject producing 3 L/day reaches 2.77 L/h. Central volume is 69.0 L with no interindividual variability, because only trough concentrations were available. Age, sex, body weight, BMI, daily dose, serum creatinine, blood urea nitrogen, and CRRT modality were screened but not retained.
  • Article: https://doi.org/10.3389/fphar.2023.1298397
  • Supplement (Table S1, covariate screening; Figure S1, covariate correlation): https://www.frontiersin.org/articles/10.3389/fphar.2023.1298397/full#supplementary-material

Yu 2023 asked a narrow question: when a critically ill patient is on continuous renal replacement therapy (CRRT), does whatever urine they still make matter for vancomycin clearance? Earlier CRRT models largely ignored residual diuresis, and the paper’s argument is that this is why they extrapolate poorly. Serum creatinine and blood urea nitrogen are unreliable in patients with acute kidney injury on CRRT, whereas a 24-hour urine volume is free, routinely charted, and directly interpretable. The final model keeps exactly one covariate: the logarithm of the 24-hour urine volume, on clearance.

Population

Four Chinese ICUs (Sir Run Run Shaw Hospital and the Second Affiliated Hospital, Zhejiang University School of Medicine; Affiliated Xiaoshan Hospital, Hangzhou Normal University; Zhejiang Zhoushan Hospital) contributed patients retrospectively between January 2019 and October 2022. Of 191 trough concentrations from 101 patients, 113 troughs from 71 patients survived the exclusions (age below 18, missing weight / height / dosing / laboratory data, non-continuous CRRT, or sampling more than 48 h after the last dose).

The cohort was mostly male (44/71, 62.0%), mean age 61.6 years (SD 14.6, range 24-87), mean weight 64.2 kg (SD 14.7, range 32.1-100), mean BMI 23.3 kg/m^2 (SD 3.97). CRRT modality was CVVH in 40/71 (56.3%) and CVVHDF in 31/71 (43.7%); CRRT settings themselves were not recorded, which the paper lists as a limitation. The median total daily vancomycin dose was 1000 mg (IQR 500-2000, range 500-3000), i.e. a median 15.4 mg/kg/day.

The covariate that matters is heavily skewed toward anuria: median 24-hour urine volume 160 mL (IQR 7.00-780, range 0.00-6220), with 33/71 (46.5%) at or below 100 mL and only 6/71 (8.45%) at or above 2500 mL (Yu 2023 Table 1). Nearly half the cohort is therefore at or near the anuric end of the covariate range, and the model’s +10 mL offset exists to keep the log term finite for them.

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

Source trace

The per-parameter origin is also recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Yu_2023_vancomycin.R.

Equation / parameter Value Source location
Structural model: one compartment, first-order elimination, IV n/a Yu 2023 Results, “PopPK model development” (AIC 800.588 one-compartment vs 803.682 two-compartment)
CL (L/h) = 1.05 * 1.90^(LOG(UV + 10) / 2.3) n/a Yu 2023 Table 2 footnote (“The final model was as follows”)
lcl = log(1.05) 1.05 L/h Yu 2023 Table 2, row “CL (L/h)”, Estimate 1.05 (RSE 17.2%); bootstrap median 1.07, 95% CI 0.721-1.53
lvc = log(69.0) 69.0 L Yu 2023 Table 2, row “V (L)”, Estimate 69.0 (RSE 6.61%); bootstrap median 68.6, 95% CI 59.9-77.7
e_urine_vol_24h_cl 1.90 Yu 2023 Table 2, row “theta UV-CL”, Estimate 1.90 (RSE 16.5%); bootstrap median 1.91, 95% CI 1.31-2.71
etalcl ~ 0.121 omega^2 = 0.121 Yu 2023 Table 2, row “omega (%) for CL”, 12.1 (RSE 47.9%); the Table 2 abbreviation footnote defines omega as the “interindividual variance for CL”
propSd = sqrt(0.0978) sigma = 0.3127 Yu 2023 Table 2, row “sigma (%)”, 9.78 (RSE 29.1%), same variance-as-percent scale
No IIV on V n/a Yu 2023 Table 2 footnote: “Intre-individual variability for V was not estimated”
Covariate retained: 24-h urine volume on CL only n/a Yu 2023 Results; Supplementary Table S1 (forward inclusion / backward elimination)
Screened but not retained: BMI, SCr, CRRT modality (on CL); weight, BMI, sex (on V) n/a Yu 2023 Supplementary Table S1, steps 3-11
Simulated dosing regimens and urine volumes n/a Yu 2023 Methods, “Simulation and dosing regimen optimization”
Published PTA curves n/a Yu 2023 Figure 3 (A: AUC >= 400 mgh/L; B: AUC 400-600 mgh/L)

Which logarithm?

The Table 2 footnote prints the final model as

CL (L/h) = 1.05 * 1.90^(LOG(UV+10)/2.3), V (L) = 69.0

LOG is not defined anywhere in the paper, and the two candidate readings are not close to each other:

  • Base-10. 2.3 is then the median-normalising constant log10(160 + 10) = 2.230, rounded to two significant figures. This matches Yu 2023 Methods, which says “the effects of continuous covariates were modeled using a median-normalized model”. Written as a plain power of the offset urine volume this is CL = 1.05 * (UV + 10)^0.1212.
  • Natural log, the NONMEM convention, in which case 2.3 is ln(10) = 2.3026 and the exponent collapses to log10(UV + 10) with no normalisation. Written as a plain power this is CL = 1.05 * (UV + 10)^0.2791.

Both are “power” covariate models, so Supplementary Table S1’s “power” functional form does not separate them. The clearances they imply do:

cl_log10 <- function(uv) 1.05 * 1.90^(log10(uv + 10) / 2.3)  # as implemented
cl_ln    <- function(uv) 1.05 * 1.90^(log(uv + 10) / 2.3)    # natural-log reading

uv_grid <- c(0, 100, 500, 1500, 2000, 3000)
tibble(
  `UV (mL/24h)`            = uv_grid,
  `CL, base-10 (L/h)`      = round(cl_log10(uv_grid), 3),
  `CL, natural log (L/h)`  = round(cl_ln(uv_grid), 3)
) |>
  knitr::kable(caption = "Typical clearance under the two readings of `LOG`.")
Typical clearance under the two readings of LOG.
UV (mL/24h) CL, base-10 (L/h) CL, natural log (L/h)
0 1.388 1.996
100 1.856 3.898
500 2.235 5.981
1500 2.550 8.097
2000 2.640 8.770
3000 2.772 9.816

The paper adjudicates this itself. Yu 2023 Figure 3A plots the probability of attaining an AUC of at least 400 mgh/L for five regimens across the same six urine volumes, and the 1500 mg q12h and 1000 mg q12h curves sit at or above the 90% line at every* urine volume. Because this model is linear and one-compartment, the steady-state AUC0-24 is exactly daily dose / CL (the PKNCA block below confirms this to within 0.03%), so the counterfactual is a one-liner:

omega <- sqrt(0.121)
n_cohort <- 200L
# Deterministic stratified representation of the log-normal CL distribution:
# the i-th of n equally spaced probability points. No RNG is involved anywhere
# in this vignette, so every number below is reproducible on any machine and any
# thread count (see known-vignette-failure-patterns.md pattern 12).
etas <- stats::qnorm((seq_len(n_cohort) - 0.5) / n_cohort, mean = 0, sd = omega)

pta_closed <- function(cl_fun, daily_dose, uv, lo = 400, hi = Inf) {
  auc <- daily_dose / (cl_fun(uv) * exp(etas))
  100 * mean(auc >= lo & auc <= hi)
}

counterfactual <- tidyr::expand_grid(
  regimen = c("1000 mg q12h", "1500 mg q12h"),
  uv      = uv_grid
) |>
  mutate(
    daily = if_else(regimen == "1000 mg q12h", 2000, 3000),
    `PTA, base-10 (%)`     = mapply(function(d, u) pta_closed(cl_log10, d, u), daily, uv),
    `PTA, natural log (%)` = mapply(function(d, u) pta_closed(cl_ln,    d, u), daily, uv)
  ) |>
  select(regimen, `UV (mL/24h)` = uv, `PTA, base-10 (%)`, `PTA, natural log (%)`)

knitr::kable(
  counterfactual,
  caption = "PTA for AUC >= 400 mg*h/L. Yu 2023 Figure 3A puts both of these regimens at or above 90% at every urine volume."
)
PTA for AUC >= 400 mg*h/L. Yu 2023 Figure 3A puts both of these regimens at or above 90% at every urine volume.
regimen UV (mL/24h) PTA, base-10 (%) PTA, natural log (%)
1000 mg q12h 0 100.0 99.5
1000 mg q12h 100 100.0 76.5
1000 mg q12h 500 99.0 30.5
1000 mg q12h 1500 97.5 8.5
1000 mg q12h 2000 96.5 5.5
1000 mg q12h 3000 95.5 2.5
1500 mg q12h 0 100.0 100.0
1500 mg q12h 100 100.0 97.0
1500 mg q12h 500 100.0 74.0
1500 mg q12h 1500 100.0 41.5
1500 mg q12h 2000 100.0 32.5
1500 mg q12h 3000 100.0 22.0

The natural-log reading puts 1500 mg q12h – the largest regimen the paper simulated – at roughly a fifth of its published attainment in a patient making 3 L of urine a day, and 1000 mg q12h below 5%. The base-10 reading reproduces the published saturation. The model file therefore implements base-10, and the gate below asserts both halves of this (the implemented reading clears 90%, the rejected reading does not).

lo_impl <- min(counterfactual$`PTA, base-10 (%)`)
lo_ln   <- min(counterfactual$`PTA, natural log (%)`)
stopifnot(
  # Yu 2023 Figure 3A: 1000 and 1500 mg q12h are above the 90% dashed line at
  # every simulated urine volume.
  lo_impl >= 90,
  # ... and the rejected reading is nowhere near it, so the check can go red.
  lo_ln < 50
)
c(`implemented (base-10) minimum PTA %` = round(lo_impl, 1),
  `rejected (natural log) minimum PTA %` = round(lo_ln, 1))
#>  implemented (base-10) minimum PTA % rejected (natural log) minimum PTA % 
#>                                 95.5                                  2.5

What scale are omega and sigma on?

Table 2 reports omega (%) for CL = 12.1 and sigma (%) = 9.78, and its abbreviation footnote calls omega the “interindividual variance for CL”. Read as variances scaled by 100 that gives omega = sqrt(0.121) = 0.348 (35.9% CV) and sigma = sqrt(0.0978) = 0.313 (31.3% proportional error); read as standard deviations it would give 0.121 and 0.0978. The difference is nearly threefold, so it has to be settled rather than assumed. Two features of the paper’s own figures settle it, both toward the variance reading.

Figure 2I, the ETA histogram, spans roughly -0.9 to +0.5 with a peak density near 1.5. An eta with SD 0.348 has a density of 1 / (0.348 * sqrt(2 * pi)) = 1.15 at zero and reaches +/-0.9 at 2.6 SD; an eta with SD 0.121 would have a peak density of 3.3 and would never reach 0.9 (7.4 SD).

Figure 3A’s 500 mg q12h curve gives a quantitative estimate. That curve is the only one of the five that stays well inside (0, 100) across the whole urine range, so qnorm(PTA) is well conditioned along it. Since AUC = daily dose / CL and CL is log-normal, qnorm(PTA) = [ln(1000 / 400) - ln(CL_typ)] / omega, and regressing the published qnorm(PTA) on ln(1000 / (400 * CL_typ)) recovers 1 / omega as the slope.

# Digitised by the operator from Yu 2023 Figure 3A, 500 mg q12h series (filled
# circles). Read against the 20-unit gridlines; roughly +/-3 percentage points.
# Used ONLY as a comparison reference and for this slope check -- no parameter
# value in the model file comes from a figure.
fig3a_500q12h <- tibble(
  uv  = uv_grid,
  pta = c(88, 64, 41, 29, 21.5, 22)
)

slope_fit <- fig3a_500q12h |>
  mutate(
    x = log(1000 / (400 * cl_log10(uv))),
    z = stats::qnorm(pta / 100)
  ) |>
  with(stats::lm(z ~ x))

omega_recovered <- 1 / stats::coef(slope_fit)[["x"]]

tibble(
  Quantity = c("omega implied by Figure 3A slope",
               "omega = sqrt(0.121)  (variance reading, implemented)",
               "omega = 0.121        (standard-deviation reading, rejected)"),
  Value = round(c(omega_recovered, sqrt(0.121), 0.121), 4)
) |>
  knitr::kable(caption = "Recovering omega from the published PTA curve.")
Recovering omega from the published PTA curve.
Quantity Value
omega implied by Figure 3A slope 0.3434
omega = sqrt(0.121) (variance reading, implemented) 0.3479
omega = 0.121 (standard-deviation reading, rejected) 0.1210 

stopifnot(
  # The digitised curve recovers the variance reading and excludes the
  # standard-deviation reading by a wide margin. The window is deliberately
  # loose (digitisation is +/-3 PTA points) but 0.121 is nowhere near it.
  omega_recovered > 0.25, omega_recovered < 0.45
)

Both parameters are therefore encoded as variances: etalcl ~ 0.121 and propSd <- sqrt(0.0978).

Virtual cohort

Original observed data are not publicly available. The cohort below reproduces the design of the paper’s own Monte Carlo simulation (Yu 2023 Methods, “Simulation and dosing regimen optimization”): five maintenance regimens crossed with six residual-diuresis scenarios, 200 subjects per arm.

The only random effect in this model is etalcl, so instead of drawing it the cohort uses the 200 equally spaced quantiles of N(0, 0.121) defined in the counterfactual chunk above. That is a deterministic stratified representation of the same log-normal clearance distribution: it removes Monte Carlo noise from every PTA below and, because rxode2’s parallel RNG streams are partitioned per solver thread, it also makes every result here identical on a 2-thread CI runner and a 16-thread workstation.

mod <- readModelDb("Yu_2023_vancomycin")

regimens <- tibble::tribble(
  ~regimen,        ~amt,  ~ii,
  "500 mg qd",      500,   24,
  "500 mg q12h",    500,   12,
  "750 mg q12h",    750,   12,
  "1000 mg q12h",  1000,   12,
  "1500 mg q12h",  1500,   12
) |>
  mutate(daily = amt * 24 / ii)

arms <- tidyr::expand_grid(regimens, uv = uv_grid) |>
  mutate(
    treatment = paste0(regimen, " | UV ", uv),
    id_offset = (dplyr::row_number() - 1L) * n_cohort
  )

# Steady state is imposed with ss = 1 rather than by simulating out to it. The
# slowest subject here (anuric, eta at the 0.25th percentile) has a half-life of
# about 70 h, so an explicit run-in would need three weeks of dosing; ss = 1 is
# exact and free. Doses are given as 1-hour infusions, the standard vancomycin
# administration -- AUC0-24 at steady state is identical for a bolus, so nothing
# below depends on that choice.
make_arm <- function(regimen, amt, ii, daily, uv, treatment, id_offset) {
  ev <- as.data.frame(
    rxode2::et(amt = amt, ii = ii, ss = 1, addl = 24 / ii - 1,
               dur = 1, cmt = "central") |>
      rxode2::et(seq(0, 24, by = 0.5), cmt = "central")
  )
  # rxode2::et() omits the `addl` column entirely when addl is 0 (the qd arm),
  # which would make the per-arm frames non-conformable for rbind(). Restore it.
  if (is.null(ev$addl)) ev$addl <- 0L
  ev <- ev[, c("time", "cmt", "amt", "ii", "addl", "evid", "ss", "dur")]
  out <- ev[rep(seq_len(nrow(ev)), n_cohort), ]
  out$id            <- rep(seq_len(n_cohort), each = nrow(ev)) + id_offset
  out$etalcl        <- rep(etas, each = nrow(ev))
  out$URINE_VOL_24H <- uv
  out$treatment     <- treatment
  out
}

events <- do.call(
  rbind,
  lapply(seq_len(nrow(arms)), function(i) do.call(make_arm, as.list(arms[i, ])))
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

# One dose record per administered dose, for PKNCA.
doses <- do.call(rbind, lapply(seq_len(nrow(arms)), function(i) {
  a <- arms[i, ]
  tidyr::expand_grid(
    id   = seq_len(n_cohort) + a$id_offset,
    time = seq(0, 24 - a$ii, by = a$ii)
  ) |>
    mutate(amt = a$amt, treatment = a$treatment)
}))

c(arms = nrow(arms), subjects = nrow(arms) * n_cohort, event_rows = nrow(events))
#>       arms   subjects event_rows 
#>         30       6000     300000

Simulation 

sim <- rxode2::rxSolve(
  mod, events = events, returnType = "data.frame",
  keep = "treatment"
)
#> ℹ parameter labels from comments will be replaced by 'label()'
sim |>
  filter(!is.na(Cc), grepl("^1000 mg q12h", treatment)) |>
  group_by(URINE_VOL_24H, 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() +
  facet_wrap(~ URINE_VOL_24H, nrow = 2,
             labeller = labeller(URINE_VOL_24H = function(x) paste0("UV ", x, " mL/24h"))) +
  labs(
    x = "Time within the steady-state day (h)", y = "Vancomycin (mg/L)",
    title = "1000 mg q12h at steady state, by residual diuresis",
    caption = "Median and 5th-95th percentile of the 200-subject deterministic cohort."
  )

The visible spread of the steady-state trough across urine volumes is the whole point of the paper: an anuric patient on 1000 mg q12h sits around 53.6 mg/L at trough, while the same regimen in a patient making 3 L/day sits around 23.9 mg/L.

PKNCA validation

AUC0-24 at steady state is computed with PKNCA for each of the 30 arms. Both dosing intervals used here (12 h and 24 h) divide 24 h exactly, so AUC0-24 is phase-independent.

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

# PKNCA anchors AUC0-24 at time 0 and warns once per subject if the
# concentration frame has no such row. The observation grid produces one for
# every subject; assert it rather than assume it.
stopifnot(
  sim_nca |>
    dplyr::group_by(treatment, id) |>
    dplyr::summarise(has0 = any(time == 0), .groups = "drop") |>
    dplyr::pull(has0) |>
    all()
)
sim_nca <- dplyr::arrange(sim_nca, id, time)

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

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

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

Instrument check: the NCA against the closed form

For a linear one-compartment model at steady state, AUC0-tau = dose / CL exactly, so AUC0-24 = daily dose / CL. Every element of the pipeline – dose amount, ss = 1 handling, the compartment the dose lands in, the volume, the observation grid, and PKNCA’s trapezoidal integration – has to be right for that identity to hold. This is a deterministic numerical check, so it is asserted tightly.

auc_subject <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD == "auclast") |>
  select(treatment, id, auc = PPORRES)

cl_subject <- sim |>
  distinct(id, treatment, URINE_VOL_24H, cl)

daily_by_arm <- doses |>
  group_by(treatment, id) |>
  summarise(daily = sum(amt), .groups = "drop")

auc_subject <- auc_subject |>
  left_join(cl_subject, by = c("id", "treatment")) |>
  left_join(daily_by_arm, by = c("id", "treatment")) |>
  mutate(auc_closed = daily / cl,
         pct_diff = 100 * (auc - auc_closed) / auc_closed)

worst_pct <- max(abs(auc_subject$pct_diff))
stopifnot(
  nrow(auc_subject) == nrow(arms) * n_cohort,
  !anyNA(auc_subject$auc),
  worst_pct < 0.5
)
c(`subjects checked` = nrow(auc_subject),
  `worst |AUC0-24 vs dose/CL| (%)` = round(worst_pct, 4))
#>               subjects checked worst |AUC0-24 vs dose/CL| (%) 
#>                       6.00e+03                       2.23e-02

Steady-state exposure by arm 

as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "cmin", "cav", "auclast")) |>
  group_by(treatment, PPTESTCD) |>
  summarise(median = median(PPORRES), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = median) |>
  left_join(arms |> select(treatment, regimen, uv), by = "treatment") |>
  arrange(match(regimen, regimens$regimen), uv) |>
  transmute(
    Regimen              = regimen,
    `UV (mL/24h)`        = uv,
    `Cmax,ss (mg/L)`     = round(cmax, 1),
    `Cmin,ss (mg/L)`     = round(cmin, 1),
    `Cav,ss (mg/L)`      = round(cav, 1),
    `AUC0-24 (mg*h/L)`   = round(auclast, 0)
  ) |>
  knitr::kable(caption = "Median steady-state NCA parameters from PKNCA, by regimen and residual diuresis.")
Median steady-state NCA parameters from PKNCA, by regimen and residual diuresis.
Regimen UV (mL/24h) Cmax,ss (mg/L) Cmin,ss (mg/L) Cav,ss (mg/L) AUC0-24 (mg*h/L)
500 mg qd 0 18.7 11.8 15.0 360
500 mg qd 100 15.0 8.1 11.2 269
500 mg qd 500 13.2 6.3 9.3 224
500 mg qd 1500 12.1 5.2 8.2 196
500 mg qd 2000 11.8 4.9 7.9 189
500 mg qd 3000 11.5 4.6 7.5 180
500 mg q12h 0 33.5 26.8 30.0 720
500 mg q12h 100 25.9 19.3 22.4 539
500 mg q12h 500 22.1 15.5 18.6 447
500 mg q12h 1500 19.9 13.2 16.3 392
500 mg q12h 2000 19.3 12.7 15.8 379
500 mg q12h 3000 18.6 11.9 15.0 361
750 mg q12h 0 50.2 40.2 45.0 1081
750 mg q12h 100 38.9 28.9 33.7 808
750 mg q12h 500 33.2 23.3 28.0 671
750 mg q12h 1500 29.8 19.8 24.5 588
750 mg q12h 2000 29.0 19.0 23.7 568
750 mg q12h 3000 27.9 17.9 22.5 541
1000 mg q12h 0 66.9 53.6 60.0 1441
1000 mg q12h 100 51.8 38.6 44.9 1078
1000 mg q12h 500 44.3 31.0 37.3 895
1000 mg q12h 1500 39.7 26.5 32.7 784
1000 mg q12h 2000 38.6 25.4 31.6 758
1000 mg q12h 3000 37.1 23.9 30.1 721
1500 mg q12h 0 100.4 80.4 90.1 2161
1500 mg q12h 100 77.7 57.8 67.3 1616
1500 mg q12h 500 66.4 46.5 55.9 1342
1500 mg q12h 1500 59.6 39.7 49.0 1177
1500 mg q12h 2000 57.9 38.0 47.4 1137
1500 mg q12h 3000 55.7 35.8 45.1 1082

Replicating Figure 3 

pta <- auc_subject |>
  group_by(treatment) |>
  summarise(
    pta_400   = 100 * mean(auc >= 400),
    pta_400600 = 100 * mean(auc >= 400 & auc <= 600),
    .groups = "drop"
  ) |>
  left_join(arms |> select(treatment, regimen, uv), by = "treatment") |>
  mutate(regimen = factor(regimen, levels = regimens$regimen)) |>
  arrange(regimen, uv)
# Replicates Figure 3 of Yu 2023: PTA vs 24-hour urine volume for five
# maintenance regimens. (A) AUC >= 400 mg*h/L; (B) AUC between 400 and 600.
pta |>
  tidyr::pivot_longer(c(pta_400, pta_400600), names_to = "panel", values_to = "pta") |>
  mutate(panel = recode(panel,
                        pta_400    = "A: AUC >= 400 mg*h/L",
                        pta_400600 = "B: AUC 400-600 mg*h/L")) |>
  ggplot(aes(factor(uv), pta, group = regimen, shape = regimen, linetype = regimen)) +
  geom_hline(yintercept = 90, linetype = "dotted") +
  geom_line() +
  geom_point(size = 2) +
  facet_wrap(~ panel) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(
    x = "24-hour urine volume (mL)", y = "PTA (%)",
    shape = NULL, linetype = NULL,
    title = "Probability of target attainment",
    caption = "Replicates Figure 3 of Yu 2023. Dotted line: 90% PTA."
  ) +
  theme(legend.position = "bottom")

Comparison against the published Figure 3A 

# Digitised by the operator from Yu 2023 Figure 3A against the 20-unit
# gridlines; roughly +/-3 percentage points, and the curves that sit on the
# 100% ceiling are recorded as 100.
fig3a_published <- tibble::tribble(
  ~regimen,       ~uv,  ~published,
  "500 mg qd",       0,   24,
  "500 mg qd",     100,  4.5,
  "500 mg qd",     500,  1.5,
  "500 mg qd",    1500,  0.5,
  "500 mg qd",    2000,  0.5,
  "500 mg qd",    3000,  0.5,
  "500 mg q12h",     0,   88,
  "500 mg q12h",   100,   64,
  "500 mg q12h",   500,   41,
  "500 mg q12h",  1500,   29,
  "500 mg q12h",  2000, 21.5,
  "500 mg q12h",  3000,   22,
  "750 mg q12h",     0,  100,
  "750 mg q12h",   100,   99,
  "750 mg q12h",   500,   88,
  "750 mg q12h",  1500,   82,
  "750 mg q12h",  2000,   80,
  "750 mg q12h",  3000,   75,
  "1000 mg q12h",    0,  100,
  "1000 mg q12h",  100,  100,
  "1000 mg q12h",  500,   99,
  "1000 mg q12h", 1500,   98,
  "1000 mg q12h", 2000, 96.5,
  "1000 mg q12h", 3000,   96,
  "1500 mg q12h",    0,  100,
  "1500 mg q12h",  100,  100,
  "1500 mg q12h",  500,  100,
  "1500 mg q12h", 1500,  100,
  "1500 mg q12h", 2000,  100,
  "1500 mg q12h", 3000,  100
)

cmp <- pta |>
  select(regimen, uv, simulated = pta_400) |>
  mutate(regimen = as.character(regimen)) |>
  left_join(fig3a_published, by = c("regimen", "uv")) |>
  mutate(diff = simulated - published)

cmp |>
  transmute(
    Regimen                = regimen,
    `UV (mL/24h)`          = uv,
    `Simulated PTA (%)`    = round(simulated, 1),
    `Yu 2023 Fig 3A (%)`   = published,
    `Difference (points)`  = round(diff, 1)
  ) |>
  knitr::kable(
    caption = "Simulated PTA for AUC >= 400 mg*h/L against the digitised Figure 3A."
  )
Simulated PTA for AUC >= 400 mg*h/L against the digitised Figure 3A.
Regimen UV (mL/24h) Simulated PTA (%) Yu 2023 Fig 3A (%) Difference (points)
500 mg qd 0 38.0 24.0 14.0
500 mg qd 100 13.0 4.5 8.5
500 mg qd 500 4.5 1.5 3.0
500 mg qd 1500 2.0 0.5 1.5
500 mg qd 2000 1.5 0.5 1.0
500 mg qd 3000 1.0 0.5 0.5
500 mg q12h 0 95.5 88.0 7.5
500 mg q12h 100 80.5 64.0 16.5
500 mg q12h 500 62.5 41.0 21.5
500 mg q12h 1500 47.5 29.0 18.5
500 mg q12h 2000 44.0 21.5 22.5
500 mg q12h 3000 38.5 22.0 16.5
750 mg q12h 0 100.0 100.0 0.0
750 mg q12h 100 98.0 99.0 -1.0
750 mg q12h 500 93.0 88.0 5.0
750 mg q12h 1500 86.5 82.0 4.5
750 mg q12h 2000 84.5 80.0 4.5
750 mg q12h 3000 80.5 75.0 5.5
1000 mg q12h 0 100.0 100.0 0.0
1000 mg q12h 100 100.0 100.0 0.0
1000 mg q12h 500 99.0 99.0 0.0
1000 mg q12h 1500 97.5 98.0 -0.5
1000 mg q12h 2000 96.5 96.5 0.0
1000 mg q12h 3000 95.5 96.0 -0.5
1500 mg q12h 0 100.0 100.0 0.0
1500 mg q12h 100 100.0 100.0 0.0
1500 mg q12h 500 100.0 100.0 0.0
1500 mg q12h 1500 100.0 100.0 0.0
1500 mg q12h 2000 100.0 100.0 0.0
1500 mg q12h 3000 100.0 100.0 0.0

The three larger regimens reproduce closely: across the eighteen 750 / 1000 / 1500 mg q12h points the largest gap is 5.5 percentage points. The two 500 mg regimens sit systematically above the published curve, by up to 22.5 points. This is a real, reproducible disagreement between the paper’s printed final-model equation and the paper’s own Monte Carlo, and it is recorded as a deviation rather than tuned away – see “Assumptions and deviations”.

Published claims 

claim <- function(text, achieved, pass, deviation = FALSE) {
  tibble(Claim = text, Achieved = achieved, Pass = pass, Deviation = deviation)
}

pta_a <- function(reg, uvs = uv_grid) {
  v <- pta$pta_400[pta$regimen == reg & pta$uv %in% uvs]
  if (length(v) != length(uvs)) stop("no unique rows for regimen '", reg, "'")
  v
}

claims <- dplyr::bind_rows(
  claim(
    "Fig 3A: 1500 mg q12h attains AUC >= 400 at every urine volume (curve on the 100% ceiling)",
    sprintf("min %.1f%%", min(pta_a("1500 mg q12h"))),
    min(pta_a("1500 mg q12h")) >= 99
  ),
  claim(
    "Fig 3A / Discussion: 1000 mg q12h stays above the 90% line under various residual diuresis",
    sprintf("min %.1f%%", min(pta_a("1000 mg q12h"))),
    min(pta_a("1000 mg q12h")) >= 90
  ),
  claim(
    "Abstract / Discussion: 750 mg q12h is suitable for anuria or oliguria (UV <= 500 mL/24h)",
    sprintf("min %.1f%%", min(pta_a("750 mg q12h", c(0, 100, 500)))),
    min(pta_a("750 mg q12h", c(0, 100, 500))) >= 90
  ),
  claim(
    "Results: PTA decreases as daily urine volume increases",
    sprintf("750 mg q12h: %.1f%% at UV 0 vs %.1f%% at UV 3000",
            pta_a("750 mg q12h", 0), pta_a("750 mg q12h", 3000)),
    pta_a("750 mg q12h", 3000) < pta_a("750 mg q12h", 0)
  ),
  claim(
    "Discussion: no regimen has a satisfactory PTA for the AUC 400-600 target (Fig 3B peaks near 57%)",
    sprintf("max %.1f%%", max(pta$pta_400600)),
    max(pta$pta_400600) < 70
  ),
  claim(
    "Fig 3A: 500 mg regimens reproduce the published curve",
    sprintf("up to %.1f points above Fig 3A",
            max(cmp$diff[cmp$regimen %in% c("500 mg qd", "500 mg q12h")])),
    FALSE, deviation = TRUE
  )
)

knitr::kable(claims, caption = "Published claims checked against the packaged model.")
Published claims checked against the packaged model.
Claim Achieved Pass Deviation
Fig 3A: 1500 mg q12h attains AUC >= 400 at every urine volume (curve on the 100% ceiling) min 100.0% TRUE FALSE
Fig 3A / Discussion: 1000 mg q12h stays above the 90% line under various residual diuresis min 95.5% TRUE FALSE
Abstract / Discussion: 750 mg q12h is suitable for anuria or oliguria (UV <= 500 mL/24h) min 93.0% TRUE FALSE
Results: PTA decreases as daily urine volume increases 750 mg q12h: 100.0% at UV 0 vs 80.5% at UV 3000 TRUE FALSE
Discussion: no regimen has a satisfactory PTA for the AUC 400-600 target (Fig 3B peaks near 57%) max 42.5% TRUE FALSE
Fig 3A: 500 mg regimens reproduce the published curve up to 22.5 points above Fig 3A FALSE TRUE 
stopifnot(all(claims$Pass[!claims$Deviation]))

Assumptions and deviations

  • LOG in the Table 2 footnote is read as base-10. The paper does not define it. Base-10 makes the divisor 2.3 the cohort median-normalising constant log10(160 + 10) = 2.230, which is what Yu 2023 Methods describes (“continuous covariates were modeled using a median-normalized model”), and it is the only reading consistent with Figure 3A: the natural-log reading puts 1500 mg q12h at 22% attainment at 3000 mL/24h where the paper shows 100%. Worked through in “Which logarithm?” above.

  • omega and sigma in Table 2 are read as variances scaled by 100, giving omega = 0.348 and a 31.3% proportional residual error, rather than as standard deviations. Table 2’s own abbreviation footnote calls omega the “interindividual variance for CL”; Figure 2I’s ETA histogram and the slope of Figure 3A’s 500 mg q12h curve both agree. Worked through in “What scale are omega and sigma on?” above.

  • lcl is not a typical clearance. Because the covariate term is not normalised to 1 at the cohort median, exp(lcl) = 1.05 L/h is the coefficient of the power term, i.e. the clearance that would apply at UV + 10 = 1 mL. The typical anuric subject has CL = 1.39 L/h and the typical subject at the cohort median urine volume has CL = 1.96 L/h. The abstract’s “final typical clearance was 1.05 L/h” is quoting the theta, not a clearance any patient in the cohort has.

  • KNOWN DEVIATION: the 500 mg regimens over-attain relative to Figure 3A by 8 to 20 percentage points, while 750 / 1000 / 1500 mg q12h reproduce within 5.5 points. The pattern is what a uniform exposure offset looks like: the larger regimens sit on the saturated part of the PTA curve where a shift is invisible, and the 500 mg regimens sit on the steep part where the same shift is worth 10-20 points. Fitting a scale factor and a spread to the digitised 500 mg q12h curve reproduces the paper’s omega to within 1.5% but implies a typical clearance about 18% higher than the printed equation gives. No reading of the printed equation, and no offset other than +10 or divisor other than 2.3, produces an 18% uniform shift, so the discrepancy is between Table 2’s printed equation and the paper’s own Monte Carlo rather than in this transcription. The packaged model implements the printed equation; the Figure 3A comparison table above keeps the disagreement visible, and the claims gate excludes only this row.

  • Dosing is simulated as 1-hour intravenous infusions. The paper says only “intermittent intravenous vancomycin” and does not report an infusion duration. Steady-state AUC0-24 – the quantity every published check here depends on – is identical for any infusion duration, so this affects only the Cmax,ss column of the NCA summary.

  • Steady state is imposed with ss = 1 rather than reached by simulating a run-in. The slowest subject in the cohort (anuric, eta at the 0.25th percentile) has a half-life near 70 h, so a run-in would need roughly three weeks of dosing to be honest; ss = 1 is exact.

  • The cohort is deterministic, not random. etalcl is set to the 200 equally spaced quantiles of N(0, 0.121) and supplied per subject, so no RNG is used anywhere in this vignette and every number above is identical on any machine and any solver-thread count.

  • Covariate distributions are not simulated. The paper’s own Monte Carlo fixes the urine volume at six values rather than sampling it, and no other covariate enters the model, so the cohort varies only in etalcl and in the fixed urine-volume scenario – exactly as in Yu 2023 Figure 3.

  • Figure-derived values. The fig3a_published and fig3a_500q12h tables were digitised by the operator from Yu 2023 Figure 3A (roughly +/-3 percentage points). They are used only as comparison references and for the omega slope check; no parameter value in inst/modeldb/specificDrugs/Yu_2023_vancomycin.R comes from a figure.

  • Screened but unretained covariates (age, sex, weight, BMI, daily dose, serum creatinine, blood urea nitrogen, CRRT modality) are recorded in the model file’s covariatesDataExcluded metadata with their Supplementary Table S1 delta-OFVs, so the paper’s covariate screen is preserved without carrying unused covariateData entries.

  • No published NCA table to compare against. Yu 2023 reports no Cmax / Tmax / AUC / half-life summary for the observed data, so there is nothing for nlmixr2lib::ncaComparisonTable() to consume. The PKNCA output is instead validated against the closed-form steady-state identity and used to reproduce the published PTA curves.