Skip to contents

Tigecycline (Su 2024)

Source: vignettes/articles/Su_2024_tigecycline.Rmd
Su_2024_tigecycline.Rmd

Model and source

Population

Su and colleagues ran a single-centre prospective study at Sir Run Run Shaw Hospital (Hangzhou, China) between July 2019 and July 2023, enrolling 98 critically ill adults who received intermittent intravenous tigecycline as a 30-minute infusion every 12 h. Sampling was intensive rather than opportunistic: eight samples per patient, drawn immediately before the seventh dose and at 0.5, 1, 2, 3, 4, 6 and 12 h afterwards, so every subject contributed a complete steady-state dosing interval. That yielded 751 concentrations for model building.

The cohort was 66.3% male with a mean age of 63.4 years (SD 18.0, range 20-92) and a mean weight of 62.3 kg (SD 12.4, range 38-92.5). Renal function spanned almost two orders of magnitude - Cockcroft-Gault creatinine clearance median 77.0 mL/min, IQR 51.1-142, range 6.40-338 - while patients on haemodialysis, peritoneal dialysis or CRRT were excluded. Liver chemistry was likewise dispersed: albumin median 27.9 g/L (IQR 25.8-30.3), GGT median 53 U/L (IQR 30-103) and total bilirubin median 15.7 umol/L (IQR 10.0-30.4). Infections were predominantly intra-abdominal (37.8%) and pulmonary (31.6%), most often caused by carbapenem-resistant A. baumannii (46.9%) or K. pneumoniae (24.5%). Almost all patients (92.9%) received 100 mg/day; 7.14% received 200 mg/day. Baseline demographics are Table 1 of the paper.

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

Source trace

Every ini() entry carries an in-file comment pointing at its source location in inst/modeldb/specificDrugs/Su_2024_tigecycline.R. The table below collects them for review.

Equation / parameter Value Source location
Two-compartment structure, first-order elimination n/a Results, “PPK model development” (AIC 8242.914 two-compartment vs 9325.681 one-compartment)
30-minute IV infusion, q12h n/a Methods, “Patient inclusion, drug administration and sample collection”
lcl (CL intercept) 3.09 L/h Table 2, “CL (L/h)”; final-model equation CL (L/h) = 3.09+(CCr/77) x 3.28
e_crcl_cl 3.28 L/h per (CCr/77) Table 2, theta CLCR-CL
lvc (V1 at 61 kg) 32.1 L Table 2, “V1 (L)”; final-model equation V1 (L) = 32.1x(BW/61)^1.95
e_wt_vc 1.95 Table 2, theta BW-V1
lq (Q at reference liver chemistry) 39.7 L/h Table 2, “Q (L/h)”; final-model equation Q (L/h) = 39.7x(logGGT/1.7)^0.956 x (logTBIL/1.2)^-0.912
e_ggt_q 0.956 Table 2, theta GGT-Q
e_tbili_q -0.912 Table 2, theta TBIL-Q
lvp (V2 at 61 kg, reference albumin) 113 L Table 2, “V2 (L)”; final-model equation V2 (L) = 113x(BW/61)^1.61 x (logALB/1.4)^4.52
e_wt_vp 1.61 Table 2, theta BW-V2
e_alb_vp 4.52 Table 2, theta ALB-V2
Linear covariate form theta_pop + theta_cov * cov/cov_median n/a Methods Eq. 1 (used for CCr on CL)
Power covariate form theta_pop * (cov/cov_median)^theta_cov n/a Methods Eq. 3 (used for the four remaining effects)
etalcl 27.0% CV Table 2, omega CL (%)
etalvc 72.5% CV Table 2, omega V1 (%)
etalq 18.8% CV Table 2, omega Q (%)
No IIV on V2 n/a Results, “it was not possible to estimate the interindividual variability for V2”
propSd 0.142 Table 2, sigma (%) 2.02, read as a variance (see Assumptions)
Covariate reference values 77 mL/min, 61 kg, 1.7, 1.2, 1.4 n/a Divisors printed inside the four final-model equations
PK/PD targets AUC/MIC >= 6.96 / 12.8 / 17.9 n/a Methods, “Simulation and dosing regimen optimization”
PTA strata CCr 30 / 80 / 130 mL/min n/a Methods, “Simulation and dosing regimen optimization”

Virtual cohort

