Skip to contents

Model and source

Dimelow 2018 developed two independent population PK models – one per drug – each pairing a three-compartment intravenous plasma model with a direct response (instantaneous-equilibrium) link to the epithelial lining fluid (ELF) of the lung. Because the two drugs were modelled separately, they are packaged as two model files sharing this one vignette.

cef <- rxode2::rxode(readModelDb("Dimelow_2018_ceftazidime"))
#> ℹ parameter labels from comments will be replaced by 'label()'
avi <- rxode2::rxode(readModelDb("Dimelow_2018_avibactam"))
#> ℹ parameter labels from comments will be replaced by 'label()'
  • Citation: Dimelow R, Wright JG, MacPherson M, Newell P, Das S. Population pharmacokinetic modelling of ceftazidime and avibactam in the plasma and epithelial lining fluid of healthy volunteers. Drugs R D. 2018;18(3):221-230. doi:10.1007/s40268-018-0241-0 (Sects. 2.2, 2.2.1, 3.1, 3.2; Tables 1 and 2). Subject data are from the phase I open-label ELF study NCT01395420.

  • Article: https://doi.org/10.1007/s40268-018-0241-0

  • Ceftazidime (Dimelow_2018_ceftazidime): Three-compartment intravenous population PK model for ceftazidime in the plasma of healthy male volunteers, jointly fitted with a direct-response (instantaneous-equilibrium) saturable plasma-to- epithelial-lining-fluid (ELF) link. Plasma disposition is parameterized by CL, V1, V2, V3, Q12 and Q13 with log-normal between-subject variability on every structural parameter. The ELF observable Celf is algebraic, not a compartment: an effect-site equilibration model was fitted and rejected (ELF half-life 13 min, no OFV improvement), so ELF concentration is an instantaneous Michaelis-Menten function of the total plasma concentration, Celf = emax_elf * Cc / (km_elf + Cc). The ELF:plasma penetration ratio therefore falls as plasma concentration rises: emax_elf / km_elf = 63.3% as Cc approaches zero, 52% at the efficacy-relevant Cc of 15.3 mg/L, and 32.1% at a typical Cmax of 70 mg/L. Ceftazidime plasma proportional residual noise itself carries between-subject variability (median 0.101, 54% CV); the ELF residual could not be identified from one ELF sample per subject and was FIXED to the plasma median. Penetration ratios are relative to TOTAL plasma and TOTAL ELF concentrations.

  • Avibactam (Dimelow_2018_avibactam): Three-compartment intravenous population PK model for avibactam in the plasma of healthy male volunteers, jointly fitted with a direct-response (instantaneous-equilibrium) power-function plasma-to-epithelial-lining-fluid (ELF) link. Plasma disposition is parameterized by CL, V1, V2, V3, Q12 and Q13 with log-normal between-subject variability on every structural parameter. The ELF observable Celf is algebraic, not a compartment: an effect-site equilibration model was fitted and rejected (ELF half-life 8 min, no OFV improvement), so ELF concentration is an instantaneous power function of the total plasma concentration, Celf = relf * Cc^pow_elf, where relf is the ELF:plasma penetration ratio at a plasma concentration of 1 mg/L. Because pow_elf is below 1, penetration falls as plasma concentration rises: 47.2% at Cc = 1 mg/L, 42% at the efficacy-relevant Cc of 2.4 mg/L, and 33.2% at a typical Cmax of 12 mg/L. Avibactam plasma proportional residual noise itself carries between-subject variability (median 0.117, 45% CV); the ELF residual could not be identified from one ELF sample per subject and was FIXED to the plasma median. Penetration ratios are relative to TOTAL plasma and TOTAL ELF concentrations.

The distinguishing feature of this analysis is the link function between plasma and ELF. Section 2.2.1.1 evaluated three forms for each drug:

Link Equation Selected for
Proportional Celf = EPR * Cp neither
Power Celf = EPR(1 mg/L) * Cp^POW avibactam (22 OFV units better than proportional, 6 better than saturable)
Saturable Celf = EMAX * Cp / (KM + Cp) ceftazidime (35 OFV units better than proportional, 8 better than power)

An effect-site equilibration model was fitted for both drugs and rejected: it gave no OFV improvement, and the estimated ELF half-lives were 13 min (ceftazidime, Sect. 3.2) and 8 min (avibactam, Sect. 3.3). Equilibration is therefore treated as instantaneous, and Celf is an algebraic observable rather than an ODE compartment in both model files.

Population

Both models were fitted to the phase I open-label study NCT01395420 (Sect. 2.1). Forty-three healthy male volunteers received either ceftazidime 2000 mg + avibactam 500 mg (cohort A, n = 22) or ceftazidime 3000 mg + avibactam 1000 mg (cohort B, n = 21) as a 2 h intravenous infusion every 8 h for 3 days – nine doses per subject. PK data were available for 42 subjects (22 cohort A, 20 cohort B).

