Skip to contents

Cefotiam (Shah 2019)

Source: vignettes/articles/Shah_2019_cefotiam.Rmd
Shah_2019_cefotiam.Rmd

Model and source

  • Citation: Shah NR, Bulitta JB, Kinzig M, Landersdorfer CB, Jiao Y, Sutaria DS, Tao X, Hohl R, Holzgrabe U, Kees F, Stephan U, Sorgel F. Novel Population Pharmacokinetic Approach to Explain the Differences between Cystic Fibrosis Patients and Healthy Volunteers via Protein Binding. Pharmaceutics. 2019 Jun 18;11(6):286. doi:10.3390/pharmaceutics11060286
  • Description: Three-compartment population PK model for IV cefotiam with a simultaneous fit of total plasma concentrations and the fraction of dose excreted unchanged in urine. Built from 14 Caucasian adults (8 cystic fibrosis [CF] patients, 6 healthy volunteers [HVs]) each given a single 1027.5 mg cefotiam dose as a 3 min IV infusion. The novel feature of this model is that every disposition parameter is referenced to the UNBOUND cefotiam concentration, and the observed TOTAL plasma concentration is reconstructed algebraically as Cc = Cunbound / fu. Body size and body composition are captured by allometric scaling on lean body mass (LBM) with fixed exponents 0.75 on all clearance terms and 1.0 on all volumes (reference LBM = 53 kg, equivalent to a standard 70 kg total body weight). Unbound total clearance is split into a renal arm (CL_R,u, whose mass is tracked in the canonical urine compartment) and a non-renal arm (CL_NR,u). Because the unbound parameters are shared across all four cohorts, the entire remaining CF-vs-HV and female-vs-male difference is carried by the unbound fraction fu, which is estimated separately for each combination of the DIS_CF and SEXF indicators: fu = 0.500 (fixed, HV males), 0.545 (HV females), 0.563 (CF males) and 0.744 (CF females). This replaces the disease-specific scale factors (FCYF) used by the same group in Bulitta_2011_cefpirome.R and Bulitta_2007_piperacillin.R; the authors removed all FCYF terms when the unbound fractions were estimated.
  • Article: https://doi.org/10.3390/pharmaceutics11060286

Shah 2019 is the protein-binding counterpart to the disease-scale-factor models this library already carries from the same group (Bulitta_2007_piperacillin, Bulitta_2011_cefpirome). The question in all three papers is the same: why do patients with cystic fibrosis (CF) appear to clear beta-lactams faster than healthy volunteers (HVs)? The earlier papers answered it with explicit disease-specific scale factors (FCYF) multiplying clearance and volume. This paper answers it a different way. Every disposition parameter is referenced to the unbound cefotiam concentration and is shared across all four subject groups; the entire residual CF-vs-HV and female-vs-male difference is absorbed into the plasma unbound fraction fu, which is estimated per group. Paper Section 2.6.5 is explicit that “all disease specific scale factors FCYF were removed from the model” once the unbound fractions were estimated.

The practical consequence, and the thing this vignette sets out to demonstrate, is that two subjects of identical lean body mass but different fu have identical unbound pharmacokinetics and different total pharmacokinetics.

Population

Paper Table 1 describes 14 Caucasian adults, all of whom received a single 1027.5 mg cefotiam dose as a 3 min IV infusion:

  • 8 patients with cystic fibrosis (4 female, 4 male): total body weight 33.0-59.0 kg (median 45.5), lean body mass (LBM) 28.8-46.2 kg (median 40.3), age 17-24 years (median 19), BMI 13.4-19.9 kg/m^2 (median 17.0).
  • 6 healthy volunteers (3 female, 3 male): total body weight 58.0-80.0 kg (median 68.5), LBM 44.6-65.4 kg (median 50.6), age 21-26 years (median 23.5), BMI 20.3-27.9 kg/m^2 (median 22.5).

LBM was computed with the formula of Cheymol and James (Table 1 footnote a) and is the size descriptor the model scales on. The per-sex LBM medians in Table 1 matter for this model because they separate the two mechanisms: female and male CF patients have essentially the same LBM (38.8 vs 40.4 kg) yet female CF patients had consistently larger clearances in the noncompartmental analysis (Table 2). Size cannot explain that; the unbound fraction can.

Plasma was sampled pre-dose, at the end of the 3 min infusion, and at 5, 10, 15, 20, 30, 45, 60 and 90 min plus 2, 3, 4, 5, 6, 8, 12 and 24 h after the end of infusion. Urine was collected over 0-1, 1-2, 2-3, 3-4, 4-5, 5-6, 6-8, 8-12 and 12-24 h from the start of infusion. Plasma concentrations and urinary fractions were fitted simultaneously, which is what identifies the renal and non-renal clearance arms separately.

Source trace

The final model is paper Table 4. The authors also deposited the Berkeley Madonna Monte Carlo source listing for the final model in the Supplementary Materials; that listing prints the estimates to six significant figures and states the ODE system verbatim, so it is used here as the authoritative implementation with the Table 4 rounding quoted alongside.

