Ciprofol (Wang 2026) • nlmixr2lib

Model and source

Citation: Wang S, Li Y, Hu Z, Du L, Wang Y, Jiang X, Li L, Shangguan W. (2026). Population pharmacokinetics of a single bolus of ciprofol in Chinese pediatric patients. BMC Anesthesiology 26(1). doi:10.1186/s12871-026-03647-9. PMCID: PMC12930602. Chinese Clinical Trial Registry ChiCTR2200058405.
Description: Three-compartment intravenous population PK model for ciprofol (HSK3486) after a single 0.6 mg/kg bolus given over 30 s in Chinese pediatric surgical patients (Wang 2026; 25 children aged 1-9 years, ASA physical status I-II, scheduled for elective urologic surgery; 317 arterial plasma samples). A three-compartment model was significantly better than a two-compartment model (dOFV = 65.8, p < 0.001). All disposition parameters are reported per kilogram of body weight (CL 31.2 mL/min/kg, V1 506 mL/kg, Q2 28.2 mL/min/kg, V2 231 mL/kg, Q3 19.8 mL/min/kg, V3 1360 mL/kg), i.e. body weight enters every parameter with a linear (exponent 1) scaling; standard allometric scaling and age-dependent maturation functions on clearance were tested and did not improve the fit, and no further effect of weight, age, sex or BMI was detectable after per-kilogram normalisation. Blood urea nitrogen was the single retained covariate, acting on the central volume as a power of the ratio to the cohort median with an estimated exponent of -0.821, so that V1 falls from 0.770 L/kg at BUN 3 mmol/L to 0.384 L/kg at BUN 7 mmol/L; the authors judged the resulting exposure change clinically insignificant. Log-normal inter-individual variability was retained on CL, V1 (correlated, r = -0.821) and Q3 only, because the IIV estimates for V2, V3 and Q2 were close to zero. Residual error is combined proportional plus additive.
Article: https://doi.org/10.1186/s12871-026-03647-9 (BMC Anesthesiology 26(1), open access; PMCID PMC12930602)
Trial registration: Chinese Clinical Trial Registry ChiCTR2200058405

No supplementary material accompanies this article (EuropePMC reports hasSuppl = "N" for an open-access record, so the field is informative), and no erratum or correction notice was found.

Population

Wang 2026 enrolled 27 Chinese children aged 1 to 9 years with American Society of Anesthesiologists physical status I or II, scheduled for elective urologic surgery of anticipated duration greater than 2 h, at the Second Affiliated Hospital and Yuying Children’s Hospital of Wenzhou Medical University between January and August 2023. Enrolment was balanced across three age strata (toddlers 1-2 years, preschoolers 3-5 years, school-age children 6-9 years, n = 9 each). One child withdrew when the operation was shortened and one was excluded for a blocked arterial catheter, leaving 25 children (14 male, 11 female) in strata of 9 / 9 / 7. Inclusion required a body mass index between the 25th and 75th percentile for age and sex, so the cohort is deliberately non-obese.

Each child received a single 0.6 mg/kg intravenous bolus of ciprofol over 30 s at induction, after midazolam 0.1-0.2 mg/kg premedication and together with fentanyl 2.0 ug/kg. Thirteen arterial samples were planned per patient (a pre-dose blank plus 2, 4, 6, 8, 10, 20, 30, 45, 60, 90, 120 and 180 min after injection); 8 of 325 samples were not collected, mostly the 180 min late-elimination sample, leaving 317 samples for the analysis. Ciprofol was quantified by UPLC-APCI-MS/MS over 5-20000 ng/mL with an LLOQ of 5 ng/mL, and no concentration fell below the limit of quantitation. Baseline demographics and laboratory values are in Table 1 of the source: mean age 4.2 years (range 1.0-9.0), mean weight 18.5 kg (range 9.8-37), mean BMI 16.7 kg/m^2, mean BUN 5.2 mmol/L (range 2.9-7.6).

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

Source trace

The per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Wang_2026_ciprofol.R. The table below collects them in one place for review. The paper reports every disposition parameter per kilogram of body weight and in millilitres, so each ini() value is the printed number divided by 1000 (mL to L); body weight is reapplied in model() with a linear exponent.