Each subject contributed exactly one bronchoalveolar-lavage ELF sample, which is the reason the ELF residual error is not separately identifiable (Sect. 2.2, Sect. 4).

The paper reports no age, weight, sex-balance or race distribution table and screened no covariates, so neither model file carries a covariateData entry.

str(cef$population)
#> List of 7
#>  $ species       : chr "human"
#>  $ n_subjects    : num 42
#>  $ n_studies     : num 1
#>  $ sex_female_pct: num 0
#>  $ disease_state : chr "healthy volunteers"
#>  $ dose_range    : chr "ceftazidime 2000 mg + avibactam 500 mg (cohort A, n = 22) or ceftazidime 3000 mg + avibactam 1000 mg (cohort B,"| __truncated__
#>  $ notes         : chr "Sect. 2.1: a previously reported phase I open-label study (NCT01395420) enrolled 43 healthy male volunteers; PK"| __truncated__

Source trace

Every ini() entry in both model files carries an in-file comment pointing at its source location. They are collected here for review.

Parameter Ceftazidime Avibactam Source location
lcl (CL, L/h) 6.55 (RSE 2%) 12.5 (RSE 1%) Table 1, CL
lvc (V1, L) 10.32 (RSE 8%) 15.10 (RSE 10%) Table 1, V 1
lvp (V2, L) 5.82 (RSE 9%) 6.52 (RSE 15%) Table 1, V 2
lvp2 (V3, L) 0.64 (RSE 8%) 1.58 (RSE 5%) Table 1, V 3
lq (Q12, L/h) 6.87 (RSE 24%) 5.43 (RSE 42%) Table 1, Q 12
lq2 (Q13, L/h) 0.040 (RSE 33%) 0.14 (RSE 41%) Table 1, Q 13
propSd (RESMp) 0.101 (RSE 9%) 0.117 (RSE 7%) Table 1, RESM p
lemax_elf (EMAX, mg/L) 45.4 (RSE 12%) n/a Table 2, ceftazidime E MAX
lkm_elf (KM, mg/L) 71.7 (RSE 22%) n/a Table 2, ceftazidime K M
lrelf (EPR at 1 mg/L) n/a 0.472 (RSE 11%) Table 2, avibactam EPR (1 mg/l)
lpow_elf (POW) n/a 0.860 (RSE 5%) Table 2, avibactam POW
propSd_Celf (RESM_ELF) 0.101, fixed 0.117, fixed Table 2, RESM ELF; Sect. 2.2 states it was fixed to the plasma median
BSV on CL / V1 / V2 / V3 / Q12 / Q13 CV 10 / 18 / 18 / 25 / 66 / 61% CV 7 / 9 / 25 / 20 / 60 / 75% Table 1, CV % column
BSV on RESMp CV 54% CV 45% Table 1, RESM p row, CV % column
BSV on ELF link parameters EMAX CV 24%, KM CV 97% EPR CV 68%, POW CV 23% Table 2, CV % column
Plasma ODE system (3-compartment, CL/V/Q) n/a n/a Sect. 2.2 (“coded as a series of differential equations”)
Celf saturable link n/a n/a Sect. 2.2.1.1, C ELF (t) = E MAX . C p (t) / (K M + C p (t))
Celf power link n/a n/a Sect. 2.2.1.1, C ELF (t) = EPR(1 mg/L) . [C p (t)]^POW

Between-subject variability is reported in Table 1 / Table 2 as a %CV of a log-normal distribution (Sect. 2.2). The model files therefore carry log(1 + CV^2) as the variance, with the arithmetic left inline so the published CV stays legible. The conversion is confirmed by the tables’ own 5th/95th percentile columns:

# Table 1, ceftazidime V1: median 10.32 L, CV 18%, published 5th-95th 7.73-13.78
om <- sqrt(log(1 + 0.18^2))
round(10.32 * exp(qnorm(c(0.05, 0.95)) * om), 2)
#> [1]  7.69 13.84

# Table 2, ceftazidime KM: median 71.7 mg/L, CV 97%, published 5th-95th 18.8-273.3
om_km <- sqrt(log(1 + 0.97^2))
round(71.7 * exp(qnorm(c(0.05, 0.95)) * om_km), 1)
#> [1]  18.8 273.7

Sects. 3.2 and 3.3 quote six derived quantities. These are deterministic functions of the Table 2 estimates – no simulation, no cohort draw – so they are checked to tight tolerance.

link_cef <- function(cp) 45.4 * cp / (71.7 + cp)          # Sect. 2.2.1.1 saturable
link_avi <- function(cp) 0.472 * cp^0.860                  # Sect. 2.2.1.1 power