Model element Source location Value
Three-compartment disposition Figure 1; Results, Structural Model chosen over 2-compartment on likelihood ratio (p < 0.001)
ODE system, urine state Supplement Madonna listing, d/dt(Cent), d/dt(Shal), d/dt(Deep), d/dt(Urin) transcribed verbatim
FSize,V = LBM / 53 Equation 1; Section 2.6.3 exponent fixed at 1.0
FSize,CL = (LBM / 53)^0.75 Equation 2; Section 2.6.3 exponent fixed at 0.75
Individual parameter form Equation 4 (CLru,i = CLru * FSize,CL,i * exp(eta)) no FCYF term
Cc = Cunbound / fu Section 2.6.5; Madonna C1 = Cent/V1cov/FU total from unbound
CLr,u = 23.8 L/h Table 4; Madonna Mean_CLR = 23.7911 SE 6.9%
CLnr,u = 11.0 L/h Table 4; Madonna Mean_CLNR = 10.9943 SE 7.0%
V1u = 15.6 L Table 4; Madonna Mean_V1 = 15.5927 SE 6.5%
V2u = 6.91 L Table 4; Madonna Mean_V2 = 6.90932 SE 14.1%
V3u = 4.56 L Table 4; Madonna Mean_V3 = 4.5646 SE 16.4%
CLd_shallow,u = 13.8 L/h Table 4; Madonna Mean_CLD = 13.8043 SE 15.0%
CLd_deep,u = 1.84 L/h Table 4; Madonna Mean_CLD3 = 1.83622 SE 26.1%
fu HV male = 0.50 Table 4 (fixed); Madonna FU_HVM fixed from literature refs 15-18
fu HV female = 0.545 Table 4; Madonna Mean_FU_HVF = 0.544842 SE 13.6%
fu CF male = 0.563 Table 4; Madonna Mean_FU_CFM = 0.562527 SE 13.5%
fu CF female = 0.744 Table 4; Madonna Mean_FU_CFF = 0.74359 SE 4.5%
BSV on CL/V terms Table 4 BSV column; Madonna CV_* log-scale SD, so omega^2 = BSV^2
BSV on the three estimated fu Table 4 footnote d fixed at 5% CV
SDin = 0.0186 mg/L, SDsl = 0.166 Table 4 additive + proportional on total plasma
UDin = 0.384% Table 4 additive on percent of dose in urine

Two alternative parameterisations are reported in the Supplementary Materials and are deliberately not encoded: Table S1 (FCYF scale factors with fu fixed to 0.5 everywhere) and Table S2 (renal clearance split into a fixed 7.2 L/h glomerular filtration arm plus a binding-independent tubular secretion arm). The paper presents neither as final; Table S2 had a -2x log-likelihood worse by 7.1.

Virtual cohort

Four arms, one per combination of DIS_CF and SEXF, 50 subjects each. LBM is drawn from a log-normal with the 15% coefficient of variation the paper used for its own Monte Carlo covariate model (Section 2.6.7), centred on the per-sex, per-cohort median LBM from Table 1 and truncated to the range Table 1 observed.

rxode2::rxSetSeed(20190618)
set.seed(20190618)

n_per_arm <- 50L

# Table 1: per-cohort, per-sex lean body mass (median [range], kg).
arm_spec <- tibble::tribble(
  ~cohort,          ~DIS_CF, ~SEXF, ~lbm_med, ~lbm_lo, ~lbm_hi,
  "CF female",      1L,      1L,    38.8,     28.8,    45.7,
  "CF male",        1L,      0L,    40.4,     39.6,    46.2,
  "HV female",      0L,      1L,    44.6,     44.6,    45.4,
  "HV male",        0L,      0L,    62.8,     55.8,    65.4
)

draw_lbm <- function(n, med, lo, hi, cv = 0.15) {
  # Log-normal around the Table 1 median, truncated to the Table 1 range.
  out <- pmin(pmax(med * exp(stats::rnorm(n, 0, cv)), lo), hi)
  out
}

cohorts <- arm_spec |>
  tidyr::uncount(n_per_arm) |>
  dplyr::group_by(cohort) |>
  dplyr::mutate(LBM = draw_lbm(dplyr::n(), lbm_med[1], lbm_lo[1], lbm_hi[1])) |>
  dplyr::ungroup() |>
  dplyr::mutate(
    id = dplyr::row_number(),
    DOSE_CEFOTIAM_MG = 1027.5
  ) |>
  dplyr::select(id, cohort, DIS_CF, SEXF, LBM, DOSE_CEFOTIAM_MG)

cohorts |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(
    n = dplyr::n(),
    `LBM median` = round(stats::median(LBM), 1),
    `LBM min` = round(min(LBM), 1),
    `LBM max` = round(max(LBM), 1),
    .groups = "drop"
  ) |>
  knitr::kable(caption = "Simulated cohort lean body mass by arm, against Table 1.")
Simulated cohort lean body mass by arm, against Table 1.
cohort n LBM median LBM min LBM max
CF female 50 38.8 28.8 45.7
CF male 50 41.4 39.6 46.2
HV female 50 44.6 44.6 45.4
HV male 50 61.9 55.8 65.4

Events

Each subject receives the study dose: 1027.5 mg as a 3 min IV infusion. The observation grid reproduces the paper’s sampling schedule exactly, with plasma times expressed relative to the end of infusion as the paper states them. Plasma observations are placed on cmt = "Cc" and urine observations on cmt = "urinePct", because this model declares two residual endpoints.

dose_amt_mg <- 1027.5
infusion_dur_h <- 3 / 60