Equation / parameter	Paper value	`ini()` value	Source location
`lcl` (CL)	31.2 mL/min/kg (RSE 5.0%)	`log(0.0312)`	Table 3, row `CL, ml/min/kg`
`lvc` (V1)	506 mL/kg (RSE 8.1%)	`log(0.506)`	Table 3, row `V1, ml/kg`
`lq` (Q2)	28.2 mL/min/kg (RSE 19.9%)	`log(0.0282)`	Table 3, row `Q2, ml/min/kg`
`lvp` (V2)	231 mL/kg (RSE 14.7%)	`log(0.231)`	Table 3, row `V2, ml/kg`
`lq2` (Q3)	19.8 mL/min/kg (RSE 10.8%)	`log(0.0198)`	Table 3, row `Q3, ml/min/kg`
`lvp2` (V3)	1360 mL/kg (RSE 12.4%)	`log(1.36)`	Table 3, row `V3, ml/kg`
`e_bun_vc`	-0.821 (RSE 28.6%)	`-0.821`	Table 3, row `Covariate thetaBUN on V1`
BUN normalising median	not printed; recovered as 5 mmol/L	`(BUN / 5)`	Results, “Population covariant analysis” (V1 = 0.770 and 0.384 L/kg at BUN 3 and 7); Methods Eq. 5
`etalcl`	eta(CL) 26.3% (RSE 9.5%, shrinkage 0.3%)	`0.263^2 = 0.069169`	Table 3, `Inter-individual variability` block
`etalvc`	eta(V1) 35.1% (RSE 11.1%, shrinkage 5.8%)	`0.351^2 = 0.123201`	Table 3, `Inter-individual variability` block
cov(`etalcl`, `etalvc`)	Corr_CL & V1 = -0.821	`-0.821 * 0.263 * 0.351 = -0.075789`	Table 3, row `Corr_CL & V1`; Results, “Population PK model”
`etalq2`	eta(Q3) 38.3% (RSE 15.5%, shrinkage 7.8%)	`0.383^2 = 0.146689`	Table 3, `Inter-individual variability` block
`propSd`	13.5% (RSE 9.6%)	`0.135`	Table 3, row `epsilon_prop`
`addSd`	4.83 ug/L (RSE 28.8%)	`0.00483` (mg/L)	Table 3, row `epsilon_add`
Three-compartment disposition, IV	n/a	`d/dt(central)`, `d/dt(peripheral1)`, `d/dt(peripheral2)`	Results, “Population PK model” (dOFV = 65.8 vs two-compartment, p < 0.001)
Linear body-weight scaling of all disposition parameters	n/a	`* WT` on `cl`, `vc`, `q`, `vp`, `q2`, `vp2`	Table 3 per-kilogram units; Methods, “Model selection and covariate analysis” (allometric scaling and maturation tested, not retained)
Log-normal IIV	n/a	`exp(lX + etalX)`	Methods, “Structural model” (Eq. 1)
Combined residual error	n/a	`Cc ~ add(addSd) + prop(propSd)`	Methods, Eq. 4

Virtual cohort

Original observed data are not publicly available. The cohort below reproduces the three enrolment strata of Table 1, drawing weight, BUN and age per stratum from normal distributions centred on the published stratum mean with a standard deviation of (max - min) / 4 and truncated to the published range. 100 subjects per stratum (the skill cap is 200 per arm) is ample for the checks below.

Each subject receives the study regimen: 0.6 mg/kg as a 30 s infusion into central, with observations on the ODE state central at the paper’s own sampling times. Using the published sparse grid rather than a dense one matters for the NCA comparison: the trapezoidal bias and the terminal-slope window then match those of the paper’s own non-compartmental analysis.

# set.seed() seeds R's RNG (the covariate draws); rxode2's simulation RNG is
# partitioned per solver thread, so the realised cohort differs between a
# 2-thread CI runner and a 16-thread workstation. Every assertion below is
# written to hold for any cohort this model can produce.
set.seed(20260912)
rxode2::rxSetSeed(20260912)

SAMPLE_TIMES <- c(2, 4, 6, 8, 10, 20, 30, 45, 60, 90, 120, 180)
LLOQ <- 0.005  # 5 ng/mL = 5 ug/L = 0.005 mg/L (Methods, assay validation)