link_claims <- tibble::tribble(
  ~Drug,         ~Claim,                                        ~Published, ~Model,
  "Ceftazidime", "ELF:plasma ratio as Cp -> 0 (EMAX/KM, %)",     63.3,  100 * 45.4 / 71.7,
  "Ceftazidime", "Celf at Cp = 70 mg/L (mg/L)",                  22.5,  link_cef(70),
  "Ceftazidime", "Penetration at Cp = 70 mg/L (%)",              32.1,  100 * link_cef(70) / 70,
  "Ceftazidime", "Cp giving Celf = 8 mg/L -> Celf there (mg/L)",  8.0,  link_cef(15.3),
  "Ceftazidime", "Penetration at Cp = 15.3 mg/L (%)",            52.0,  100 * link_cef(15.3) / 15.3,
  "Avibactam",   "Penetration at Cp = 1 mg/L (%)",               47.0,  100 * link_avi(1),
  "Avibactam",   "Celf at Cp = 2.4 mg/L (mg/L)",                  1.0,  link_avi(2.4),
  "Avibactam",   "Penetration at Cp = 2.4 mg/L (%)",             42.0,  100 * link_avi(2.4) / 2.4,
  "Avibactam",   "Celf at Cp = 12 mg/L (mg/L)",                   4.0,  link_avi(12),
  "Avibactam",   "Penetration at Cp = 12 mg/L (%)",              33.2,  100 * link_avi(12) / 12
) |>
  mutate(`Difference (%)` = 100 * (Model - Published) / Published)

link_claims |>
  mutate(across(c(Published, Model, `Difference (%)`), \(x) round(x, 2))) |>
  knitr::kable(caption = "Sects. 3.2 / 3.3 derived quantities vs. the packaged link functions.")
Sects. 3.2 / 3.3 derived quantities vs. the packaged link functions.
Drug Claim Published Model Difference (%)
Ceftazidime ELF:plasma ratio as Cp -> 0 (EMAX/KM, %) 63.3 63.32 0.03
Ceftazidime Celf at Cp = 70 mg/L (mg/L) 22.5 22.43 -0.32
Ceftazidime Penetration at Cp = 70 mg/L (%) 32.1 32.04 -0.19
Ceftazidime Cp giving Celf = 8 mg/L -> Celf there (mg/L) 8.0 7.98 -0.20
Ceftazidime Penetration at Cp = 15.3 mg/L (%) 52.0 52.18 0.35
Avibactam Penetration at Cp = 1 mg/L (%) 47.0 47.20 0.43
Avibactam Celf at Cp = 2.4 mg/L (mg/L) 1.0 1.00 0.21
Avibactam Penetration at Cp = 2.4 mg/L (%) 42.0 41.76 -0.58
Avibactam Celf at Cp = 12 mg/L (mg/L) 4.0 4.00 -0.01
Avibactam Penetration at Cp = 12 mg/L (%) 33.2 33.33 0.40

# These are algebra on the published Table 2 estimates, not a simulated cohort,
# so the tolerance is tight. The residual 0-1% is the paper's own rounding of
# its printed intermediates (e.g. 22.5 for an exact 22.43).
stopifnot(max(abs(link_claims$`Difference (%)`)) < 1.5)

The model files must give the same numbers as the hand-written link functions above. Checked against a typical-value solve so the identity is exact:

solve_typical <- function(mod, amt) {
  ev <- rxode2::et(amt = amt, dur = 2, cmt = "central") |>
    rxode2::et(seq(0, 8, by = 0.02), cmt = "Cc")
  as.data.frame(rxode2::rxSolve(rxode2::zeroRe(mod), ev, useLinCmt = FALSE))
}
tv_cef <- solve_typical(cef, 2000)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalemax_elf', 'etalkm_elf', 'etapropSd'
tv_avi <- solve_typical(avi, 500)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalrelf', 'etalpow_elf', 'etapropSd'

# Celf computed inside the model must equal the published link function applied
# to the model's own Cc, at every time point.
stopifnot(
  max(abs(tv_cef$Celf - link_cef(tv_cef$Cc))) < 1e-8,
  max(abs(tv_avi$Celf - link_avi(tv_avi$Cc))) < 1e-8
)

Virtual cohort

Original observed data are not publicly available. The cohort below mirrors the NCT01395420 design (Sect. 2.1): two dose groups, a 2 h intravenous infusion every 8 h, nine doses. 150 subjects per arm – below the 200/arm cap, and ample for the percentile bands the paper plots.

# rxode2's RNG streams are partitioned per solver thread, so this cohort is
# reproducible on a given machine but NOT across thread counts. Every assertion
# below is written to hold for any cohort these models can produce.
rxode2::rxSetSeed(20180727)

n_per_arm <- 150L
tau <- 8            # dosing interval (h)
n_doses <- 9        # Sect. 2.1: nine doses per subject
t_last <- tau * (n_doses - 1)   # 64 h: start of the final infusion

# Observe densely over the FIRST and the FINAL dosing interval: the first gives
# a single-dose NCA comparison, the final gives the steady-state one the paper
# simulated (Fig. 1 is plotted "versus time from start of infusion of the last
# dose").
obs_grid <- sort(unique(c(seq(0, tau, by = 0.1), seq(t_last, t_last + tau, by = 0.1))))