# Section 2.3: times are stated post END of infusion.
post_eoi_h <- c(c(0, 5, 10, 15, 20, 30, 45, 60, 90) / 60, 2, 3, 4, 5, 6, 8, 12, 24)
plasma_times <- sort(unique(c(0, infusion_dur_h + post_eoi_h)))

# Section 2.4: urine collection interval ends, from the START of infusion.
urine_times <- c(1, 2, 3, 4, 5, 6, 8, 12, 24)

dose_rows <- cohorts |>
  dplyr::mutate(
    time = 0,
    evid = 1L,
    amt = dose_amt_mg,
    cmt = "central",
    rate = dose_amt_mg / infusion_dur_h
  )

obs_rows <- cohorts |>
  tidyr::crossing(time = plasma_times) |>
  dplyr::mutate(evid = 0L, amt = NA_real_, cmt = "Cc", rate = 0)

urine_rows <- cohorts |>
  tidyr::crossing(time = urine_times) |>
  dplyr::mutate(evid = 0L, amt = NA_real_, cmt = "urinePct", rate = 0)

events <- dplyr::bind_rows(dose_rows, obs_rows, urine_rows) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

stopifnot(!anyDuplicated(events[, c("id", "time", "evid", "cmt")]))

Simulation

useLinCmt = FALSE is required: a model with more than one declared endpoint has its dvid table corrupted by the linear-compartment solver, which makes every cmt choice fail.

mod <- readModelDb("Shah_2019_cefotiam")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep = c("cohort", "LBM", "DIS_CF", "SEXF"),
  useLinCmt = FALSE
) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'

sim_unique <- sim |>
  dplyr::distinct(id, time, .keep_all = TRUE)

Exact structural checks

These three checks compare a deterministic solve against a closed form derived from the published parameters. Both sides use the same numbers, so the only difference is integration error and a tight bound is the correct assertion. Random effects are suppressed with omega = NA, sigma = NA rather than zeroRe(), which segfaults on multi-endpoint models.

typ_cov <- arm_spec |>
  dplyr::mutate(
    id = dplyr::row_number(),
    LBM = lbm_med,
    DOSE_CEFOTIAM_MG = dose_amt_mg
  ) |>
  dplyr::select(id, cohort, DIS_CF, SEXF, LBM, DOSE_CEFOTIAM_MG)

typ_times <- sort(unique(c(
  seq(0, 0.2, by = 0.005),
  seq(0.25, 6, by = 0.05),
  seq(6.25, 48, by = 0.25)
)))

typ_events <- dplyr::bind_rows(
  typ_cov |>
    dplyr::mutate(
      time = 0, evid = 1L, amt = dose_amt_mg,
      cmt = "central", rate = dose_amt_mg / infusion_dur_h
    ),
  typ_cov |>
    tidyr::crossing(time = typ_times) |>
    dplyr::mutate(evid = 0L, amt = NA_real_, cmt = "Cc", rate = 0),
  typ_cov |>
    tidyr::crossing(time = typ_times) |>
    dplyr::mutate(evid = 0L, amt = NA_real_, cmt = "urinePct", rate = 0)
) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

typ <- rxode2::rxSolve(
  mod, typ_events,
  keep = c("cohort", "LBM", "DIS_CF", "SEXF"),
  omega = NA, sigma = NA,
  useLinCmt = FALSE
) |>
  as.data.frame() |>
  dplyr::distinct(id, time, .keep_all = TRUE)

Check 1: total-scale clearance equals fu * CL_unbound * (LBM/53)^0.75

Because the volumes and clearances are both unbound-referenced, the classical total-concentration clearance recovered by NCA must be fu times the unbound clearance. This is the identity the whole paper rests on.

cl_ru <- 23.7911
cl_nru <- 10.9943
v1u <- 15.5927
v2u <- 6.90932
v3u <- 4.5646

fu_lookup <- c(
  "CF female" = 0.74359, "CF male" = 0.562527,
  "HV female" = 0.544842, "HV male" = 0.5
)

closed_form <- typ_cov |>
  dplyr::mutate(
    fu = unname(fu_lookup[cohort]),
    cl_total_pred = fu * (cl_ru + cl_nru) * (LBM / 53)^0.75,
    vss_total_pred = fu * (v1u + v2u + v3u) * (LBM / 53)
  )

# Dose / AUCinf is exactly CL for any linear system, whatever the input
# function, so PKNCA's cl.obs on the deterministic profile is a direct read of
# the model's total-concentration clearance.
typ_nca_in <- typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort) |>
  dplyr::arrange(id, time)

typ_conc <- PKNCA::PKNCAconc(
  typ_nca_in, Cc ~ time | cohort + id,
  concu = "mg/L", timeu = "h"
)
typ_dose <- PKNCA::PKNCAdose(
  typ_cov |> dplyr::mutate(time = 0, amt = dose_amt_mg) |>
    dplyr::select(id, cohort, time, amt),
  amt ~ time | cohort + id,
  doseu = "mg"
)

cl_obs_typ <- suppressWarnings(
  PKNCA::pk.nca(PKNCA::PKNCAdata(
    typ_conc, typ_dose,
    intervals = data.frame(start = 0, end = Inf, cl.obs = TRUE)
  ))
)$result |>
  as.data.frame() |>
  dplyr::filter(PPTESTCD == "cl.obs") |>
  dplyr::transmute(cohort, cl_total_sim = PPORRES)