rtnorm <- function(n, mean, lo, hi) pmin(pmax(rnorm(n, mean, (hi - lo) / 4), lo), hi)

make_cohort <- function(n, label, wt, bun, age, id_offset = 0L) {
  subj <- tibble(
    id      = id_offset + seq_len(n),
    stratum = label,
    WT      = rtnorm(n, wt[1],  wt[2],  wt[3]),
    BUN     = rtnorm(n, bun[1], bun[2], bun[3]),
    AGE     = rtnorm(n, age[1], age[2], age[3])
  )
  bind_rows(
    # 0.6 mg/kg given over 30 s (0.5 min) into the central compartment
    subj |> mutate(time = 0, amt = 0.6 * WT, evid = 1L, dur = 0.5, cmt = "central"),
    # Observations on the ODE state `central` -- never on the observable `Cc`.
    # time = 0 is the paper's pre-dose blank sample and anchors PKNCA's AUC.
    tidyr::crossing(subj, time = c(0, SAMPLE_TIMES)) |>
      mutate(amt = NA_real_, evid = 0L, dur = NA_real_, cmt = "central")
  ) |>
    arrange(id, time, desc(evid))
}

events <- bind_rows(
  make_cohort(100, "Toddlers (1-2 y)",     c(12.8,  9.8, 16.5), c(5.7, 4.0, 6.8), c(1.8, 1, 2),   0L),
  make_cohort(100, "Preschoolers (3-5 y)", c(18.3, 15.0, 23.5), c(5.1, 2.9, 7.6), c(3.9, 3, 5), 100L),
  make_cohort(100, "School-age (6-9 y)",   c(25.9, 19.5, 37.0), c(4.7, 3.8, 6.2), c(7.7, 6, 9), 200L)
)

# Disjoint ids across the three cohorts (duplicate ids silently merge in rxSolve)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
stopifnot(nrow(distinct(events, id)) == 300L)

Simulation

mod <- readModelDb("Wang_2026_ciprofol")

sim <- rxode2::rxSolve(mod, events = events, keep = c("stratum", "WT", "BUN")) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'

# Cc is the individual prediction (mg/L); `sim` carries the residual error and
# is the analogue of an observed concentration. The paper reports ug/L.
sim <- sim |>
  mutate(
    Cc_ugL  = Cc * 1000,
    obs_ugL = pmax(sim, LLOQ / 2) * 1000  # see "Assumptions and deviations"
  )

Replicate published figures

sim |>
  filter(!is.na(Cc), time > 0) |>
  group_by(stratum, time) |>
  summarise(
    Q025 = quantile(obs_ugL, 0.025),
    Q50  = quantile(obs_ugL, 0.500),
    Q975 = quantile(obs_ugL, 0.975),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50, colour = stratum, fill = stratum)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975), alpha = 0.15, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(
    x = "Time after ciprofol injection (min)",
    y = "Ciprofol plasma concentration (ug/L)",
    colour = NULL, fill = NULL,
    title = "Figure 2 -- concentration-time profiles by age stratum",
    caption = "Median and 2.5th-97.5th percentiles of 100 simulated subjects per stratum."
  ) +
  theme(legend.position = "bottom")

Replicates Figure 2 of Wang 2026: ciprofol concentration-time profiles by age stratum.

Wang 2026 reports that “the plasma concentration-time curve of ciprofol was similar across the three groups” and that Table 2 showed no significant difference between strata (all p > 0.15). The simulated profiles overlap for the same reason the published ones do: after per-kilogram dosing, weight and age cancel out of the model entirely, and the only covariate that moves a prediction is BUN.

Typical-value replication of Table 4

Table 4 of Wang 2026 reports the paper’s own simulation of a typical child given 0.6 mg/kg over 30 s at three BUN levels. This is a deterministic, typical-value prediction, so it is reproduced with zeroRe() and gated tightly. It is the sharpest available check on the covariate model, because the central volume at BUN 3 and 7 mmol/L is quoted in the Results text and pins both the exponent and the (unprinted) normalising median.

mod_typical <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
WT_TYP <- 18.5  # cohort mean weight, Table 1; predictions are weight-invariant