make_events <- function(amt, n, id_offset, cohort) {
  rxode2::et(amt = amt, dur = 2, ii = tau, addl = n_doses - 1L, cmt = "central") |>
    rxode2::et(obs_grid, cmt = "Cc") |>
    rxode2::et(id = id_offset + seq_len(n)) |>
    as.data.frame() |>
    dplyr::mutate(cohort = cohort)
}

ev_cef <- dplyr::bind_rows(
  make_events(2000, n_per_arm, 0L,           "A: ceftazidime 2000 mg"),
  make_events(3000, n_per_arm, n_per_arm,    "B: ceftazidime 3000 mg")
)
ev_avi <- dplyr::bind_rows(
  make_events(500,  n_per_arm, 0L,           "A: avibactam 500 mg"),
  make_events(1000, n_per_arm, n_per_arm,    "B: avibactam 1000 mg")
)

# Disjoint IDs across arms are mandatory -- duplicated IDs silently merge into
# one subject receiving the summed dose. Both checks are on the RAW frame: a
# `unique()` before `anyDuplicated()` would remove the very duplicates being
# looked for and give a gate that can never go red.
stopifnot(
  anyDuplicated(ev_cef[, c("id", "time", "evid")]) == 0L,
  anyDuplicated(ev_avi[, c("id", "time", "evid")]) == 0L,
  dplyr::n_distinct(ev_cef$id) == 2L * n_per_arm,
  dplyr::n_distinct(ev_avi$id) == 2L * n_per_arm
)

Simulation

Two rxode2 details are load-bearing for these models.

useLinCmt = FALSE is required. The default ODE-to-linCmt() conversion silently drops the second peripheral compartment of a k12/k21/k13/k31 parameterisation and mis-maps the two endpoints onto compartment slots.

Observation rows use cmt = "Cc", not cmt = "central". Both models declare two endpoints (Cc ~ prop(...) and Celf ~ prop(...)), so rxode2 builds a dvid-to-cmt map onto the endpoint slots that follow the three ODE states; an observation row pointing at central is then rejected outright. Because Cc and Celf are already declared endpoints, naming one on an observation row appends nothing and renumbers nothing – the usual hazard of writing an algebraic observable into cmt does not apply here. The assertion below pins that: the ODE states must still be the three named compartments in order, and peripheral2 must actually carry drug.

sim_cef <- rxode2::rxSolve(cef, events = ev_cef, keep = "cohort", useLinCmt = FALSE) |>
  as.data.frame()
sim_avi <- rxode2::rxSolve(avi, events = ev_avi, keep = "cohort", useLinCmt = FALSE) |>
  as.data.frame()

stopifnot(
  # No compartment renumbering: the states, in order, are exactly the three
  # declared by d/dt(), and the deep compartment is populated (it would be
  # absent or identically zero if slot 3 had been displaced).
  identical(cef$state, c("central", "peripheral1", "peripheral2")),
  identical(avi$state, c("central", "peripheral1", "peripheral2")),
  max(sim_cef$peripheral2, na.rm = TRUE) > 0,
  max(sim_avi$peripheral2, na.rm = TRUE) > 0,
  all(c("Cc", "Celf") %in% names(sim_cef)),
  dplyr::n_distinct(sim_cef$id) == 2L * n_per_arm,
  all(sim_cef$Cc >= 0), all(sim_avi$Cc >= 0)
)

Replicate published figures

Figure 2 – simulated ELF versus plasma concentration

Figure 2 plots simulated ELF against plasma concentration for 1000 subjects receiving ceftazidime/avibactam 2000-500 mg q8h, with the population median and the 5th / 95th percentiles. The ceftazidime panel is visibly saturable above about 100 mg/L; the avibactam panel departs from linearity only mildly, above about 20 mg/L (Sect. 3.5).

# Replicates Figure 2 of Dimelow 2018 (cohort A dose only, as in the paper).
fig2 <- dplyr::bind_rows(
  sim_cef |> filter(cohort == "A: ceftazidime 2000 mg") |> mutate(drug = "a  Ceftazidime"),
  sim_avi |> filter(cohort == "A: avibactam 500 mg")    |> mutate(drug = "b  Avibactam")
) |>
  filter(!is.na(Cc), Cc > 0) |>
  group_by(drug) |>
  mutate(cp_bin = cut(Cc, breaks = 10^seq(-2, 2.6, by = 0.1), labels = FALSE)) |>
  group_by(drug, cp_bin) |>
  summarise(
    Cp = median(Cc),
    Q05 = quantile(Celf, 0.05), Q50 = quantile(Celf, 0.50), Q95 = quantile(Celf, 0.95),
    .groups = "drop"
  )