cl_check <- closed_form |>
  dplyr::select(cohort, fu, LBM, cl_total_pred, vss_total_pred) |>
  dplyr::left_join(cl_obs_typ, by = "cohort") |>
  dplyr::mutate(pct_diff = 100 * (cl_total_sim - cl_total_pred) / cl_total_pred)

cl_check |>
  dplyr::transmute(
    Cohort = cohort,
    `fu` = fu,
    `LBM (kg)` = round(LBM, 1),
    `CL closed form (L/h)` = round(cl_total_pred, 3),
    `CL from solve (L/h)` = round(cl_total_sim, 3),
    `% diff` = round(pct_diff, 3)
  ) |>
  knitr::kable(
    caption = "Total-concentration clearance: closed form vs deterministic solve."
  )
Total-concentration clearance: closed form vs deterministic solve.
Cohort fu LBM (kg) CL closed form (L/h) CL from solve (L/h) % diff
CF female 0.743590 38.8 20.471 20.467 -0.022
CF male 0.562527 40.4 15.963 15.960 -0.022
HV female 0.544842 44.6 16.652 16.648 -0.021
HV male 0.500000 62.8 19.753 19.749 -0.018 

# Pure numerical agreement -- both sides use the same published numbers.
stopifnot(all(abs(cl_check$pct_diff) < 0.5))

Check 2: unbound PK is identical across all four groups at equal LBM

This is the paper’s central claim (Discussion: “the unbound clearances and unbound volumes of distribution were the same in all subject groups when subjects had the same body size”), and it is what makes the simulated PK/PD breakpoints in Table 5 identical between female and male CF patients.

eq_cov <- arm_spec |>
  dplyr::mutate(
    id = dplyr::row_number(),
    LBM = 53, # every arm at the reference LBM
    DOSE_CEFOTIAM_MG = dose_amt_mg
  ) |>
  dplyr::select(id, cohort, DIS_CF, SEXF, LBM, DOSE_CEFOTIAM_MG)

eq_events <- dplyr::bind_rows(
  eq_cov |>
    dplyr::mutate(
      time = 0, evid = 1L, amt = dose_amt_mg,
      cmt = "central", rate = dose_amt_mg / infusion_dur_h
    ),
  eq_cov |>
    tidyr::crossing(time = typ_times) |>
    dplyr::mutate(evid = 0L, amt = NA_real_, cmt = "Cc", rate = 0)
) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

eq <- rxode2::rxSolve(
  mod, eq_events,
  keep = c("cohort", "SEXF", "DIS_CF"),
  omega = NA, sigma = NA, useLinCmt = FALSE
) |>
  as.data.frame() |>
  dplyr::distinct(id, time, .keep_all = TRUE)

unbound_spread <- eq |>
  dplyr::filter(!is.na(Cunbound), Cunbound > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    rel_range = (max(Cunbound) - min(Cunbound)) / stats::median(Cunbound),
    .groups = "drop"
  )

total_spread <- eq |>
  dplyr::filter(!is.na(Cc), Cc > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    rel_range = (max(Cc) - min(Cc)) / stats::median(Cc),
    .groups = "drop"
  )

cat(sprintf(
  "Max relative spread across the four arms at equal LBM:\n  unbound Cc: %.3g\n  total   Cc: %.3g\n",
  max(unbound_spread$rel_range), max(total_spread$rel_range)
))
#> Max relative spread across the four arms at equal LBM:
#>   unbound Cc: 0
#>   total   Cc: 0.363

# Unbound profiles must superimpose exactly; total profiles must not.
stopifnot(max(unbound_spread$rel_range) < 1e-8)
stopifnot(max(total_spread$rel_range) > 0.3)

Check 3: fraction excreted unchanged is CLr,u / (CLr,u + CLnr,u)

Both clearance arms carry the same fu and the same allometric exponent, so their ratio is invariant to cohort, sex and body size. The model therefore predicts a single urinary recovery for everyone, which is a sharp, fully deterministic prediction to hold against the observed data.

fe_pred <- 100 * cl_ru / (cl_ru + cl_nru)

fe_sim <- typ |>
  dplyr::filter(!is.na(urinePct)) |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(fe_48h = max(urinePct), .groups = "drop")

cat(sprintf("Closed-form asymptotic recovery: %.2f%% of dose\n", fe_pred))
#> Closed-form asymptotic recovery: 68.39% of dose
print(as.data.frame(fe_sim))
#>      cohort   fe_48h
#> 1 CF female 68.39392
#> 2   CF male 68.39392
#> 3 HV female 68.39392
#> 4   HV male 68.39392

stopifnot(all(abs(fe_sim$fe_48h - fe_pred) / fe_pred < 0.01))

Paper Table 2 reports observed median recoveries of 70.3% (CF, range 47.6-77.8%) and 66.3% (HV, range 59.4-72.7%). The model’s single value of 68.4% sits between the two cohort medians and inside both observed ranges, which is the correct behaviour for a model that deliberately shares both clearance arms across cohorts.

Replicating the published figures

Figure 3: visual predictive check by cohort 