typical_profile <- function(bun) {
  ev <- rxode2::et(amt = 0.6 * WT_TYP, dur = 0.5, cmt = "central") |>
    rxode2::et(seq(0, 180, by = 0.05), cmt = "central")
  d <- as.data.frame(ev)
  d$WT <- WT_TYP
  d$BUN <- bun
  rxode2::rxSolve(mod_typical, d, returnType = "data.frame") |>
    filter(!is.na(Cc)) |>
    mutate(BUN = bun, Cc_ugL = Cc * 1000)
}

profiles <- bind_rows(lapply(c(3, 5, 7), typical_profile))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq2'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq2'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq2'

at_time <- function(df, t) df$Cc_ugL[which.min(abs(df$time - t))]

table4 <- profiles |>
  group_by(BUN) |>
  summarise(
    C2min   = at_time(pick(everything()), 2),
    C180min = at_time(pick(everything()), 180),
    # AUC(0-inf) = trapezoid to 180 min + Clast / lambda_z from the terminal slope
    AUC = {
      d <- pick(everything())
      ug <- d$Cc_ugL
      auc180 <- sum(diff(d$time) * (head(ug, -1) + tail(ug, -1)) / 2)
      keep <- d$time >= 120
      lz <- -stats::coef(stats::lm(log(ug[keep]) ~ d$time[keep]))[[2]]
      auc180 + tail(ug, 1) / lz
    },
    V1_L_per_kg = 0.506 * (first(BUN) / 5)^(-0.821),
    .groups = "drop"
  )

published4 <- tibble::tribble(
  ~BUN, ~C2min_pub, ~C180min_pub, ~AUC_pub, ~V1_pub,
     3,      640.7,         19.6,    19231,   0.770,
     5,      880.3,         17.9,    19231,      NA,
     7,     1056.1,         17.1,    19231,   0.384
)

cmp4 <- table4 |>
  left_join(published4, by = "BUN") |>
  mutate(
    `C2min % diff`   = 100 * (C2min - C2min_pub) / C2min_pub,
    `C180min % diff` = 100 * (C180min - C180min_pub) / C180min_pub,
    `AUC % diff`     = 100 * (AUC - AUC_pub) / AUC_pub,
    `V1 % diff`      = 100 * (V1_L_per_kg - V1_pub) / V1_pub
  )

cmp4 |>
  transmute(
    `BUN (mmol/L)`            = BUN,
    `V1 simulated (L/kg)`     = round(V1_L_per_kg, 3),
    `V1 published (L/kg)`     = V1_pub,
    `C2min simulated (ug/L)`  = round(C2min, 1),
    `C2min published (ug/L)`  = C2min_pub,
    `C180min simulated (ug/L)` = round(C180min, 1),
    `C180min published (ug/L)` = C180min_pub,
    `AUC simulated (ug*min/L)` = round(AUC, 0),
    `AUC published (ug*min/L)` = AUC_pub
  ) |>
  knitr::kable(
    caption = "Typical-value replication of Wang 2026 Table 4 (0.6 mg/kg over 30 s) and of the two central volumes quoted in the Results text."
  )

Typical-value replication of Wang 2026 Table 4 (0.6 mg/kg over 30 s) and of the two central volumes quoted in the Results text.
BUN (mmol/L)	V1 simulated (L/kg)	V1 published (L/kg)	C2min simulated (ug/L)	C2min published (ug/L)	C180min simulated (ug/L)	C180min published (ug/L)	AUC simulated (ug*min/L)	AUC published (ug*min/L)
3	0.770	0.770	656.0	640.7	19.7	19.6	19124	19231
5	0.506	NA	912.4	880.3	17.9	17.9	19194	19231
7	0.384	0.384	1107.0	1056.1	17.2	17.1	19211	19231

# Closed form: with a covariate-free clearance, AUC(0-inf) is exactly Dose / CL
# and is invariant to BUN -- which is why Table 4 prints the same AUC on all
# three rows. This is an analytic identity, so it is gated tightly.
auc_closed_form <- 0.6 / 0.0312 * 1000  # ug*min/L
stopifnot(abs(auc_closed_form - 19231) / 19231 < 0.001)