ggplot(fig2, aes(Cp, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~drug, scales = "free") +
  scale_x_log10() + scale_y_log10() +
  labs(
    x = "Total plasma concentration (mg/L)", y = "Total ELF concentration (mg/L)",
    title = "Figure 2 -- simulated ELF vs. plasma concentration",
    caption = "Replicates Figure 2 of Dimelow 2018 (2000-500 mg q8h). Line = median, band = 5th-95th percentile."
  )

Figure 3 – simulated ELF concentration-time profiles vs. the PK/PD targets

Figure 3 overlays the simulated ELF profiles on the plasma-derived PK/PD targets (50% fT > 8 mg/L for ceftazidime, 50% fT > 1 mg/L for avibactam) with a vertical line at the 4 h midpoint of the 8 h dosing interval.

# Replicates Figure 3 of Dimelow 2018: final dosing interval, cohort A dose.
targets <- c("a  Ceftazidime" = 8, "b  Avibactam" = 1)

fig3 <- dplyr::bind_rows(
  sim_cef |> filter(cohort == "A: ceftazidime 2000 mg") |> mutate(drug = "a  Ceftazidime"),
  sim_avi |> filter(cohort == "A: avibactam 500 mg")    |> mutate(drug = "b  Avibactam")
) |>
  filter(time >= t_last, !is.na(Celf)) |>
  mutate(tad = time - t_last) |>
  group_by(drug, tad) |>
  summarise(
    Q05 = quantile(Celf, 0.05), Q50 = quantile(Celf, 0.50), Q95 = quantile(Celf, 0.95),
    .groups = "drop"
  ) |>
  mutate(target = targets[drug])

ggplot(fig3, aes(tad, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2) +
  geom_line(linewidth = 0.8) +
  geom_hline(aes(yintercept = target), linetype = "dashed", colour = "firebrick") +
  geom_vline(xintercept = 4, linetype = "dotted") +
  facet_wrap(~drug, scales = "free_y") +
  scale_y_log10() +
  labs(
    x = "Time from start of the final infusion (h)", y = "Total ELF concentration (mg/L)",
    title = "Figure 3 -- simulated ELF profiles vs. plasma PK/PD targets",
    caption = paste(
      "Replicates Figure 3 of Dimelow 2018 (2000-500 mg q8h). Dashed line = PK/PD target,",
      "dotted line = 4 h dosing-interval midpoint."
    )
  )

Section 3.5 – ELF concentration at the dosing-interval midpoint

Section 3.5 reports the one fully quantitative simulation result in the paper: at 4 h post dose, median (5th-95th percentile) simulated ELF concentrations were 14.7 (7.1-23.4) mg/L for ceftazidime and 1.4 (0.75-2.85) mg/L for avibactam, for 1000 subjects on 2000-500 mg q8h.

midpoint <- function(sim, arm) {
  x <- sim |> filter(cohort == arm, abs(time - (t_last + 4)) < 1e-6)
  c(median = median(x$Celf), p05 = unname(quantile(x$Celf, 0.05)),
    p95 = unname(quantile(x$Celf, 0.95)), n = nrow(x))
}
mid_cef <- midpoint(sim_cef, "A: ceftazidime 2000 mg")
mid_avi <- midpoint(sim_avi, "A: avibactam 500 mg")
stopifnot(mid_cef[["n"]] == n_per_arm, mid_avi[["n"]] == n_per_arm)  # the gate had rows to test

tibble::tibble(
  Drug = c("Ceftazidime", "Avibactam"),
  `Published median` = c(14.7, 1.4),
  `Model median` = c(mid_cef[["median"]], mid_avi[["median"]]),
  `Published 5th-95th` = c("7.1-23.4", "0.75-2.85"),
  `Model 5th-95th` = sprintf(
    "%.2f-%.2f",
    c(mid_cef[["p05"]], mid_avi[["p05"]]), c(mid_cef[["p95"]], mid_avi[["p95"]])
  )
) |>
  mutate(`Difference in median (%)` = round(100 * (`Model median` - `Published median`) / `Published median`, 1),
         `Model median` = round(`Model median`, 2)) |>
  knitr::kable(caption = "Sect. 3.5, ELF concentration 4 h after the dose at steady state.")
Sect. 3.5, ELF concentration 4 h after the dose at steady state.
Drug Published median Model median Published 5th-95th Model 5th-95th Difference in median (%)
Ceftazidime 14.7 14.22 7.1-23.4 4.60-34.32 -3.2
Avibactam 1.4 1.53 0.75-2.85 0.54-4.48 9.6

# The MEDIAN is the structural quantity: a mis-transcribed clearance, volume,
# dose or link parameter moves it by tens of percent. It is checked; the
# 5th-95th spread deliberately is NOT (see Assumptions -- the published
# plasma OMEGA is a full covariance matrix whose off-diagonals were not
# printed, so the packaged diagonal block over-disperses the tails).
#
# Bounds are set OUTSIDE the range seen across two independent cohort draws,
# not tightened onto one run: ceftazidime realised -0.2% and +5.9%, avibactam
# +6.6% and -3.6%. Do not tighten these back -- rxode2 draws a different
# cohort on a machine with a different solver-thread count.
stopifnot(
  abs(mid_cef[["median"]] - 14.7) / 14.7 < 0.20,
  abs(mid_avi[["median"]] -  1.4) /  1.4 < 0.25
)

Section 3.5 also claims that “most subjects still achieved ceftazidime and avibactam ELF exposures exceeding their respective plasma PK/PD targets beyond the midpoint of the dosing interval.”

pct_over <- function(sim, arm, target) {
  x <- sim |> filter(cohort == arm, abs(time - (t_last + 4)) < 1e-6)
  stopifnot(nrow(x) == n_per_arm)
  100 * mean(x$Celf > target)
}
attain <- c(
  Ceftazidime = pct_over(sim_cef, "A: ceftazidime 2000 mg", 8),
  Avibactam   = pct_over(sim_avi, "A: avibactam 500 mg",    1)
)
round(attain, 1)
#> Ceftazidime   Avibactam 
#>        76.0        66.7

# "Most subjects" = a majority, which is the paper's own claim and therefore
# the bound. Realised across two draws: ceftazidime 80-84%, avibactam 66-71%.
# With n = 150 the binomial SE is about 4 points, so the 50% floor sits several
# SE below either arm, and the gate still goes red if a link or disposition
# parameter is wrong (halving the penetration ratio would drop ceftazidime
# below 50%).
stopifnot(all(attain > 50))

PKNCA validation

NCA is computed with PKNCA over two intervals per subject: the first dosing interval (0-8 h) and the final, steady-state interval (64-72 h). Both outputs are run separately – plasma Cc and ELF Celf.

tau_intervals <- data.frame(
  start   = c(0, t_last),
  end     = c(tau, t_last + tau),
  cmax    = TRUE,
  tmax    = TRUE,
  cmin    = TRUE,
  auclast = TRUE,
  cav     = TRUE
)

dose_frame <- function(ev) {
  ev |>
    filter(evid != 0) |>
    distinct(id, cohort) |>
    tidyr::crossing(time = tau * seq(0, n_doses - 1L)) |>
    mutate(amt = 1) |>          # amount is irrelevant to the parameters requested
    arrange(id, time)
}

run_nca <- function(sim, ev, col) {
  conc <- sim |>
    dplyr::filter(!is.na(.data[[col]])) |>       # ONLY !is.na -- never time > 0 / conc > 0
    dplyr::select(id, cohort, time, conc = dplyr::all_of(col))
  stopifnot(nrow(conc) > 0, any(conc$time == 0), any(conc$time == t_last))
  co <- PKNCA::PKNCAconc(conc, conc ~ time | cohort + id)
  do <- PKNCA::PKNCAdose(as.data.frame(dose_frame(ev)), amt ~ time | cohort + id)
  suppressWarnings(PKNCA::pk.nca(PKNCA::PKNCAdata(co, do, intervals = tau_intervals)))
}

nca_cef_plasma <- run_nca(sim_cef, ev_cef, "Cc")
nca_cef_elf    <- run_nca(sim_cef, ev_cef, "Celf")
nca_avi_plasma <- run_nca(sim_avi, ev_avi, "Cc")
nca_avi_elf    <- run_nca(sim_avi, ev_avi, "Celf")
tidy_nca <- function(res, drug, matrix) {
  as.data.frame(res$result) |>
    mutate(drug = drug, matrix = matrix,
           interval = ifelse(start == 0, "First dose (0-8 h)", "Steady state (64-72 h)"))
}
nca_all <- dplyr::bind_rows(
  tidy_nca(nca_cef_plasma, "Ceftazidime", "Plasma"),
  tidy_nca(nca_cef_elf,    "Ceftazidime", "ELF"),
  tidy_nca(nca_avi_plasma, "Avibactam",   "Plasma"),
  tidy_nca(nca_avi_elf,    "Avibactam",   "ELF")
)

nca_all |>
  filter(PPTESTCD %in% c("cmax", "tmax", "cmin", "auclast", "cav")) |>
  group_by(drug, matrix, cohort, interval, PPTESTCD) |>
  summarise(value = median(PPORRES, na.rm = TRUE), .groups = "drop") |>
  mutate(PPTESTCD = nlmixr2lib::ncaParamLabel(PPTESTCD)) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = value) |>
  mutate(across(where(is.numeric), \(x) signif(x, 3))) |>
  dplyr::rename("Drug" = drug, "Matrix" = matrix, "Arm" = cohort, "Interval" = interval) |>
  knitr::kable(caption = "Median per-subject NCA of the simulated cohort.")
Median per-subject NCA of the simulated cohort.
Drug Matrix Arm Interval AUClast Cavg Cmax Cmin Tmax
Avibactam ELF A: avibactam 500 mg First dose (0-8 h) 15.1 1.89 5.07 0.000 2
Avibactam ELF A: avibactam 500 mg Steady state (64-72 h) 15.9 1.98 5.21 0.300 2
Avibactam ELF B: avibactam 1000 mg First dose (0-8 h) 28.1 3.52 8.94 0.000 2
Avibactam ELF B: avibactam 1000 mg Steady state (64-72 h) 29.1 3.64 9.05 0.610 2
Avibactam Plasma A: avibactam 500 mg First dose (0-8 h) 38.4 4.80 14.00 0.000 2
Avibactam Plasma A: avibactam 500 mg Steady state (64-72 h) 39.9 4.98 14.30 0.597 2
Avibactam Plasma B: avibactam 1000 mg First dose (0-8 h) 77.0 9.62 28.40 0.000 2
Avibactam Plasma B: avibactam 1000 mg Steady state (64-72 h) 80.6 10.10 29.10 1.230 2
Ceftazidime ELF A: ceftazidime 2000 mg First dose (0-8 h) 106.0 13.30 24.30 0.000 2
Ceftazidime ELF A: ceftazidime 2000 mg Steady state (64-72 h) 111.0 13.80 24.90 4.730 2
Ceftazidime ELF B: ceftazidime 3000 mg First dose (0-8 h) 143.0 17.90 29.70 0.000 2
Ceftazidime ELF B: ceftazidime 3000 mg Steady state (64-72 h) 150.0 18.70 30.10 7.540 2
Ceftazidime Plasma A: ceftazidime 2000 mg First dose (0-8 h) 282.0 35.30 90.60 0.000 2
Ceftazidime Plasma A: ceftazidime 2000 mg Steady state (64-72 h) 308.0 38.50 94.20 8.300 2
Ceftazidime Plasma B: ceftazidime 3000 mg First dose (0-8 h) 422.0 52.80 135.00 0.000 2
Ceftazidime Plasma B: ceftazidime 3000 mg Steady state (64-72 h) 458.0 57.30 142.00 12.800 2

ELF:plasma AUC ratio vs. the published non-compartmental result

The phase I study’s own non-compartmental analysis – the comparator this paper was written to improve on – found the ELF AUC to be 31-35% of the plasma AUC for each drug (Key Points; Sect. 4 quotes about 32% for ceftazidime and about 35% for avibactam). That ratio is the single best end-to-end check on the packaged models: it exercises the plasma disposition and the link function together, over a whole dosing interval.

auc_ratio <- nca_all |>
  filter(PPTESTCD == "auclast", interval == "Steady state (64-72 h)") |>
  select(drug, matrix, cohort, id, PPORRES) |>
  tidyr::pivot_wider(names_from = matrix, values_from = PPORRES) |>
  mutate(ratio = 100 * ELF / Plasma)

auc_summary <- auc_ratio |>
  group_by(drug, cohort) |>
  summarise(`Median ELF:plasma AUC ratio (%)` = round(median(ratio), 1), .groups = "drop") |>
  mutate(`Published NCA ratio (%)` = "31-35")

auc_summary |>
  dplyr::rename("Drug" = drug, "Arm" = cohort) |>
  knitr::kable(caption = "Steady-state AUC0-tau ratio, ELF vs. plasma, against the phase I NCA result.")
Steady-state AUC0-tau ratio, ELF vs. plasma, against the phase I NCA result.
Drug Arm Median ELF:plasma AUC ratio (%) Published NCA ratio (%)
Avibactam A: avibactam 500 mg 40.0 31-35
Avibactam B: avibactam 1000 mg 34.1 31-35
Ceftazidime A: ceftazidime 2000 mg 37.9 31-35
Ceftazidime B: ceftazidime 3000 mg 31.6 31-35

# The published band is 31-35%; the gate allows 25-45% so that it tolerates the
# cohort draw while still going red on a mis-transcribed EMAX, KM, EPR, POW, CL
# or V1. Realised across two independent draws of the four arms: 31.3-38.2%.
# Only the MEDIAN is gated -- the per-subject spread is inflated by the dropped
# OMEGA off-diagonals, so its tails are not a reproducible quantity.
stopifnot(all(auc_summary$`Median ELF:plasma AUC ratio (%)` > 25),
          all(auc_summary$`Median ELF:plasma AUC ratio (%)` < 45))

Both drugs sit inside the published band, and the ratio falls from cohort A to cohort B for each drug – the paper’s central finding that penetration is concentration-dependent and therefore higher at the lower, efficacy-relevant plasma concentrations.

Comparison against published NCA

The paper reports no NCA table of its own. It does quote two literature plasma Cmax values while interpreting the link functions (Sects. 3.2 and 3.3, reference [26]): about 70 mg/L for ceftazidime 2000 mg and about 12 mg/L for avibactam 500 mg. Those are compared here against the first-dose simulated Cmax of the packaged models.

sim_long <- nca_all |>
  filter(matrix == "Plasma", interval == "First dose (0-8 h)",
         cohort %in% c("A: ceftazidime 2000 mg", "A: avibactam 500 mg")) |>
  mutate(drug = ifelse(drug == "Ceftazidime", "Ceftazidime 2000 mg", "Avibactam 500 mg")) |>
  select(drug, PPTESTCD, PPORRES)

published <- tibble::tribble(
  ~drug,                  ~cmax,
  "Ceftazidime 2000 mg",  70,
  "Avibactam 500 mg",     12
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = sim_long,
  reference     = published,
  by            = "drug",
  params        = "cmax",
  units         = c(cmax = "mg/L"),
  tolerance_pct = 20
)
cmp |>
  dplyr::rename("Drug" = drug) |>
  knitr::kable(caption = "Simulated first-dose plasma Cmax vs. the literature values quoted by Dimelow 2018. * differs by >20%.")
Simulated first-dose plasma Cmax vs. the literature values quoted by Dimelow 2018. * differs by >20%.
NCA parameter Drug Reference Simulated % diff
Cmax (mg/L) Ceftazidime 2000 mg 70 90.6 +29.4%*
Cmax (mg/L) Avibactam 500 mg 12 14 +16.8%

The ceftazidime row is starred and is not tuned away. The discrepancy is internal to the paper: Table 1 gives V1 = 10.32 L and CL = 6.55 L/h, which for a 2000 mg 2 h infusion imply an end-of-infusion concentration near 90 mg/L, while Sect. 3.2 quotes about 70 mg/L from an external reference to argue that KM (71.7 mg/L) is “comparable to the typical Cmax”. Both statements are in the source; only the Table 1 parameters are the model, so the packaged model follows Table 1. The avibactam row agrees within tolerance.

Assumptions and deviations

  • Plasma OMEGA is carried as diagonal. Sect. 3.1 states that between-subject variability “was adequately characterized using a full covariance matrix”, but no off-diagonal element is printed anywhere in the paper, and inventing one is not an option. Each parameter’s published %CV is therefore carried on the diagonal and all correlations are set to zero. The consequence is visible and one-directional: medians are unaffected (they are reproduced to within 6% in the Sect. 3.5 check above), but the simulated 5th-95th percentile band is wider than the published one – 5.2-33.1 mg/L against the paper’s 7.1-23.4 for ceftazidime. Every assertion in this vignette is therefore written on medians, never on the tails.
  • Between-subject variability on the residual error is encoded, and is unusual. Table 1 reports RESMp as a population median (0.101 ceftazidime, 0.117 avibactam) with its own %CV (54% and 45%), and its 5th/95th percentile columns confirm a log-normal spread. The models therefore scale the proportional residual SD per subject (propSdi <- propSd * exp(etapropSd)) rather than treating it as a single population constant.
  • ELF residual error is fixed, not estimated. Sect. 2.2: with one ELF observation per individual the ELF residual could not be identified and was fixed to the plasma median. Table 2 confirms RESM ELF carries CV 0%. Sect. 4 notes the corollary the authors themselves flag: the ELF between-subject variance estimates are conditional on that assumed residual, and would shrink if a larger residual had been assumed.
  • Penetration ratios are relative to TOTAL concentrations in both matrices. Sect. 4 is explicit about this and about why it is conservative: plasma free fraction is about 85% (ceftazidime) and 92% (avibactam), ELF protein content is much lower than plasma, so the free fraction in ELF is likely higher than in plasma. Neither model file applies a protein-binding correction, and Celf should not be read as a free concentration without one. The PK/PD targets plotted in Figure 3 are plasma-derived (mouse data, references [10-12]) and the paper states they likely overestimate what is required in ELF.
  • Celf is an algebraic observable, not a compartment. This follows the paper’s model selection (Sect. 3.2, Sect. 3.3), not a simplification: the effect-site equilibration model was fitted, gave no OFV improvement, and implied ELF half-lives of 13 min (ceftazidime) and 8 min (avibactam). The KELF estimates behind those half-lives are quoted in prose only and are not tabulated, so the rejected effect-site model is not packaged.
  • No covariates. The paper screened none and reports no demographic table beyond “43 healthy male volunteers”, so neither model file declares covariateData and the virtual cohort carries no covariate columns.
  • Dosing history is reconstructed, not published per subject. Sect. 2.1 gives nine 2 h infusions every 8 h; the vignette’s event table reproduces exactly that. The paper’s own simulations (Sect. 3.5, Figs. 2-3) used 1000 subjects on the 2000-500 mg regimen; this vignette uses 150 per arm, which is sufficient for a median comparison and stays under the 200/arm cap.
  • Literature Cmax comparison is external to this paper. The 70 mg/L and 12 mg/L values in the NCA comparison table come from reference [26] as quoted by Dimelow 2018, not from an NCA the authors performed. See the discussion of the starred row above.
  • useLinCmt = FALSE is mandatory for every rxSolve() call against these models. The default ODE-to-linCmt() conversion drops peripheral2 from a k13/k31 parameterisation and mis-maps the two endpoints onto compartment slots.