vpc <- sim_unique |>
  dplyr::filter(!is.na(Cc), Cc > 0) |>
  dplyr::mutate(group = ifelse(DIS_CF == 1L, "Patients with CF", "Healthy volunteers")) |>
  dplyr::group_by(group, time) |>
  dplyr::summarise(
    p10 = stats::quantile(Cc, 0.10),
    p25 = stats::quantile(Cc, 0.25),
    p50 = stats::quantile(Cc, 0.50),
    p75 = stats::quantile(Cc, 0.75),
    p90 = stats::quantile(Cc, 0.90),
    .groups = "drop"
  )

ggplot2::ggplot(vpc, ggplot2::aes(x = time)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = p10, ymax = p90), alpha = 0.18) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = p25, ymax = p75), alpha = 0.30) +
  ggplot2::geom_line(ggplot2::aes(y = p50), linewidth = 0.7) +
  ggplot2::facet_wrap(~group) +
  ggplot2::scale_y_log10() +
  ggplot2::labs(
    x = "Time (h)", y = "Total cefotiam in plasma (mg/L)",
    title = "Replicates Figure 3 of Shah 2019",
    subtitle = "Median with interquartile range and 80% prediction interval"
  ) +
  ggplot2::theme_minimal()

The protein-binding mechanism, total vs unbound

Figure 4 of the paper shows that the probability of target attainment is near-identical between female and male CF patients, because target attainment is driven by unbound concentration. The panel below is the direct visualisation of that result: at equal LBM the four arms separate on total concentration and superimpose completely on unbound concentration.

eq |>
  dplyr::filter(time > 0, time <= 8, !is.na(Cc)) |>
  dplyr::select(cohort, time, Total = Cc, Unbound = Cunbound) |>
  tidyr::pivot_longer(c(Total, Unbound), names_to = "scale", values_to = "conc") |>
  ggplot2::ggplot(ggplot2::aes(time, conc, colour = cohort)) +
  ggplot2::geom_line(linewidth = 0.8) +
  ggplot2::facet_wrap(~scale) +
  ggplot2::scale_y_log10() +
  ggplot2::labs(
    x = "Time (h)", y = "Cefotiam in plasma (mg/L)", colour = NULL,
    title = "All four groups at LBM = 53 kg",
    subtitle = "Total concentrations separate by group; unbound concentrations superimpose"
  ) +
  ggplot2::theme_minimal()

PKNCA validation

The paper’s Table 2 reports noncompartmental parameters from the observed data, stratified by cohort and by sex. Those are the reference values. Simulated profiles are censored below 1 mg/L before the NCA: the paper does not state an LLOQ, but 1 mg/L is the lowest concentration for which it reports assay recovery (Section 2.5), and an uncensored simulated profile lets lambda.z run down into concentrations no assay could have measured, which biases the terminal half-life upward.

lloq_mgL <- 1.0

sim_nca <- sim_unique |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(Cc = ifelse(time > 0 & Cc < lloq_mgL, NA_real_, Cc)) |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort)

# Defensive time-zero record so PKNCA never reports an AUC range starting
# before the first measurement.
sim_nca <- sim_nca |>
  dplyr::bind_rows(
    sim_nca |>
      dplyr::distinct(id, cohort) |>
      dplyr::mutate(time = 0, Cc = 0)
  ) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::arrange(id, time)

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

dose_df <- cohorts |>
  dplyr::mutate(time = 0, amt = dose_amt_mg) |>
  dplyr::select(id, cohort, time, amt)

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

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, aucinf.obs = TRUE,
  half.life = TRUE, cl.obs = TRUE, mrt.obs = TRUE, vss.obs = TRUE
)

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

nca_long <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c(
    "cmax", "tmax", "aucinf.obs", "half.life", "cl.obs", "mrt.obs", "vss.obs"
  ))

Comparison against the published noncompartmental analysis

Paper Table 2 reports medians by cohort and by sex within cohort. The reference rows below are the per-sex medians, which is the finest stratification the paper publishes and the one that exercises the fu split.

# Shah 2019 Table 2, per-cohort per-sex medians.
reference_nca <- tibble::tribble(
  ~cohort,     ~cl.obs, ~vss.obs, ~cmax, ~half.life, ~mrt.obs,
  "CF female", 22.1,    13.3,     NA,    NA,         NA,
  "CF male",   15.9,    12.3,     NA,    NA,         NA,
  "HV female", 16.2,    10.7,     NA,    NA,         NA,
  "HV male",   19.1,    16.7,     NA,    NA,         NA
) |>
  # Cmax, half-life and MRT are reported by cohort only (not split by sex),
  # so the cohort-level value is carried onto both sexes of that cohort.
  dplyr::mutate(
    cmax = c(124, 124, 111, 111),
    half.life = c(0.931, 0.931, 1.08, 1.08),
    mrt.obs = c(0.699, 0.699, 0.707, 0.707)
  )

nca_tbl <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_long,
  reference = reference_nca,
  by = "cohort",
  params = c("cmax", "cl.obs", "vss.obs", "half.life", "mrt.obs"),
  units = c(
    cmax = "mg/L", cl.obs = "L/h", vss.obs = "L",
    half.life = "h", mrt.obs = "h"
  )
)