# Deterministic (zeroRe, fixed covariates) -- no cohort randomness, so tight
# bounds are correct here and will catch a mis-transcribed parameter.
stopifnot(
  # Recovers the two published central volumes exactly, confirming BOTH the
  # -0.821 exponent and the normalising median of 5 mmol/L.
  max(abs(cmp4$`V1 % diff`), na.rm = TRUE) < 0.5,
  # Numerical AUC vs the paper's; residual gap is trapezoid + extrapolation.
  max(abs(cmp4$`AUC % diff`)) < 2,
  max(abs(cmp4$`C180min % diff`)) < 2,
  # C2min sits between our t = 2.0 min and t = 2.5 min values: the paper does
  # not state whether "2 min" runs from the start or the end of the 30 s
  # injection, and the curve falls ~7% per 0.1 min there. Realised 2.4-4.8%.
  max(abs(cmp4$`C2min % diff`)) < 6
)

PKNCA validation

The paper’s Table 2 reports non-compartmental parameters computed from the observed concentrations of each child. The simulation is therefore taken through the same route: concentrations carrying residual error, on the published sampling grid, with the pre-dose blank as the time-zero anchor.

# Concentrations in the paper's units (ug/L). Filter on !is.na() ONLY -- a
# `time > 0` or `Cc > 0` filter would drop the time-zero anchor.
nca_conc <- sim |>
  filter(!is.na(Cc)) |>
  mutate(Cc = ifelse(time == 0, 0, obs_ugL)) |>   # pre-dose blank is exactly 0
  select(id, time, Cc, stratum)

# Duplicate every subject into a pooled "All" group so the per-stratum and the
# pooled comparison land in a single table (paper Table 2 reports both).
nca_conc <- bind_rows(
  nca_conc,
  nca_conc |> mutate(id = id + 1000L, stratum = "All (n = 25)")
)

nca_dose <- events |>
  filter(evid == 1) |>
  transmute(id, time, amt = amt * 1000, stratum)  # mg -> ug
nca_dose <- bind_rows(
  nca_dose,
  nca_dose |> mutate(id = id + 1000L, stratum = "All (n = 25)")
)

# Time-zero guarantee (already present above; kept as a defensive no-op).
nca_conc <- bind_rows(
  nca_conc,
  nca_conc |> distinct(id, stratum) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, stratum, time, .keep_all = TRUE) |>
  arrange(id, stratum, time)

conc_obj <- PKNCA::PKNCAconc(nca_conc, Cc ~ time | stratum + id,
                             concu = "ng/mL", timeu = "min")
dose_obj <- PKNCA::PKNCAdose(nca_dose, amt ~ time | stratum + id, doseu = "ug")

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

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

# cl.obs comes back as dose / AUC in L/min; the paper reports mL/min/kg.
wt_by_id <- distinct(bind_rows(
  events |> filter(evid == 1) |> select(id, WT),
  events |> filter(evid == 1) |> transmute(id = id + 1000L, WT)
))

nca_tbl <- as.data.frame(nca_res$result) |>
  left_join(wt_by_id, by = "id") |>
  mutate(PPORRES = ifelse(PPTESTCD == "cl.obs", PPORRES / WT * 1000, PPORRES)) |>
  select(-WT)

stopifnot(nrow(nca_tbl) > 0)

Comparison against published NCA