The original concentrations are not public, so the checks below use virtual cohorts whose covariate distributions reproduce the Table 1 marginals. Skewed laboratory covariates are drawn log-normally with the log-scale SD implied by the published inter-quartile range (sd = log(IQR ratio) / 1.349); body weight is drawn normally from its published mean and SD. Every draw is truncated to the observed range in Table 1.

Five cohorts of 200 subjects each are built. Four are dosed to exact steady state (rxode2’s ss = 1): one with the full covariate distribution and three with creatinine clearance pinned at 30, 80 and 130 mL/min (the strata of Figure 4). These carry the NCA and PTA work, where steady state is the intended condition. The fifth replicates the paper’s actual sampling occasion - the seventh q12h dose, i.e. 72 h of therapy, which is not steady state for the slower subjects - and is used only for the Figure 3 comparison, so that simulated and observed profiles are compared at the same accumulation state.

# set.seed() seeds R's RNG. It does NOT seed rxode2's simulation RNG, and
# rxode2's streams are partitioned PER SOLVER THREAD, so the cohort below is
# reproducible on this machine and different on a machine with a different
# thread count. Every assertion downstream is written to hold for ANY cohort
# the model can produce; see pattern 12 of the skill's
# known-vignette-failure-patterns reference.
set.seed(20240129)
rxode2::rxSetSeed(20240129)

n_per_arm <- 200L

# Table 1 marginals. sd_log is recovered from the published IQR.
sd_log <- function(q1, q3) log(q3 / q1) / 1.349

draw_covariates <- function(n, crcl = NULL) {
  tibble::tibble(
    CRCL  = if (is.null(crcl)) {
      pmax(6.4, pmin(338, exp(rnorm(n, log(77), sd_log(51.1, 142)))))
    } else {
      rep(crcl, n)
    },
    WT    = pmax(38, pmin(92.5, rnorm(n, 62.3, 12.4))),
    ALB   = pmax(17.3, pmin(38.3, exp(rnorm(n, log(27.9), sd_log(25.8, 30.3))))),
    GGT   = pmax(8,    pmin(1088, exp(rnorm(n, log(53),   sd_log(30, 103))))),
    TBILI = pmax(4,    pmin(144,  exp(rnorm(n, log(15.7), sd_log(10.0, 30.4)))))
  )
}

# Steady state is imposed EXACTLY, with rxode2's `ss = 1` dosing flag, rather
# than approached by stacking a finite number of doses. That matters here: the
# peripheral volume is large (113 L at the reference covariates) and IIV on V1
# is 72.5%, so a subject in the tail of the eta distribution has a terminal
# half-life of several days and is still ~50% short of steady state after nine
# q12h doses. Accumulation shortfall is a PHYSICAL, per-subject difference, so
# it would silently corrupt the AUC identity below and vary with whichever
# cohort the machine happens to draw. `ss = 1` removes it entirely.
tad_obs <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 10, 12)

# `n_doses = NULL` gives exact steady state; an integer instead stacks that many
# q12h doses and samples around the last one, which is what the paper did.
make_cohort <- function(n, stratum, crcl = NULL, id_offset = 0L, n_doses = NULL) {
  cov <- draw_covariates(n, crcl) |>
    dplyr::mutate(id = id_offset + seq_len(n), stratum = stratum)
  t_ref <- if (is.null(n_doses)) 0 else (n_doses - 1) * 12
  doses <- cov |>
    dplyr::mutate(time = 0, amt = 50, evid = 1L, dur = 0.5, ii = 12,
                  addl = if (is.null(n_doses)) NA_integer_ else n_doses - 1L,
                  ss   = if (is.null(n_doses)) 1L else 0L,
                  cmt  = "central")
  obs <- cov |>
    tidyr::crossing(time = t_ref + tad_obs) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, dur = NA_real_,
                  ii = NA_real_, addl = NA_integer_, ss = NA_integer_,
                  cmt = "central")
  dplyr::bind_rows(doses, obs) |>
    dplyr::mutate(t_ref = t_ref) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

ss_strata   <- c("Full covariate distribution", "CCr 30 mL/min",
                 "CCr 80 mL/min", "CCr 130 mL/min")
crcl_strata <- c("CCr 30 mL/min", "CCr 80 mL/min", "CCr 130 mL/min")
dose7_stratum <- "Seventh dose (paper sampling)"