knitr::kable(
  nca_tbl,
  caption = "Simulated vs published (Shah 2019 Table 2) noncompartmental parameters."
)
Simulated vs published (Shah 2019 Table 2) noncompartmental parameters.
NCA parameter cohort Reference Simulated % diff
Cmax (mg/L) CF female 124 111 -10.7%
Cmax (mg/L) CF male 124 135 +8.7%
Cmax (mg/L) HV female 111 129 +16.3%
Cmax (mg/L) HV male 111 108 -2.8%
t½ (h) CF female 0.931 0.687 -26.3%*
t½ (h) CF male 0.931 0.654 -29.8%*
t½ (h) HV female 1.08 0.778 -28.0%*
t½ (h) HV male 1.08 0.716 -33.7%*
CL/F (L/h) CF female 22.1 20.6 -6.8%
CL/F (L/h) CF male 15.9 17.9 +12.5%
CL/F (L/h) HV female 16.2 17.1 +5.7%
CL/F (L/h) HV male 19.1 20.3 +6.4%
Vss/F (L) CF female 13.3 13.8 +3.9%
Vss/F (L) CF male 12.3 11.2 -8.9%
Vss/F (L) HV female 10.7 11.9 +10.9%
Vss/F (L) HV male 16.7 13.8 -17.2%
MRT (h) CF female 0.699 0.651 -6.9%
MRT (h) CF male 0.699 0.621 -11.2%
MRT (h) HV female 0.707 0.705 -0.3%
MRT (h) HV male 0.707 0.687 -2.8% 
attr(nca_tbl, "footnote")
#> [1] "* differs from reference by more than ±20%."

The % diff column of the rendered table is a formatted string (it carries the * tolerance flag), so the assertions below recompute the differences numerically from the same medians the table is built from.

sim_med <- nca_long |>
  dplyr::group_by(cohort, PPTESTCD) |>
  dplyr::summarise(simulated = stats::median(PPORRES, na.rm = TRUE), .groups = "drop")

ref_long <- reference_nca |>
  tidyr::pivot_longer(-cohort, names_to = "PPTESTCD", values_to = "reference")

diffs <- dplyr::inner_join(sim_med, ref_long, by = c("cohort", "PPTESTCD")) |>
  dplyr::filter(!is.na(simulated), !is.na(reference)) |>
  dplyr::mutate(pct_diff = 100 * (simulated - reference) / reference)

knitr::kable(
  diffs |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), \(x) round(x, 2))),
  caption = "Numeric simulated-vs-published differences used by the assertions below."
)
Numeric simulated-vs-published differences used by the assertions below.
cohort PPTESTCD simulated reference pct_diff
CF female cl.obs 20.61 22.10 -6.76
CF female cmax 110.68 124.00 -10.74
CF female half.life 0.69 0.93 -26.26
CF female mrt.obs 0.65 0.70 -6.86
CF female vss.obs 13.81 13.30 3.86
CF male cl.obs 17.88 15.90 12.48
CF male cmax 134.77 124.00 8.68
CF male half.life 0.65 0.93 -29.76
CF male mrt.obs 0.62 0.70 -11.19
CF male vss.obs 11.20 12.30 -8.90
HV female cl.obs 17.12 16.20 5.71
HV female cmax 129.05 111.00 16.26
HV female half.life 0.78 1.08 -27.97
HV female mrt.obs 0.70 0.71 -0.35
HV female vss.obs 11.86 10.70 10.86
HV male cl.obs 20.33 19.10 6.42
HV male cmax 107.95 111.00 -2.75
HV male half.life 0.72 1.08 -33.73
HV male mrt.obs 0.69 0.71 -2.76
HV male vss.obs 13.82 16.70 -17.22 

# Terminal half-life is deliberately excluded from the pooled gate; it is not
# a property of the model alone but of the assay's quantification limit, and
# the paper publishes no LLOQ. The section below quantifies that and gates it
# separately on a quantity that does not depend on an assumed LLOQ.
pct <- abs(diffs$pct_diff[diffs$PPTESTCD != "half.life"])

# Structural gate: a mis-transcribed clearance, dose or unit would move the
# whole distribution by tens of percent. Assert on the CENTRE and on a robust
# quantile, not on the extreme -- the extreme of a random cohort is not
# reproducible across rxode2 builds.
stopifnot(stats::median(pct) < 12)
stopifnot(stats::quantile(pct, 0.75) < 20)

# Clearance is the parameter the model is really making a claim about, and the
# one Table 2 reports per sex. Hold its centre tighter.
cl_pct <- abs(diffs$pct_diff[diffs$PPTESTCD == "cl.obs"])
stopifnot(stats::median(cl_pct) < 12)

Terminal half-life is window-limited, not mis-specified

Every simulated half-life above is 26-34% below the published value, uniformly across all four arms. That uniformity is the clue: a genuine parameter error would not bias all four arms by the same fraction while leaving clearance, volume and MRT within 12%. The cause is that a noncompartmental terminal half-life is not a property of the model at all – it is a property of how far down the concentration profile the assay can still see, and Shah 2019 does not publish a quantification limit.

The model’s own terminal half-life is unambiguous: it is log(2) over the slowest eigenvalue of the three-compartment rate matrix, and because fu cancels out of the ODE system entirely it depends only on LBM.

kel <- (cl_ru + cl_nru) / v1u
k12 <- 13.8043 / v1u
k21 <- 13.8043 / v2u
k13 <- 1.83622 / v1u
k31 <- 1.83622 / v3u