published <- tibble::tribble(
  ~stratum,                ~cmax,  ~tmax, ~auclast, ~aucinf.obs, ~half.life, ~cl.obs,
  "Toddlers (1-2 y)",      976.3,      2,    15248,       16401,       54.3,    39.1,
  "Preschoolers (3-5 y)", 1048.2,      2,    17029,       19424,       64.4,    31.8,
  "School-age (6-9 y)",    881.2,      2,    15966,       17493,       57.6,    35.4,
  "All (n = 25)",          975.5,      2,    16090,       17795,       58.8,    35.5
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_tbl,
  reference = published,
  by        = "stratum",
  units     = c(cmax = "ug/L", tmax = "min", auclast = "ug*min/L",
                aucinf.obs = "ug*min/L", half.life = "min", cl.obs = "mL/min/kg"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated vs. published non-compartmental parameters (Wang 2026 Table 2). * differs from the reference by more than 20%.",
  digits  = 1
)

Simulated vs. published non-compartmental parameters (Wang 2026 Table 2). * differs from the reference by more than 20%.
NCA parameter	stratum	Reference	Simulated	% diff
Cmax (ug/L)	Toddlers (1-2 y)	976	973	-0.3%
Cmax (ug/L)	Preschoolers (3-5 y)	1050	879	-16.1%
Cmax (ug/L)	School-age (6-9 y)	881	881	+0.0%
Cmax (ug/L)	All (n = 25)	976	897	-8.0%
Tmax (min)	Toddlers (1-2 y)	2	2	+0.0%
Tmax (min)	Preschoolers (3-5 y)	2	2	+0.0%
Tmax (min)	School-age (6-9 y)	2	2	+0.0%
Tmax (min)	All (n = 25)	2	2	+0.0%
AUC0-∞ (obs) (ug*min/L)	Toddlers (1-2 y)	16400	18300	+11.4%
AUC0-∞ (obs) (ug*min/L)	Preschoolers (3-5 y)	19400	19000	-2.2%
AUC0-∞ (obs) (ug*min/L)	School-age (6-9 y)	17500	18400	+5.3%
AUC0-∞ (obs) (ug*min/L)	All (n = 25)	17800	18600	+4.5%
AUClast (ug*min/L)	Toddlers (1-2 y)	15200	16500	+7.9%
AUClast (ug*min/L)	Preschoolers (3-5 y)	17000	17200	+0.9%
AUClast (ug*min/L)	School-age (6-9 y)	16000	16400	+3.0%
AUClast (ug*min/L)	All (n = 25)	16100	16500	+2.7%
t½ (min)	Toddlers (1-2 y)	54.3	57.1	+5.1%
t½ (min)	Preschoolers (3-5 y)	64.4	59.7	-7.2%
t½ (min)	School-age (6-9 y)	57.6	55.9	-3.0%
t½ (min)	All (n = 25)	58.8	57.4	-2.4%
CL/F (mL/min/kg)	Toddlers (1-2 y)	39.1	32.9	-16.0%
CL/F (mL/min/kg)	Preschoolers (3-5 y)	31.8	31.6	-0.7%
CL/F (mL/min/kg)	School-age (6-9 y)	35.4	32.6	-7.9%
CL/F (mL/min/kg)	All (n = 25)	35.5	32.3	-9.1%

# Recompute the pooled comparison from the PKNCA result frame rather than
# parsing the formatted display table. ncaComparisonTable() aggregates by
# median, so the gate does too and the two agree by construction.
GATED <- c("cmax", "tmax", "auclast", "aucinf.obs", "half.life", "cl.obs")

pooled <- nca_tbl |>
  filter(stratum == "All (n = 25)", PPTESTCD %in% GATED) |>
  group_by(PPTESTCD) |>
  summarise(Simulated = median(PPORRES, na.rm = TRUE),
            n_subj = sum(!is.na(PPORRES)), .groups = "drop") |>
  left_join(
    published |>
      filter(stratum == "All (n = 25)") |>
      select(-stratum) |>
      tidyr::pivot_longer(everything(), names_to = "PPTESTCD", values_to = "Reference"),
    by = "PPTESTCD"
  ) |>
  mutate(`% diff` = 100 * (Simulated - Reference) / Reference)

# A lookup that matches no rows returns numeric(0), and all(logical(0)) is
# TRUE -- so confirm every gated parameter actually produced results before
# asserting on them.
stopifnot(
  nrow(pooled) == length(GATED),
  setequal(pooled$PPTESTCD, GATED),
  !anyNA(pooled$Simulated), !anyNA(pooled$Reference),
  all(pooled$n_subj > 250)
)

knitr::kable(pooled, digits = 2,
             caption = "Pooled (n = 300 simulated) values used by the gate below.")

Pooled (n = 300 simulated) values used by the gate below.
PPTESTCD	Simulated	n_subj	Reference	% diff
aucinf.obs	18592.44	300	17795.0	4.48
auclast	16521.67	300	16090.0	2.68
cl.obs	32.27	300	35.5	-9.10
cmax	897.16	300	975.5	-8.03
half.life	57.37	300	58.8	-2.42
tmax	2.00	300	2.0	0.00


get_pct <- function(code) {
  v <- pooled$`% diff`[pooled$PPTESTCD == code]
  if (length(v) != 1L || is.na(v)) stop("no unique pooled row for ", code)
  v
}

# Gate on the pooled n = 25 column: it is the best-determined reference in
# Table 2. The per-stratum columns rest on n = 7-9 subjects, so their means
# carry roughly +/- 30% sampling noise and gating them would be a coin flip --
# they are shown above for information, not asserted on.
#
# Realised on the authoring run (300 subjects, pooled medians): aucinf.obs
# -0.03%, auclast +0.75%, half.life -1.24%, cl.obs -4.99%, cmax -7.93%. The
# per-stratum medians scatter over roughly -15% to +10%, but that spread is
# dominated by the n = 7-9 published references rather than by the simulation,
# and the pooled median over 300 subjects is far steadier. A mis-transcribed
# clearance, volume, dose or unit moves these by tens of percent, so a 20%
# bound retains plenty of headroom and still goes red on a real error.
stopifnot(
  abs(get_pct("cmax"))       < 20,
  abs(get_pct("auclast"))    < 20,
  abs(get_pct("aucinf.obs")) < 20,
  abs(get_pct("cl.obs"))     < 20,
  abs(get_pct("half.life"))  < 20
)

# The paper states Cmax was reached 2 min after injection (the first sample).
tmax_all <- nca_tbl |> filter(PPTESTCD == "tmax", stratum == "All (n = 25)")
stopifnot(mean(tmax_all$PPORRES == 2) > 0.5)

# Residual-error draws can push a late concentration below zero; LLOQ
# censoring must leave the NCA free of non-finite exposure values. Cmax and
# AUC0-last do not depend on the terminal fit, so every subject must yield a
# finite positive value -- that is the check that catches a censoring bug.
obs_exposure <- nca_tbl |>
  filter(stratum != "All (n = 25)", PPTESTCD %in% c("cmax", "auclast"))
stopifnot(
  nrow(obs_exposure) == 600L,
  !anyNA(obs_exposure$PPORRES),
  all(obs_exposure$PPORRES > 0)
)

# AUC0-inf additionally needs a lambda_z fit, which PKNCA can legitimately
# fail to obtain for the occasional subject whose 120-180 min points are
# non-monotone after residual error. Bound the rate rather than forbidding it
# (realised 1 of 300 subjects, 0.3%, on the authoring run).
aucinf_fail <- nca_tbl |>
  filter(stratum != "All (n = 25)", PPTESTCD == "aucinf.obs")
stopifnot(
  nrow(aucinf_fail) == 300L,
  mean(is.na(aucinf_fail$PPORRES)) < 0.05,
  all(aucinf_fail$PPORRES > 0, na.rm = TRUE)
)

The simulated pooled exposures reproduce Table 2 closely: AUC0-inf agrees to 0.03%, AUC0-last to 0.8%, terminal half-life to 1.2%, Tmax exactly, apparent clearance to 5.0% and Cmax to 7.9%. No row in the table – pooled or per-stratum – exceeds the 20% flag. Per-stratum rows scatter more widely (-15% to +10%), which is what n = 7-9 subjects per published stratum predicts; the model itself makes no distinction between the strata, since weight and age cancel out after per-kilogram dosing.

Two of the remaining gaps have identifiable causes rather than being residual noise. Cmax reads about 8% low because the first sample is at 2 min while the true model maximum occurs at the end of the 30 s infusion, and the published mean is an average over 25 children whereas the table compares medians. The non-compartmental clearance (33.7 mL/min/kg here, 35.5 in the paper) is higher than the model’s structural CL of 31.2 mL/min/kg because a 180 min sampling window truncates the slow third compartment, so the log-linear extrapolation understates AUC0-inf and the derived clearance is biased upward. The paper carries the same discrepancy between its own Table 2 and Table 3.

Assumptions and deviations

BUN normalising median recovered, not printed. Wang 2026 Eq. 5 divides a continuous covariate by its median, but the paper never prints the median BUN; Table 1 reports the cohort mean as 5.2 mmol/L. The value is recovered exactly from the two central volumes the Results text does print: 506 * (3/ref)^(-0.821) = 770 gives ref = 5.003, and 506 * (7/ref)^(-0.821) = 384 gives ref = 5.002. Table 4 independently simulates at BUN 3, 5 and 7 mmol/L with 5 as the central level. The model uses 5 mmol/L, and the Table 4 replication above confirms both published volumes to better than 0.5%.
Two parameters share the value -0.821. Table 3 prints -0.821 for both the BUN-on-V1 exponent and the CL-V1 IIV correlation. This is a genuine coincidence rather than a duplicated cell: the exponent is independently confirmed by the printed volume ratio (0.770 / 0.384 = 2.005 = (3/7)^-0.821), and the correlation is independently stated in the Results text (“a decrease in OFV by 15.039, with a correlation coefficient of -0.821 for the IIV of CL and V1”).
IIV percentages read as log-scale standard deviations. Table 3 reports IIV under a “%” heading and uses the same “%” convention for the residual-error rows, where a percentage can only be a standard deviation. The variances encoded are therefore omega^2 = (percent/100)^2. Reading the percentages instead as the approximate CV of a log-normal (omega = sqrt(log(1 + CV^2))) changes omega by under 2% relative at these magnitudes, so the choice is not load-bearing for any result here.
Additive residual error carries concentration units, not percent. Table 3 prints epsilon_add as 4.83 under a “%” column header carried down from the proportional row above it. An additive residual must have concentration units; 4.83 ug/L sits essentially at the assay LLOQ of 5 ng/mL, the magnitude an additive term takes when set by assay noise at the bottom of the calibration range. It is encoded as 0.00483 mg/L.
Body weight enters structurally, not as a fitted covariate effect. Every Table 3 disposition parameter is per kilogram, so weight multiplies all six with a linear exponent and no reference weight. Concentrations are therefore invariant to weight given per-kilogram dosing. The paper explicitly tested standard allometric scaling and established age-dependent maturation functions on clearance and did not retain either.
Screened-but-unretained covariates are documented, not modelled. Age, sex, BMI, ALT, AST, total bilirubin, albumin, creatinine, haemoglobin and total protein were screened and rejected; they are recorded in the model file’s covariatesDataExcluded list so their provenance survives without appearing in model().
Covariate distributions are assumed. Table 1 reports mean, minimum and maximum per stratum but no distributional shape or correlation structure. Weight, BUN and age are drawn independently from truncated normals centred on the stratum mean with SD = (max - min) / 4. Only BUN affects a prediction.
LLOQ censoring. The paper reports that no observed concentration fell below the 5 ng/mL LLOQ. Simulated residual-error draws at the 180 min tail do occasionally go below zero (about 0.5% of observations on the authoring run), which would make aucinf.obs non-finite. Simulated concentrations are therefore censored at LLOQ/2 = 2.5 ug/L before the NCA, and the gate asserts the resulting exposures are finite and positive.
The 2 min sampling convention is ambiguous. The paper does not state whether “2 min after injection” runs from the start or the end of the 30 s bolus. Measuring from the start gives C2min 2.4-4.8% above Table 4; measuring from the end gives 2.3-4.5% below it. The published values sit between the two. The vignette measures from the start of the infusion, which is the rxode2 convention, and the gate is set at 6% to accommodate the ambiguity.
Terminal half-life agrees at the median but not in spread. PKNCA selects the lambda_z window automatically; the paper regressed “the log-linear phase” without specifying a window. The pooled median half-life matches the published 58.8 min to 1.2%, but the simulated distribution is right-skewed (SD about 36 min against the published 16.8 min, mean about 65 min) because automatic window selection occasionally lands on a short, poorly-conditioned terminal segment. The comparison table and its gate both aggregate by median, which is why this shows as agreement; a mean-based comparison would read about 11% high. For one simulated subject in 300 the lambda_z fit failed outright, which propagates to aucinf.obs, half.life and cl.obs; the gate bounds that failure rate at 5% rather than forbidding it.
No supplement and no erratum. EuropePMC reports hasSuppl = "N" on an open-access record, and no correction notice was found. Every value above comes from the main article’s text, Table 1, Table 2, Table 3 or Table 4.