events <- dplyr::bind_rows(
  make_cohort(n_per_arm, "Full covariate distribution", NULL, id_offset =   0L),
  make_cohort(n_per_arm, "CCr 30 mL/min",                30, id_offset = 200L),
  make_cohort(n_per_arm, "CCr 80 mL/min",                80, id_offset = 400L),
  make_cohort(n_per_arm, "CCr 130 mL/min",              130, id_offset = 600L),
  # Replicates the paper's own sampling occasion: the seventh q12h dose, i.e.
  # 72 h of therapy, NOT steady state. Used only for the Figure 3 comparison so
  # that simulated and observed profiles carry the same accumulation state.
  make_cohort(n_per_arm, dose7_stratum, NULL, id_offset = 800L, n_doses = 7L)
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

Simulation 

mod <- readModelDb("Su_2024_tigecycline")

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("stratum", "WT", "CRCL", "ALB", "GGT", "TBILI", "t_ref")
) |>
  as.data.frame() |>
  dplyr::mutate(tad = time - t_ref,
                Cc_ugL = Cc * 1000)   # model works in mg/L; the paper plots ug/L
#> ℹ parameter labels from comments will be replaced by 'label()'

# At steady state the trough must repeat: C(0) == C(12) for every subject.
ss_check <- sim |>
  dplyr::filter(stratum %in% ss_strata, tad %in% c(0, 12)) |>
  tidyr::pivot_wider(id_cols = c(stratum, id), names_from = tad,
                     values_from = Cc, names_prefix = "t") |>
  dplyr::mutate(pct = 100 * (t12 - t0) / t0)
stopifnot(max(abs(ss_check$pct)) < 1e-6)

Structural identity check

Before any cohort-level comparison, confirm that the typical-value parameters at the paper’s reference covariates reproduce the Table 2 estimates exactly. These are deterministic, so the tolerance is tight.

ref <- tibble::tibble(CRCL = 77, WT = 61, ALB = 10^1.4, GGT = 10^1.7, TBILI = 10^1.2)

typ <- rxode2::rxSolve(
  rxode2::zeroRe(mod),
  events = tibble::tibble(id = 1L, time = 0, amt = 50, evid = 1L, dur = 0.5,
                          cmt = "central") |>
    dplyr::bind_rows(tibble::tibble(id = 1L, time = c(0, 1), amt = NA_real_,
                                    evid = 0L, dur = NA_real_, cmt = "central")) |>
    dplyr::bind_cols(ref[rep(1, 3), ]) |>
    dplyr::arrange(time, dplyr::desc(evid))
) |> as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq'

structural <- tibble::tibble(
  Parameter = c("CL at CCr 77 (L/h)", "V1 at 61 kg (L)",
                "Q at reference liver chemistry (L/h)", "V2 at 61 kg, reference ALB (L)"),
  Published = c(3.09 + 3.28, 32.1, 39.7, 113),
  Model     = c(typ$cl[1], typ$vc[1], typ$q[1], typ$vp[1])
) |>
  dplyr::mutate(`% difference` = 100 * (Model - Published) / Published)

knitr::kable(structural, digits = 3,
             caption = "Typical-value parameters at the paper's reference covariates.")
Typical-value parameters at the paper’s reference covariates.
Parameter Published Model % difference
CL at CCr 77 (L/h) 6.37 6.37 0
V1 at 61 kg (L) 32.10 32.10 0
Q at reference liver chemistry (L/h) 39.70 39.70 0
V2 at 61 kg, reference ALB (L) 113.00 113.00 0 

# Deterministic - a mis-transcribed value or exponent breaks this immediately.
stopifnot(max(abs(structural$`% difference`)) < 0.5)

Note the first row: the paper’s Abstract and Results call 3.09 L/h “the typical value of CL”, but the printed final-model equation is CL = 3.09 + (CCr/77) x 3.28, so 3.09 L/h is the intercept and the typical patient at the median CCr of 77 mL/min clears 6.37 L/h. The equation governs; see Assumptions below.

Replicate published figures

Figure 3 - steady-state concentration-time profile

Figure 3 of Su 2024 is a visual predictive check over one dosing interval, with the observed median and 5th / 95th percentiles drawn on a log concentration axis in ug/L. Because the underlying observations are not public, the published percentiles below were read off the rendered figure and are therefore accurate to roughly +/-10%. They are used only for a wide order-of-magnitude gate (a factor-of-two error in dose or volume breaks it); nothing is tuned to them.

This comparison uses the seventh-dose cohort, not the exact-steady-state one, because Su 2024 sampled after 72 h of therapy and with 72.5% IIV on V1 a sizeable minority of subjects have not accumulated fully by then. Comparing an exact-steady-state simulation against those observations would bias the model side upwards.