rate_matrix <- matrix(
  c(
    -(kel + k12 + k13), k21, k31,
    k12, -k21, 0,
    k13, 0, -k31
  ),
  nrow = 3, byrow = TRUE
)

lambda_z_ref <- -max(Re(eigen(rate_matrix)$values))
thalf_ref <- log(2) / lambda_z_ref

cat(sprintf(
  "Model terminal half-life at the reference LBM of 53 kg: %.2f h\n", thalf_ref
))
#> Model terminal half-life at the reference LBM of 53 kg: 1.85 h

# Terminal disposition scales as LBM^(1 - 0.75) = LBM^0.25.
thalf_by_arm <- arm_spec |>
  dplyr::transmute(
    cohort,
    `LBM (kg)` = lbm_med,
    `Model terminal t1/2 (h)` = round(thalf_ref * (lbm_med / 53)^0.25, 2)
  )
knitr::kable(thalf_by_arm, caption = "Model terminal half-life by arm.")
Model terminal half-life by arm.
cohort LBM (kg) Model terminal t1/2 (h)
CF female 38.8 1.71
CF male 40.4 1.72
HV female 44.6 1.77
HV male 62.8 1.93

Paper Table 2 publishes terminal half-life as a median with a range: 0.931 [0.881-1.91] h in CF patients and 1.08 [0.753-1.66] h in healthy volunteers. These are not estimates of the model’s asymptotic terminal slope – they are window-limited NCA estimates on real, LLOQ-truncated data, and truncation can only ever shorten an estimated terminal half-life. So the correct expectation is that the model’s asymptotic value exceeds both published medians, which it does. It lands inside the CF observed range and slightly above the narrower HV observed upper bound of 1.66 h, which is the expected direction rather than a defect.

# Window truncation can only shorten an NCA half-life, so the model's
# asymptotic value must exceed both published medians.
stopifnot(thalf_ref > 0.931, thalf_ref > 1.08)
# Order-of-magnitude sanity only -- deliberately loose, not a tuned bound.
stopifnot(thalf_ref < 3)

Sweeping the censoring threshold shows exactly how the NCA estimate walks from the model’s true terminal slope up to the steep distribution phase as the assumed quantification limit rises. The published medians sit between the 0.1 and 1 mg/L curves, implying the study’s real quantification limit was somewhere near 0.3-0.5 mg/L.

halflife_at_lloq <- function(profile, dose_frame, lloq) {
  d <- profile |>
    dplyr::mutate(Cc = ifelse(time > 0 & Cc < lloq, NA_real_, Cc)) |>
    dplyr::filter(!is.na(Cc)) |>
    dplyr::select(id, time, Cc, cohort)
  res <- suppressWarnings(PKNCA::pk.nca(PKNCA::PKNCAdata(
    PKNCA::PKNCAconc(d, Cc ~ time | cohort + id),
    PKNCA::PKNCAdose(dose_frame, amt ~ time | cohort + id),
    intervals = data.frame(start = 0, end = Inf, half.life = TRUE)
  )))
  as.data.frame(res$result) |>
    dplyr::filter(PPTESTCD == "half.life") |>
    dplyr::transmute(cohort, thalf = PPORRES)
}

typ_profile <- typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort)
typ_dose_frame <- typ_cov |>
  dplyr::mutate(time = 0, amt = dose_amt_mg) |>
  dplyr::select(id, cohort, time, amt)

sens <- lapply(
  c(1, 0.1, 0.01),
  \(l) halflife_at_lloq(typ_profile, typ_dose_frame, l) |> dplyr::mutate(lloq = l)
) |>
  dplyr::bind_rows() |>
  dplyr::mutate(thalf = round(thalf, 2)) |>
  tidyr::pivot_wider(names_from = lloq, values_from = thalf, names_prefix = "LLOQ ") |>
  dplyr::rename(Cohort = cohort)

knitr::kable(
  sens,
  caption = paste(
    "Terminal half-life (h) from NCA on the deterministic profile as a",
    "function of the assumed quantification limit (mg/L)."
  )
)
Terminal half-life (h) from NCA on the deterministic profile as a function of the assumed quantification limit (mg/L).
Cohort LLOQ 1 LLOQ 0.1 LLOQ 0.01
CF female 0.78 1.51 1.68
CF male 0.84 1.59 1.70
HV female 0.84 1.63 1.74
HV male 0.86 1.69 1.90 

# The published medians must be reachable within the swept window -- i.e. the
# 1 mg/L censoring used above is too aggressive and a lower LLOQ recovers the
# published value. No LLOQ was tuned to make this true.
stopifnot(min(sens$`LLOQ 1`) < 0.931, max(sens$`LLOQ 0.01`) > 1.08)

The 1 mg/L threshold used in the main NCA is retained rather than tuned downward: it is the only concentration for which the paper reports assay recovery, and clearance, volume and MRT – which are dominated by the early, high-concentration part of the profile – are insensitive to it.

The clearance rows are the substantive result. The model reproduces the sign and ordering of the Table 2 medians that motivated the paper: within the CF cohort female clearance exceeds male clearance despite near-identical LBM, whereas within the HV cohort male clearance exceeds female clearance because male HVs are 41% leaner-heavier. No covariate other than fu and LBM is doing any work.

cl_sim <- nca_long |>
  dplyr::filter(PPTESTCD == "cl.obs") |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(cl = stats::median(PPORRES), .groups = "drop")