vpc <- sim |>
  dplyr::filter(stratum == dose7_stratum) |>
  dplyr::group_by(tad) |>
  dplyr::summarise(Q05 = quantile(Cc_ugL, 0.05),
                   Q50 = quantile(Cc_ugL, 0.50),
                   Q95 = quantile(Cc_ugL, 0.95), .groups = "drop")

ggplot(vpc, aes(tad, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "firebrick", linewidth = 1) +
  scale_y_log10() +
  labs(x = "Time after dose (h)", y = "Tigecycline concentration (ug/L)",
       title = "Steady-state dosing interval, 50 mg q12h",
       caption = "Replicates Figure 3 of Su 2024 (median with 5th-95th percentile band).")

digitised <- tibble::tribble(
  ~tad, ~pub_Q05, ~pub_Q50, ~pub_Q95,
  0.5,       680,     1600,     3300,
  2,         255,      570,     1400,
  4,         210,      440,     1180,
  6,         195,      420,     1120,
  12,        100,      320,      750
)

fig3 <- digitised |>
  dplyr::left_join(vpc, by = "tad") |>
  dplyr::mutate(`Median ratio (model / paper)` = Q50 / pub_Q50,
                `Spread sd(log C), model` = log(Q95 / Q05) / 3.29,
                `Spread sd(log C), paper` = log(pub_Q95 / pub_Q05) / 3.29)

fig3 |>
  dplyr::select(tad, pub_Q50, Q50, `Median ratio (model / paper)`,
                `Spread sd(log C), paper`, `Spread sd(log C), model`) |>
  dplyr::rename("TAD (h)" = tad,
                "Paper median (ug/L)" = pub_Q50,
                "Model median (ug/L)" = Q50) |>
  knitr::kable(digits = c(1, 0, 0, 2, 3, 3),
               caption = "Digitised Figure 3 percentiles vs the simulated cohort. Published values were read off the rendered figure (+/-10%).")
Digitised Figure 3 percentiles vs the simulated cohort. Published values were read off the rendered figure (+/-10%).
TAD (h) Paper median (ug/L) Model median (ug/L) Median ratio (model / paper) Spread sd(log C), paper Spread sd(log C), model
0.5 1600 1425 0.89 0.480 0.388
2.0 570 656 1.15 0.518 0.337
4.0 440 562 1.28 0.525 0.372
6.0 420 506 1.21 0.531 0.400
12.0 320 407 1.27 0.612 0.479

The medians agree to within about a third across the interval, and the model’s dispersion is of the right order. Two systematic differences are visible and are recorded as deviations rather than tuned away. First, the simulated profile is flatter than the observed one: the model matches the peak (ratio 0.87) but over-predicts the 12 h trough by roughly 35%. Second, the simulated 5th-95th spread is narrower than the observed spread at every time point - sd(log C) of roughly 0.35-0.52 against 0.48-0.61 - with the gap largest early in the interval. The second of these bears directly on the Table 2 variability block, whose reported scale is ambiguous; note that the wider (variance) reading of the omega rows would close this particular gap while over-shooting elsewhere. See Assumptions.

# Loose, but still red-capable: a mis-transcribed volume (32.1 vs 113 L) or a
# dose error (50 vs 100 mg) moves these medians by a factor of 2 or more.
stopifnot(all(fig3$`Median ratio (model / paper)` > 0.6),
          all(fig3$`Median ratio (model / paper)` < 1.7))

Figure 4 - probability of target attainment

Figure 4 reports PTA against the three CART-derived PK/PD targets (AUC/MIC >= 6.96 for complicated intra-abdominal infection, >= 12.8 for community-acquired pneumonia and >= 17.9 for complicated skin and skin structure infection), stratified by creatinine clearance, for the 50 mg q12h and 100 mg q12h regimens.

Steady-state exposure in a linear model is exactly dose-proportional, so the 100 mg q12h arm is obtained by doubling the simulated 50 mg AUC rather than by simulating a second cohort.

auc_ss <- sim |>
  dplyr::filter(stratum %in% crcl_strata) |>
  dplyr::group_by(stratum, id) |>
  dplyr::arrange(tad, .by_group = TRUE) |>
  dplyr::summarise(
    # Linear-up / linear-down trapezoid over the steady-state interval, then
    # scaled to 24 h; the paper's targets are AUC0-24 / MIC ratios.
    auc24_50mg = 2 * sum(diff(tad) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups = "drop"
  ) |>
  dplyr::mutate(auc24_100mg = 2 * auc24_50mg)

mic_grid <- c(0.125, 0.25, 0.5, 1, 2, 4, 8)
targets  <- c(`AUC/MIC >= 6.96 (cIAI)` = 6.96,
              `AUC/MIC >= 12.8 (CAP)`  = 12.8,
              `AUC/MIC >= 17.9 (cSSSI)` = 17.9)

pta <- tidyr::expand_grid(
  stratum = unique(auc_ss$stratum),
  regimen = c("50 mg q12h", "100 mg q12h"),
  MIC     = mic_grid,
  target_label = names(targets)
) |>
  dplyr::rowwise() |>
  dplyr::mutate(
    PTA = {
      a <- auc_ss$auc24_50mg[auc_ss$stratum == stratum]
      if (regimen == "100 mg q12h") a <- 2 * a
      100 * mean(a / MIC >= targets[[target_label]])
    }
  ) |>
  dplyr::ungroup() |>
  dplyr::mutate(stratum = factor(stratum, levels = c("CCr 30 mL/min", "CCr 80 mL/min",
                                                     "CCr 130 mL/min")))

ggplot(pta, aes(MIC, PTA, colour = stratum, linetype = regimen)) +
  geom_line() + geom_point(size = 1) +
  geom_hline(yintercept = 90, linetype = "dotted") +
  scale_x_log10(breaks = mic_grid) +
  facet_wrap(~target_label) +
  labs(x = "MIC (mg/L)", y = "PTA (%)", colour = NULL, linetype = NULL,
       caption = "Replicates Figure 4 of Su 2024 (panels A-C, stratified by creatinine clearance).") +
  theme(legend.position = "bottom")

# Deterministic consequence of the CL equation: exposure falls as CCr rises.
med_auc <- auc_ss |>
  dplyr::group_by(stratum) |>
  dplyr::summarise(median_auc24 = median(auc24_50mg), .groups = "drop")

knitr::kable(med_auc, digits = 2,
             caption = "Median steady-state AUC0-24 (mg*h/L) at 50 mg q12h by CCr stratum.")
Median steady-state AUC0-24 (mg*h/L) at 50 mg q12h by CCr stratum.
stratum median_auc24
CCr 130 mL/min 11.68
CCr 30 mL/min 22.24
CCr 80 mL/min 15.28 

ratio_30_130 <- med_auc$median_auc24[med_auc$stratum == "CCr 30 mL/min"] /
  med_auc$median_auc24[med_auc$stratum == "CCr 130 mL/min"]
# CL(30)/CL(130) = (3.09 + 3.28*30/77) / (3.09 + 3.28*130/77) = 4.368 / 8.627,
# so the exposure ratio is ~1.98 by construction and only cohort noise moves it.
stopifnot(ratio_30_130 > 1.7, ratio_30_130 < 2.3)

PKNCA validation

PKNCA is used for the steady-state interval NCA. The strongest available gate here is an internal identity: for a linear model at steady state, the AUC over one dosing interval must equal dose / CL exactly. The paper reports no NCA table of its own, so this identity - plus the structural check above - is what pins the transcription.

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

conc_obj <- PKNCA::PKNCAconc(as.data.frame(sim_nca), Cc ~ time | stratum + id)

dose_df <- events |>
  dplyr::filter(evid == 1, stratum %in% ss_strata) |>
  dplyr::select(id, stratum, time, amt)

dose_obj <- PKNCA::PKNCAdose(as.data.frame(dose_df), amt ~ time | stratum + id)

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

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

nca_summary <- nca_tab |>
  dplyr::filter(PPTESTCD %in% c("auclast", "cmax", "tmax", "cmin")) |>
  tidyr::pivot_wider(id_cols = c(stratum, id), names_from = PPTESTCD,
                     values_from = PPORRES)

nca_summary |>
  dplyr::group_by(stratum) |>
  dplyr::summarise(dplyr::across(c(auclast, cmax, tmax, cmin), median), .groups = "drop") |>
  dplyr::rename("Stratum" = stratum,
                "AUC0-12,ss (mg*h/L)" = auclast,
                "Cmax,ss (mg/L)" = cmax,
                "Tmax (h)" = tmax,
                "Cmin,ss (mg/L)" = cmin) |>
  knitr::kable(digits = 3,
               caption = "Median steady-state NCA parameters at 50 mg q12h, by CCr stratum.")
Median steady-state NCA parameters at 50 mg q12h, by CCr stratum.
Stratum AUC0-12,ss (mg*h/L) Cmax,ss (mg/L) Tmax (h) Cmin,ss (mg/L)
CCr 130 mL/min 5.838 1.392 0.5 0.317
CCr 30 mL/min 11.117 1.921 0.5 0.754
CCr 80 mL/min 7.634 1.561 0.5 0.456
Full covariate distribution 7.742 1.643 0.5 0.448 
# AUC0-12,ss must equal dose / CL for every individual. Because steady state is
# imposed exactly (`ss = 1`), the only residual difference is trapezoidal
# discretisation of a convex curve, which is one-sided: every subject should
# come out slightly HIGH, by well under a percent.
cl_by_id <- sim |>
  dplyr::filter(stratum %in% crcl_strata) |>
  dplyr::group_by(stratum, id) |>
  dplyr::summarise(cl = dplyr::first(cl), .groups = "drop")

identity_chk <- nca_summary |>
  dplyr::inner_join(cl_by_id, by = c("stratum", "id")) |>
  dplyr::mutate(expected = 50 / cl,
                pct_diff = 100 * (auclast - expected) / expected)

knitr::kable(
  identity_chk |>
    dplyr::group_by(stratum) |>
    dplyr::summarise(`Median % difference` = median(pct_diff),
                     `Min % difference`    = min(pct_diff),
                     `Max % difference`    = max(pct_diff), .groups = "drop") |>
    dplyr::rename("Stratum" = stratum),
  digits = 3,
  caption = "PKNCA AUC0-12,ss against the closed-form dose / CL. Pure trapezoidal error."
)
PKNCA AUC0-12,ss against the closed-form dose / CL. Pure trapezoidal error.
Stratum Median % difference Min % difference Max % difference
CCr 130 mL/min 0.235 -0.014 0.525
CCr 30 mL/min 0.140 -0.087 0.411
CCr 80 mL/min 0.191 -0.031 0.511 

# Both sides use the same drawn parameters and steady state is exact, so this is
# numerical error only and a tight two-sided bound is correct (see the repository
# note on cohort-extreme assertions). It is not quite one-sided: the trapezoid
# over-estimates the convex post-peak decline but under-estimates the concave
# rise during the 0.5 h infusion, and which term wins depends on the subject.
stopifnot(max(abs(identity_chk$pct_diff)) < 2)

Published dosing claims

Su 2024 states three quantitative conclusions from the Monte Carlo simulation. Each is checked below against the packaged model at the median-CCr stratum, which is the closest match to the paper’s typical patient.

pta_at <- function(stratum_, regimen_, mic, target_label_) {
  v <- pta$PTA[pta$stratum == stratum_ & pta$regimen == regimen_ &
                 pta$MIC == mic & pta$target_label == target_label_]
  if (length(v) != 1L) stop("no unique PTA row for ", stratum_, " / ", regimen_,
                            " / MIC ", mic, " / ", target_label_)
  v
}

claims <- tibble::tribble(
  ~Claim, ~Source, ~Achieved, ~Pass, ~Deviation,
  "50 mg q12h attains the cIAI target (6.96) with PTA >= 90% at MIC 1",
  "Results, 'sufficient to achieve an AUC/MIC ratio of 6.96 at MICs <= 1 mg/L'",
  pta_at("CCr 80 mL/min", "50 mg q12h", 1, "AUC/MIC >= 6.96 (cIAI)"), NA, FALSE,

  "50 mg q12h does NOT attain the cIAI target at MIC 2",
  "Results, same sentence (attainment limited to MIC <= 1 mg/L)",
  pta_at("CCr 80 mL/min", "50 mg q12h", 2, "AUC/MIC >= 6.96 (cIAI)"), NA, FALSE,

  "50 mg q12h does NOT attain the CAP target (12.8) at MIC 1",
  "Results, 'may attain the suboptimal target at MICs > 0.5' for AUC/MIC 12.8",
  pta_at("CCr 80 mL/min", "50 mg q12h", 1, "AUC/MIC >= 12.8 (CAP)"), NA, FALSE,

  "100 mg q12h restores the CAP target at MIC 1 (PTA >= 90%)",
  "Conclusion, '100 mg every 12 h was needed for community-acquired pneumonia'",
  pta_at("CCr 80 mL/min", "100 mg q12h", 1, "AUC/MIC >= 12.8 (CAP)"), NA, FALSE,

  "50 mg q12h does NOT attain the cSSSI target (17.9) at MIC 0.5",
  "Results, 'may attain the suboptimal target at ... 0.25 mg/L' for AUC/MIC 17.9",
  pta_at("CCr 80 mL/min", "50 mg q12h", 0.5, "AUC/MIC >= 17.9 (cSSSI)"), NA, TRUE
)

claims$Pass <- c(
  claims$Achieved[1] >= 90,
  claims$Achieved[2] <  90,
  claims$Achieved[3] <  90,
  claims$Achieved[4] >= 90,
  claims$Achieved[5] <  90
)

claims |>
  dplyr::mutate(Achieved = sprintf("%.1f%% PTA", Achieved)) |>
  dplyr::rename("Achieved (model)" = Achieved, "Known deviation" = Deviation) |>
  knitr::kable(caption = "Published dosing conclusions checked against the packaged model.")
Published dosing conclusions checked against the packaged model.
Claim Source Achieved (model) Pass Known deviation
50 mg q12h attains the cIAI target (6.96) with PTA >= 90% at MIC 1 Results, ‘sufficient to achieve an AUC/MIC ratio of 6.96 at MICs <= 1 mg/L’ 100.0% PTA TRUE FALSE
50 mg q12h does NOT attain the cIAI target at MIC 2 Results, same sentence (attainment limited to MIC <= 1 mg/L) 66.0% PTA TRUE FALSE
50 mg q12h does NOT attain the CAP target (12.8) at MIC 1 Results, ‘may attain the suboptimal target at MICs > 0.5’ for AUC/MIC 12.8 76.0% PTA TRUE FALSE
100 mg q12h restores the CAP target at MIC 1 (PTA >= 90%) Conclusion, ‘100 mg every 12 h was needed for community-acquired pneumonia’ 100.0% PTA TRUE FALSE
50 mg q12h does NOT attain the cSSSI target (17.9) at MIC 0.5 Results, ‘may attain the suboptimal target at … 0.25 mg/L’ for AUC/MIC 17.9 98.5% PTA FALSE TRUE 

# Gate on the claims that reproduce; the flagged deviation stays visible in the
# table but is excluded rather than having the bound widened until it passes.
stopifnot(all(claims$Pass[!claims$Deviation]))

The one deviation is the cSSSI boundary. Su 2024 places the 50 mg q12h attainment ceiling for the 17.9 target at MIC 0.25 mg/L; the packaged model puts it at 0.5 mg/L, i.e. one two-fold dilution more permissive. The direction is consistent with the discrepancies noted for Figure 3: the transcribed model carries less exposure variability than the published analysis did, which lifts PTA in the tail of the MIC grid. The paper’s headline recommendation - 50 mg q12h for complicated intra-abdominal infection, 100 mg q12h for pneumonia and skin and skin structure infections - is reproduced.

Assumptions and deviations

  • CL = 3.09 L/h is an intercept, not a typical value. The Abstract, the Results and the Conclusion all describe 3.09 L/h as “the typical value of CL”, but Methods Eq. 1 defines the linear covariate model as theta_i = theta_pop + theta_cov * cov_i / cov_median and the printed final-model equation is CL = 3.09 + (CCr/77) x 3.28. Evaluated at the cohort’s median CCr of 77 mL/min this gives 6.37 L/h, and the equation is encoded as printed. Three independent lines of evidence support the equation over the prose. First, the other three parameters are internally consistent with their equations (V1, Q and V2 all reduce to their Table 2 values at the reference covariates), so CL is the only row where “Table 2 value” and “typical value” would differ - exactly what an intercept does. Second, Supplementary Table S1 tabulates ten previously published tigecycline models with CL between 4.81 and 23.1 L/h; 6.37 L/h sits at the low end, consistent with the paper’s own claim that its CL was “lower than those estimated in previously published studies”, whereas 3.09 L/h would be far below every one of them. Third, the paper’s PTA conclusions (Figure 4) are attainable only with the higher clearance: at 3.09 L/h the 50 mg q12h regimen would deliver an AUC0-24 of 32 mg*h/L and would comfortably clear the 6.96 target at MIC 2, contradicting the stated MIC <= 1 mg/L ceiling. The same reading - Table 2 reporting the intercept of an additive linear renal-clearance model - is used in Delattre_2010_amikacin, which has the identical covariate structure.
  • log in the Q and V2 equations is base 10. The paper writes logGGT, logTBIL and logALB without stating the base. The divisors printed in the equations settle it: 1.7, 1.2 and 1.4 against the Table 1 medians of 53 U/L, 15.7 umol/L and 27.9 g/L give log10 values of 1.72, 1.20 and 1.45, while the natural logarithms are 3.97, 2.75 and 3.33. All three agree on base 10, so it is not a coincidence of rounding on any single covariate. The normalizer is the median of the log-transformed covariate, not the log of the median covariate - the equations divide logGGT by 1.7, not by log10(53) written out - but for these values the two are the same number to the printed precision.
  • The residual-error scale in Table 2 is ambiguous, and the literal reading is falsified. The variability block prints omega CL (%) 27.0, omega V1 (%) 72.5, omega Q (%) 18.8 and sigma (%) 2.02. Read uniformly as coefficients of variation, the first three are ordinary popPK values but the last is a 2% proportional residual error, which is not attainable for clinical LC-MS/MS concentration data. The paper’s own Figure 2 rules it out directly: a 2.02% proportional CV puts the 95% band of DV about IPRED at +/-4%, whereas digitising the DV-vs-IPRED panel gives a band of roughly +/-20-35% in the densely sampled region - an order of magnitude wider, and far beyond the plotting-symbol size. The value is therefore encoded as the proportional-error variance expressed in percent, sigma^2 = 0.0202, i.e. propSd = sqrt(0.0202) = 0.142, which reproduces the observed scatter. The three omega rows are encoded as coefficients of variation via the repository convention omega^2 = log(CV^2 + 1), matching how the same style of table is handled in Delattre_2010_amikacin. This mixed reading is the best-supported one, but it is a reading: applying the variance convention to the omega rows as well would give SDs of 0.52 / 0.85 / 0.43 instead of 0.27 / 0.72 / 0.19. Digitising the Figure 3 VPC does not cleanly separate the two - the coefficient-of-variation reading is slightly too narrow at late times and the variance reading is clearly too wide - and digitising the Figure 4 PTA slopes gives an intermediate value. Any downstream use that depends on the magnitude of the between-subject variability (rather than on typical-value exposure) should treat this as uncertain.
  • Steady state is imposed exactly, not approached. Su 2024 sampled around the seventh dose, i.e. after 72 h of q12h therapy. Reproducing that literally by stacking doses is not reproducible across cohorts: IIV on V1 is 72.5% and V2 is large, so a subject in the tail of the eta distribution has a terminal half-life of several days and is still tens of percent short of steady state at 72 h - or even at 480 h. The vignette therefore uses rxode2’s ss = 1 dosing flag, which solves the exact periodic steady state for each subject’s own parameters. This makes the AUC identity a statement about numerical quadrature rather than about accumulation, which is what it is meant to test, and it matches the steady-state condition the paper’s own Monte Carlo simulation assumes. The one comparison against published observations - Figure 3 - deliberately does not use it, and instead replicates the seventh-dose occasion the paper actually sampled, so that both sides carry the same accumulation state.
  • No IIV on V2. The paper states plainly that “it was not possible to estimate the interindividual variability for V2”, so V2 carries no eta. This is a reported structural fact, not an omission.
  • Covariate reference values. 77 mL/min, 61 kg, 1.7, 1.2 and 1.4 are taken from the divisors inside the printed final-model equations. Table 1 reports body weight as a mean (62.3 kg) rather than a median, so 61 kg is known only from the equation; the two are close enough that the distinction does not matter numerically.
  • Screened but not retained. Supplementary Table S3 records that daily dose on V1 and total bilirubin on CL both entered forward inclusion and were then removed in backward elimination, and that age, BMI, sex, BUN, ALT, AST and ALP were screened without entering. Only the five covariates of the final model are encoded.
  • Virtual cohort. Covariate distributions are reconstructed from the Table 1 marginals (mean and SD for weight, median and IQR for the laboratory values) under independence. The real cohort’s covariates are certainly correlated - albumin, bilirubin and GGT co-vary with illness severity - so the simulated spread of V2 and Q is wider than the study’s would have been.
  • Concentration units. The model works in mg/L (mg dose over L volume); Su 2024 plots concentrations in ug/L, i.e. 1000 x the model’s Cc. The vignette carries an explicit Cc_ugL column for figure comparisons.
  • No published NCA table. Su 2024 reports no observed Cmax / AUC / half-life summary, so there is no external NCA comparison to make. The PKNCA section gates on the internal steady-state identity AUC0-12,ss = dose / CL instead.