print(as.data.frame(cl_sim))
#>      cohort       cl
#> 1 CF female 20.60530
#> 2   CF male 17.88466
#> 3 HV female 17.12495
#> 4   HV male 20.32585

cl_of <- function(x) cl_sim$cl[cl_sim$cohort == x]
stopifnot(cl_of("CF female") > cl_of("CF male")) # driven purely by fu
stopifnot(cl_of("HV male") > cl_of("HV female")) # driven purely by LBM

Urinary recovery against Table 2 

urine_obs <- sim_unique |>
  dplyr::filter(!is.na(urinePct), time == 24) |>
  dplyr::mutate(group = ifelse(DIS_CF == 1L, "Patients with CF", "Healthy volunteers")) |>
  dplyr::group_by(group) |>
  dplyr::summarise(
    `Simulated median (%)` = round(stats::median(urinePct), 1),
    `Simulated 10th-90th (%)` = sprintf(
      "%.1f-%.1f",
      stats::quantile(urinePct, 0.10), stats::quantile(urinePct, 0.90)
    ),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    `Observed median (%)` = c(66.3, 70.3),
    `Observed range (%)` = c("59.4-72.7", "47.6-77.8")
  )

knitr::kable(
  urine_obs,
  caption = "Cumulative 24 h urinary recovery of unchanged cefotiam vs Shah 2019 Table 2."
)
Cumulative 24 h urinary recovery of unchanged cefotiam vs Shah 2019 Table 2.
group Simulated median (%) Simulated 10th-90th (%) Observed median (%) Observed range (%)
Healthy volunteers 69.8 58.2-76.4 66.3 59.4-72.7
Patients with CF 68.8 58.7-76.6 70.3 47.6-77.8 

# The model shares both clearance arms across cohorts, so it predicts one
# recovery for everyone; assert only that it lands inside both observed ranges.
stopifnot(all(urine_obs$`Simulated median (%)` > 59.4))
stopifnot(all(urine_obs$`Simulated median (%)` < 77.8))

Assumptions and deviations

  • Only the final model is encoded. The Supplementary Materials also report a disease-scale-factor model (Table S1) and a glomerular-filtration / tubular-secretion model (Table S2). Neither is presented by the paper as final, and Table S2 fitted worse (-2x log-likelihood worse by 7.1), so both are excluded as model-development iterations.
  • Six-figure parameter values come from the authors’ own code. Where the Berkeley Madonna listing in the Supplementary Materials prints more significant figures than Table 4, the listing is used. Every such value is quoted against its Table 4 rounding in the model file comments; none disagrees with Table 4 beyond rounding.
  • BSV is a log-scale standard deviation. Table 4 footnote a calls the BSV an “apparent coefficient of variation of a normal distribution on natural logarithmic scale”, and Section 2.6.4 says eta had “standard deviation BSV”. The model therefore uses omega^2 = BSV^2, not log(1 + CV^2). The supplement’s normal(0, CV_CLR) draws confirm this directly.
  • No LLOQ is published, and the terminal half-life comparison is therefore not well posed. Simulated profiles are censored at 1 mg/L for the NCA only, that being the lowest concentration at which Section 2.5 reports assay recovery; the model itself has no censoring. This makes every simulated terminal half-life 26-34% lower than Table 2, uniformly across all four arms. The sensitivity sweep above shows the NCA estimate moving from 0.68 h at a 1 mg/L limit to 1.70 h uncensored, so the published 0.931 / 1.08 h medians imply a real quantification limit near 0.3-0.5 mg/L. The censoring threshold was deliberately not lowered to close the gap, because choosing an LLOQ to match a validation target is tuning. Half-life is instead excluded from the pooled numeric gate and checked separately against the published ranges, which needs no LLOQ assumption. Clearance, volume and MRT are dominated by the early profile and are insensitive to the threshold.
  • The urinary endpoint needs a dose column. The paper fitted the fraction of dose in urine and reports its residual SD on a percent scale (0.384%). To keep that published number attached to the paper’s own observation scale, the model expresses urinePct = 100 * urine / DOSE_CEFOTIAM_MG. A user who prefers a cumulative amount can read the urine state directly, in mg.
  • Cohort LBM distributions are reconstructed, not published. Table 1 gives per-sex medians and ranges but no distributional form. The cohort here is drawn log-normally with the 15% CV the paper used in its own Monte Carlo covariate model (Section 2.6.7) and truncated to the Table 1 ranges. The HV female range in Table 1 is extremely narrow (44.6-45.4 kg, n = 3), so that arm is nearly deterministic.
  • Table 2 reports Cmax, half-life and MRT by cohort only, not by sex, so the cohort value is carried onto both sexes of that cohort in the comparison table. Only clearance and Vss are genuinely sex-stratified in Table 2.
  • No albumin covariate. The Supplementary Materials back-calculate an albumin ratio consistent with the estimated fu values, but the study did not measure albumin and the paper states the effect of albumin “was not included in the model”. No albumin column is carried here.
  • fu random effects. The three estimated unbound fractions carry a BSV the authors fixed at 5% CV (Table 4 footnote d); the fixed HV-male fu carries no random effect, matching the supplement listing, which draws ETA_FU_CFF, ETA_FU_CFM and ETA_FU_HVF but no ETA_FU_HVM.