Piperacillin (Dohmann 2025)

Model and source

mod <- nlmixr2lib::readModelDb("Dohmann_2025_piperacillin")
ui <- rxode2::rxode(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
ui$reference
#> [1] "Dohmann E, Hagel S, Kurlbaum M, Schellong P, Scherf-Clavel O, Surat G. Probability of pharmacokinetic/pharmacodynamic target attainment for different piperacillin/tazobactam dosing regimens in renally impaired patients in a non-intensive care unit setting. Br J Clin Pharmacol. 2025;91(11):3070-3081. doi:10.1002/bcp.70153"

Dohmann 2025 developed a one-compartment population PK model for total serum piperacillin in 49 adults on general (non-intensive-care) wards with chronic kidney disease, including 10 patients on thrice-weekly intermittent haemodialysis. The model was then used for Monte Carlo probability-of-target- attainment (PTA) simulations against a conservative and an aggressive PK/PD target for Pseudomonas aeruginosa.

The structure, as printed in Methods 2.5 (p.3074), is

CL_ind = CL_pop * (eGFR / 21.8 mL/min)^beta_eGFR * exp(eta_CL) + (CL_HD * theta_HD)
Vd_ind = Vd_pop * (BSA  / 1.66 m^2  )^beta_BSA  * exp(eta_Vd)

with theta_HD a 0/1 regressor that switches the haemodialysis clearance on during a dialysis session. The packaged model encodes a BSA reference of 1.88 m^2 rather than the printed 1.66 m^2. That is the single deliberate deviation from the printed source; it is justified in detail, with both readings simulated side by side, in the errata section at the end of this vignette.

ui
#>  ── rxode2-based free-form 1-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>              lcl              lvc lcl_hemodialysis        e_crcl_cl 
#>         1.258461         2.735665         1.376244         0.540000 
#>         e_bsa_vc           propSd 
#>         1.460000         0.260000 
#> 
#> Omega ($omega): 
#>        etalcl etalvc
#> etalcl 0.1156 0.0000
#> etalvc 0.0000 0.0576
#> attr(,"lotriLabels")
#> [1] "Table 2: omega_CL = 0.34 (SE 0.042, RSE 12.4%)"
#> [2] "Table 2: omega_Vd = 0.24 (SE 0.056, RSE 23.2%)"
#> attr(,"lotriFix")
#>        etalcl etalvc
#> etalcl  FALSE  FALSE
#> etalvc  FALSE  FALSE
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1          central
#>  ── μ-referencing ($muRefTable): ──  
#>   theta    eta level
#> 1   lcl etalcl    id
#> 2   lvc etalvc    id
#> 
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     compartmentData <- list(central = list(analyte = "piperacillin", 
#>         units = "mg", specimen = "serum", verified = TRUE))
#>     covariateData <- list(BSA = list(description = "Body surface area (Mosteller formula)", 
#>         units = "m^2", type = "continuous", reference_category = NULL, 
#>         notes = "Dohmann 2025 Methods 2.5: BSA computed with the Mosteller formula and retained as the only covariate on Vd, entering as the power term (BSA / 1.88 m^2)^beta_BSA. IMPORTANT -- the reference value encoded here (1.88 m^2) is NOT the value printed in the paper's typeset Vd equation on p.3074, which reads 1.66 m^2. The printed 1.66 m^2 is treated as a typographical error and 1.88 m^2 is encoded instead, per operator decision (sidecar request-001 / response-001, question q1, answer A). Five independent lines of evidence inside the paper itself support 1.88 m^2: (i) Methods 2.6 states the simulation cohort drew BSA from a lognormal distribution with MEAN 1.88 m^2 and SD 0.12 m^2, i.e. the authors' own cohort-typical BSA; (ii) Table 1 mean height 170.9 cm and mean weight 76.1 kg give a Mosteller BSA of sqrt(170.9 x 76.1 / 3600) = 1.90 m^2, so a cohort median of 1.66 m^2 is not attainable; (iii) the Discussion states 'the population estimate for volume of distribution of the herein presented model is 15.4 L for an adult patient', which only holds if the 15.42 L typical value corresponds to a cohort-typical rather than a small adult; (iv) the companion eGFR reference in the SAME pair of equations (21.8 mL/min) is demonstrably the cohort central value on the individualized scale (Table 1 mean 19.9 mL/min/1.73 m^2 x 1.90 / 1.73 = 21.85 mL/min), so the same convention applied to BSA gives 1.88-1.90 m^2; (v) reimplementing the paper's own Monte Carlo simulation reproduces the six mid-range cells of Table 3 to within about 1 percentage point with 1.88 m^2 but over-predicts every cell with 1.66 m^2. See the vignette 'Assumptions, deviations and errata' section, which reports both readings side by side.", 
#>         source_name = "BSA"), CRCL = list(description = "Estimated glomerular filtration rate, MDRD formula individualized using body surface area (i.e. de-normalized to absolute mL/min, NOT mL/min/1.73 m^2)", 
#>         units = "mL/min", type = "continuous", reference_category = NULL, 
#>         notes = "Dohmann 2025 Methods 2.5 lists the available covariate as 'eGFR (Modification of Diet in Renal Disease formula individualized using BSA)', and the typeset CL equation on p.3074 divides by '21.8 mL/min' -- absolute mL/min units, not mL/min/1.73 m^2. This is the raw un-normalized variant of the CRCL canonical (same convention as Georges_2009_ceftazidime.R and Delattre_2010_amikacin.R); supplying a BSA-normalized value here would silently rescale the renal term. Consistency check on the reference value: Table 1 reports a cohort mean eGFR of 19.9 mL/min/1.73 m^2, which individualized to the cohort-typical BSA gives 19.9 x 1.90 / 1.73 = 21.85 mL/min, matching the printed 21.8 mL/min. Study eGFR was computed with MDRD at Wuerzburg and CKD-EPI at Jena (Table 1 note). Enters as the power term (CRCL / 21.8)^beta_eGFR on the non-dialysis clearance only; the haemodialysis arm is renal-function independent. In the paper's own Monte Carlo simulations (Methods 2.6) eGFR was 'set exactly to the value according to the group', i.e. the group labels 40 / 30 / 20 / 10 mL/min are covariate values on this absolute scale.", 
#>         source_name = "eGFR"), RRT_HEMODIAL_ACTIVE = list(description = "Haemodialysis-active indicator (1 while an intermittent haemodialysis session is running, 0 otherwise)", 
#>         units = "(binary)", type = "binary", reference_category = "0 (interdialytic / no dialysis running)", 
#>         notes = "Time-varying within subject. Dohmann 2025 Methods 2.5: 'Patients with iHD were modelled by including a second clearance process during the times of dialysis. The dialysis clearance was turned on or off by a regressor variable (theta_HD = 0 or 1) according to the documented dialysis times.' The paper's own symbol is theta_HD; it is a data-supplied 0/1 regressor, not an estimated parameter, so it is registered here as a covariate column. Group 3 patients received thrice-weekly intermittent haemodialysis with a mean session length of 258.8 +/- 24.8 min (Table 1); the Monte Carlo simulations (Methods 2.6) used a 4-h session once daily beginning directly after the end of the piperacillin infusion. Gates the additive cl_hemodialysis arm exactly as in Veinstein_2013_gentamicin.R and VanWart_2025_telavancin.R.", 
#>         source_name = "theta_HD"))
#>     covariatesDataExcluded <- list(AGE = list(description = "Age", 
#>         units = "years", type = "continuous", notes = "Dohmann 2025 Methods 2.5 lists age among the available covariates for model building; Results 3.2 reports that 'in the stepwise testing of covariates, only BSA and eGFR were significant covariates', so age was screened and not retained. No point estimate is reported (see Table S1 / Figure S4, not on disk)."), 
#>         BMI = list(description = "Body mass index", units = "kg/m^2", 
#>             type = "continuous", notes = "Dohmann 2025 Methods 2.5: BMI was among the covariates considered on Vd; screened and not retained in the final model (Results 3.2)."), 
#>         HT = list(description = "Height", units = "cm", type = "continuous", 
#>             notes = "Dohmann 2025 Methods 2.5: height was among the covariates considered on Vd; screened and not retained in the final model (Results 3.2). Height is nevertheless an input to the Mosteller BSA that IS retained."), 
#>         WT = list(description = "Body weight", units = "kg", 
#>             type = "continuous", notes = "Dohmann 2025 Methods 2.5: weight was among the covariates considered on Vd; screened and not retained in the final model (Results 3.2). Weight is nevertheless an input to the Mosteller BSA that IS retained."), 
#>         SEXF = list(description = "Female sex indicator (1 = female, 0 = male)", 
#>             units = "(binary)", type = "binary", notes = "Dohmann 2025 Methods 2.5 gives the categorical-covariate form explicitly -- Vd_ind = Vd_pop * exp(beta_sex * theta_sex) * exp(eta_Vd) with theta_sex = 1 for female and 0 for male -- and tested sex on both CL and Vd. Results 3.2: only BSA and eGFR were significant, so no sex effect is present in the final model and no beta_sex estimate is reported."))
#>     description <- "One-compartment population PK model for intravenous piperacillin (administered as piperacillin/tazobactam) in adult non-intensive-care ward patients with chronic kidney disease, including patients on thrice-weekly intermittent haemodialysis (n=49; 135 total-piperacillin serum samples; Dohmann 2025). Elimination is linear and is decomposed into a non-dialysis clearance carrying a power effect of BSA-individualized MDRD eGFR, plus an additive haemodialysis-arm clearance gated on/off by the time-varying RRT_HEMODIAL_ACTIVE regressor. Central volume carries a power effect of body surface area. Residual error is proportional; there is no interindividual variability on the haemodialysis clearance. The model was built in Monolix 2023R1 and was used for Monte Carlo probability-of-target-attainment simulations against a conservative (60% fT > MIC) and an aggressive (100% fT > 4 x MIC) PK/PD target for Pseudomonas aeruginosa."
#>     population <- list(species = "human", n_subjects = 49L, n_studies = 1L, 
#>         age_range = "35-88 years (mean 69.3 +/- 11.1; Table 1 and Discussion, strengths (i))", 
#>         weight_range = "mean 76.1 +/- 17.1 kg; BMI mean 26.1 +/- 5.6 kg/m^2 spanning <18.5 to >40 kg/m^2 (Table 1 and Discussion, strengths (i))", 
#>         sex_female_pct = 36.7, race_ethnicity = "Not reported (two German university hospitals)", 
#>         disease_state = "Adults on general (non-intensive-care) wards with chronic kidney disease receiving intermittent intravenous piperacillin/tazobactam. Three prespecified renal-function groups: group 1 eGFR 20-40 mL/min/1.73 m^2 (n=20), group 2 eGFR <20 mL/min/1.73 m^2 (n=19), group 3 thrice-weekly intermittent haemodialysis (n=10). Indications included empirical therapy for nosocomial infection of unclear focus, and targeted therapy for hospital-acquired pneumonia or severe skin infection with or without bacteraemia (Methods 2.3). Exclusion criteria: pregnancy or breastfeeding, age under 18 years, beta-lactam hypersensitivity, participation in another clinical trial.", 
#>         dose_range = "Piperacillin/tazobactam 4.5 g (4000 mg piperacillin / 500 mg tazobactam) as a 30-min intermittent intravenous infusion. Group 1 q8h, groups 2 and 3 q12h, per the Summary of Product Characteristics; group 3 additionally received 2.25 g after haemodialysis (Table 1). A minority of patients deviated from the recommended interval (six patients in group 1 received q12h, three in group 2 received q8h, one in group 3 received q24h). Only piperacillin was quantified.", 
#>         regions = "Germany (Wuerzburg University Hospital n=30, Jena University Hospital n=19)", 
#>         renal_function = "Chronic kidney disease with eGFR <40 mL/min/1.73 m^2, or thrice-weekly intermittent haemodialysis. Cohort mean eGFR 19.9 +/- 9.3 mL/min/1.73 m^2 (group 1 29.2 +/- 6.6, group 2 14.6 +/- 3.6, group 3 11.2 +/- 3.3). eGFR was calculated with MDRD at Wuerzburg and CKD-EPI at Jena; creatinine clearance was not measured because it is not routinely determined in non-ICU patients (Discussion, limitations (i)).", 
#>         notes = "Prospective multicentre study conducted 2 March 2022 to 17 March 2023 (ethics approvals 110/21 and 2021-2399). Sampling at steady state after 3 days of therapy at three time points for groups 1 and 2 (within 30 min pre-dose, within 30 min after end of infusion, and 4 or 6 h post-infusion) and immediately before and after a dialysis session for group 3. A total of 135 samples were analysed (80 from 30 Wuerzburg patients and 55 from 19 Jena patients); one sample at each site was excluded for an implausibly high concentration. Bioanalysis by HPLC-MS/MS (Sciex QTRAP 4500MD with Agilent 1290 UHPLC), analytical range 0.5-190.0 mg/L, validated to EMA guidance. Concentrations are reported and modelled as TOTAL piperacillin; the paper converts to unbound by multiplying by 0.8 (an 80% unbound fraction taken from the literature), so the free target of 16 mg/L corresponds to a total concentration of 19.8-20 mg/L (Methods 2.4 and 2.6). Model estimated in Monolix 2023R1; simulations in Simulx 2023R1. Supplementary Table S1 and Figures S1-S4 (covariate screening and goodness-of-fit diagnostics) are not on disk; they contain no final parameter values, all of which are in Table 2.")
#>     reference <- "Dohmann E, Hagel S, Kurlbaum M, Schellong P, Scherf-Clavel O, Surat G. Probability of pharmacokinetic/pharmacodynamic target attainment for different piperacillin/tazobactam dosing regimens in renally impaired patients in a non-intensive care unit setting. Br J Clin Pharmacol. 2025;91(11):3070-3081. doi:10.1002/bcp.70153"
#>     units <- list(time = "h", dosing = "mg", concentration = "mg/L")
#>     vignette <- "Dohmann_2025_piperacillin"
#>     ini({
#>         lcl <- 1.25846098961001
#>         label("Non-dialysis clearance CL_pop at eGFR 21.8 mL/min (L/h)")
#>         lvc <- 2.73566536813518
#>         label("Central volume of distribution Vd_pop at BSA 1.88 m^2 (L)")
#>         lcl_hemodialysis <- 1.37624402526639
#>         label("Additive haemodialysis-arm clearance CL_HD (L/h)")
#>         e_crcl_cl <- 0.54
#>         label("Power exponent of (CRCL / 21.8 mL/min) on non-dialysis CL (unitless)")
#>         e_bsa_vc <- 1.46
#>         label("Power exponent of (BSA / 1.88 m^2) on Vd (unitless)")
#>         propSd <- c(0, 0.26)
#>         label("Proportional residual error (fraction)")
#>         etalcl ~ 0.1156
#>         label("Table 2: omega_CL = 0.34 (SE 0.042, RSE 12.4%)")
#>         etalvc ~ 0.0576
#>         label("Table 2: omega_Vd = 0.24 (SE 0.056, RSE 23.2%)")
#>     })
#>     model({
#>         cl <- exp(lcl + etalcl) * (CRCL/21.8)^e_crcl_cl
#>         vc <- exp(lvc + etalvc) * (BSA/1.88)^e_bsa_vc
#>         cl_hemodialysis <- exp(lcl_hemodialysis)
#>         cl_total <- cl + RRT_HEMODIAL_ACTIVE * cl_hemodialysis
#>         kel <- cl_total/vc
#>         d/dt(central) <- -kel * central
#>         Cc <- central/vc
#>         Cc ~ prop(propSd)
#>     })
#> }

Source trace

Every value in ini() and every equation in model(), with its location in Dohmann 2025.

Source trace for Dohmann_2025_piperacillin.

Quantity	Value	Source location
1-compartment, linear elimination	structure	Methods 2.5; Results 3.2 (‘Based on the available data, a 1-compartment model was chosen’)
Second, additive clearance during dialysis, gated by a 0/1 regressor	structure	Methods 2.5 (‘Patients with iHD were modelled by including a second clearance process during the times of dialysis’)
CL covariate form: power on eGFR, reference 21.8 mL/min	equation	Methods 2.5, first typeset display equation, p.3074
Vd covariate form: power on BSA	equation	Methods 2.5, second typeset display equation, p.3074
BSA reference divisor	1.88 m^2 (paper prints 1.66 m^2)	Methods 2.6 simulation-cohort mean; see errata section
lcl = log(3.52)	CL_pop = 3.52 L/h (RSE 5.68%)	Table 2
lvc = log(15.42)	Vd_pop = 15.42 L (RSE 7.02%)	Table 2
lcl_hemodialysis = log(3.96)	CL_HD = 3.96 L/h (RSE 34.4%)	Table 2
e_crcl_cl = 0.54	beta_eGFR = 0.54 (RSE 20.0%)	Table 2
e_bsa_vc = 1.46	beta_BSA = 1.46 (RSE 28.7%)	Table 2
etalcl ~ 0.1156	omega_CL = 0.34 SD -> variance 0.34^2	Table 2 and its note (‘standard deviation of random effects’)
etalvc ~ 0.0576	omega_Vd = 0.24 SD -> variance 0.24^2	Table 2 and its note
No IIV on CL_HD	omitted	Results 3.2 (‘Interindividual variability on haemodialysis clearance was not included in the model’)
propSd = 0.26	proportional error 0.26 (RSE 9.38%)	Table 2; Results 3.2 (proportional error model selected)
Unbound fraction 0.80 (used for targets, not in the model)	0.80	Methods 2.6
Conservative target: free >16 mg/L for >=60% of tau	target	Methods 2.2 and 2.6
Aggressive target: free >4 x MIC (64 mg/L) for 100% of tau	target	Methods 2.2 and 2.6

Note that omega_CL and omega_Vd are standard deviations, not variances: the Table 2 note says so explicitly, and the reported relative standard errors corroborate it. A variance estimated from n = 49 subjects cannot have a relative standard error much below sqrt(2/49) = 20.2%, yet omega_CL is reported with an RSE of 12.4%, which is only attainable on the SD scale.

Population

Study population (Dohmann 2025 Methods 2.1 and Table 1).

Field	Value
species	human
n_subjects	49
n_studies	1
age_range	35-88 years (mean 69.3 +/- 11.1; Table 1 and Discussion, strengths (i))
weight_range	mean 76.1 +/- 17.1 kg; BMI mean 26.1 +/- 5.6 kg/m^2 spanning <18.5 to >40 kg/m^2 (Table 1 and Discussion, strengths (i))
sex_female_pct	36.7
race_ethnicity	Not reported (two German university hospitals)
disease_state	Adults on general (non-intensive-care) wards with chronic kidney disease receiving intermittent intravenous piperacillin/tazobactam. Three prespecified renal-function groups: group 1 eGFR 20-40 mL/min/1.73 m^2 (n=20), group 2 eGFR <20 mL/min/1.73 m^2 (n=19), group 3 thrice-weekly intermittent haemodialysis (n=10). Indications included empirical therapy for nosocomial infection of unclear focus, and targeted therapy for hospital-acquired pneumonia or severe skin infection with or without bacteraemia (Methods 2.3). Exclusion criteria: pregnancy or breastfeeding, age under 18 years, beta-lactam hypersensitivity, participation in another clinical trial.
dose_range	Piperacillin/tazobactam 4.5 g (4000 mg piperacillin / 500 mg tazobactam) as a 30-min intermittent intravenous infusion. Group 1 q8h, groups 2 and 3 q12h, per the Summary of Product Characteristics; group 3 additionally received 2.25 g after haemodialysis (Table 1). A minority of patients deviated from the recommended interval (six patients in group 1 received q12h, three in group 2 received q8h, one in group 3 received q24h). Only piperacillin was quantified.
regions	Germany (Wuerzburg University Hospital n=30, Jena University Hospital n=19)
renal_function	Chronic kidney disease with eGFR <40 mL/min/1.73 m^2, or thrice-weekly intermittent haemodialysis. Cohort mean eGFR 19.9 +/- 9.3 mL/min/1.73 m^2 (group 1 29.2 +/- 6.6, group 2 14.6 +/- 3.6, group 3 11.2 +/- 3.3). eGFR was calculated with MDRD at Wuerzburg and CKD-EPI at Jena; creatinine clearance was not measured because it is not routinely determined in non-ICU patients (Discussion, limitations (i)).
notes	Prospective multicentre study conducted 2 March 2022 to 17 March 2023 (ethics approvals 110/21 and 2021-2399). Sampling at steady state after 3 days of therapy at three time points for groups 1 and 2 (within 30 min pre-dose, within 30 min after end of infusion, and 4 or 6 h post-infusion) and immediately before and after a dialysis session for group 3. A total of 135 samples were analysed (80 from 30 Wuerzburg patients and 55 from 19 Jena patients); one sample at each site was excluded for an implausibly high concentration. Bioanalysis by HPLC-MS/MS (Sciex QTRAP 4500MD with Agilent 1290 UHPLC), analytical range 0.5-190.0 mg/L, validated to EMA guidance. Concentrations are reported and modelled as TOTAL piperacillin; the paper converts to unbound by multiplying by 0.8 (an 80% unbound fraction taken from the literature), so the free target of 16 mg/L corresponds to a total concentration of 19.8-20 mg/L (Methods 2.4 and 2.6). Model estimated in Monolix 2023R1; simulations in Simulx 2023R1. Supplementary Table S1 and Figures S1-S4 (covariate screening and goodness-of-fit diagnostics) are not on disk; they contain no final parameter values, all of which are in Table 2.

Table 1 of the paper reports mean height 170.9 cm and mean weight 76.1 kg, which under the Mosteller formula the Methods name gives a cohort-typical BSA of 1.90 m^2. Cohort mean eGFR is 19.9 mL/min/1.73 m^2; the model covariate is the BSA-individualized (de-normalized) eGFR in absolute mL/min, so the cohort-typical covariate value is 21.9 mL/min, matching the 21.8 mL/min divisor printed in the CL equation.

Virtual cohort and simulation machinery

The paper’s Monte Carlo design (Methods 2.6) drew BSA from a lognormal distribution with mean 1.88 m^2 and SD 0.12 m^2, set eGFR exactly to the group value, and propagated the modelled IIV together with the parameter-estimation uncertainty from the Table 2 standard errors. Residual error is not part of a PTA calculation, so the simulations below use the individual predicted concentration Cc with the random effects and residual error zeroed via rxode2::zeroRe(), and supply per-subject parameters explicitly.

uiz <- rxode2::zeroRe(ui)

# Table 2 point estimates and standard errors.
est <- c(vd = 15.42, bbsa = 1.46, cl = 3.52, begfr = 0.54, clhd = 3.96)
se  <- c(vd = 1.08,  bbsa = 0.42, cl = 0.20, begfr = 0.11, clhd = 1.36)

# 200 virtual patients per arm - the vignette cap. The paper used 10 trials of
# 100 patients each; Tables 3 and 4 report the pooled mean, which 200 patients
# estimate without bias but with more Monte Carlo noise (about 1-3 percentage
# points for a PTA near 90%).
n_sub <- 200L

# Each cohort is seeded from its own eGFR group rather than inheriting a single
# global seed. With one global seed every draw below would depend on the
# cumulative RNG consumption of every chunk above it, so merely inserting or
# reordering a chunk would silently change every downstream assertion. Seeding
# per group makes each cohort reproducible independently of what precedes it,
# and gives common random numbers to the two tables that share a group.
make_cohort <- function(n, egfr_group) {
  set.seed(20250823L + as.integer(egfr_group))
  # lognormal BSA with arithmetic mean 1.88 and SD 0.12 (Methods 2.6)
  bsa_mu <- log(1.88^2 / sqrt(0.12^2 + 1.88^2))
  bsa_sd <- sqrt(log(1 + 0.12^2 / 1.88^2))
  bsa <- exp(stats::rnorm(n, bsa_mu, bsa_sd))
  data.frame(
    id   = seq_len(n),
    # The renal-function GROUPS are defined on the body-surface-normalized
    # scale - Table 1 reports eGFR in mL/min/1.73 m^2 and the groups are the
    # bands "20-40" and "<20" on that scale - whereas the model COVARIATE is
    # the BSA-individualized (absolute) eGFR in mL/min, which is why the CL
    # equation's reference divisor is printed as 21.8 mL/min rather than
    # 19.9 mL/min/1.73 m^2. So a group label is converted to the covariate by
    # multiplying by BSA / 1.73, exactly the conversion that reconciles the
    # cohort mean 19.9 with the printed 21.8 reference. Setting CRCL to the
    # bare group label instead leaves the whole cohort about 8.7% low on
    # clearance and over-predicts both published PTA tables.
    CRCL = egfr_group * bsa / 1.73,
    BSA  = bsa,
    # Structural parameters carry BOTH the estimation uncertainty (Table 2 SE)
    # and the interindividual variability (Table 2 omega, on the log scale).
    lcl = log(pmax(0.1, est[["cl"]] + se[["cl"]] * stats::rnorm(n)) *
                exp(stats::rnorm(n, 0, 0.34))),
    lvc = log(pmax(1.0, est[["vd"]] + se[["vd"]] * stats::rnorm(n)) *
                exp(stats::rnorm(n, 0, 0.24))),
    e_crcl_cl        = est[["begfr"]] + se[["begfr"]] * stats::rnorm(n),
    e_bsa_vc         = est[["bbsa"]]  + se[["bbsa"]]  * stats::rnorm(n),
    lcl_hemodialysis = log(pmax(0.1, est[["clhd"]] + se[["clhd"]] * stats::rnorm(n)))
  )
}

par_cols <- c("id", "lcl", "lvc", "e_crcl_cl", "e_bsa_vc", "lcl_hemodialysis")

One cohort is drawn per eGFR group and reused across every regimen tested in that group (common random numbers), so regimen-to-regimen differences are not contaminated by resampling noise.

# Steady-state solve for a non-dialysis arm. rxode2's ss = 1 gives the exact
# analytic steady state, so no run-in period is needed.
solve_ss <- function(co, dose, tinf, tau, grid = 0.1, bsa_scale = 1) {
  ev <- rxode2::et(amt = dose, dur = tinf, ii = tau, ss = 1, cmt = "central") |>
    rxode2::et(seq(0, tau, by = grid), cmt = "central") |>
    rxode2::et(id = co$id)
  evdf <- as.data.frame(ev)
  evdf <- dplyr::left_join(evdf, co[, c("id", "CRCL", "BSA")], by = "id")
  # Emulating an alternative BSA reference r is exact via a rescale, because
  # (BSA / r)^b == ((BSA * 1.88 / r) / 1.88)^b.
  evdf$BSA <- evdf$BSA * bsa_scale
  evdf$RRT_HEMODIAL_ACTIVE <- 0
  suppressWarnings(
    rxode2::rxSolve(uiz, params = co[, par_cols], events = evdf,
                    returnType = "data.frame")
  )
}

# Fraction of the dosing interval spent above a TOTAL-concentration threshold,
# and the resulting PTA. Methods 2.6 evaluates on a 0.1 h grid and integrates a
# 0/1 indicator over the interval.
pta_from <- function(sol, thresh_total, frac_required) {
  per_id <- sol |>
    dplyr::group_by(id) |>
    dplyr::summarise(frac = mean(Cc > thresh_total), .groups = "drop")
  100 * mean(per_id$frac >= frac_required)
}

# Methods 2.4 and 2.6: the model predicts TOTAL concentrations and the paper
# converts with an unbound fraction of 0.80 taken from the literature. The
# paper's own stated equivalents are used here rather than a recomputed 16/0.8:
# Methods 2.2 and the Figure 3 caption give "19.8 mg/L equivalent to 16 mg/L
# unbound", and the Figure 5 caption gives 79.2 mg/L for 4 x MIC. (A literal
# 16/0.8 would give 20.0 and 80.0; the paper's numbers imply fu = 0.808. The
# 1% difference is immaterial to PTA, but quoting the paper's own thresholds
# keeps this vignette auditable against its figures.)
fu_total <- 19.8 / 16                # 1.2375; the paper's total/unbound factor
thresh_conservative <- 19.8          # Methods 2.2, Figure 3 caption
thresh_aggressive   <- 79.2          # Figure 5 caption (4 x MIC, total)

Structural check: solver against the closed form

At steady state a one-compartment model with zero-order input has an exact closed form. Both sides of this check use the same drawn per-subject parameters, so the only difference is numerical solver error and a tight bound is appropriate.

co30 <- make_cohort(n_sub, 30)
sol30 <- solve_ss(co30, dose = 4000, tinf = 0.5, tau = 8)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

cl_i <- exp(co30$lcl) * (co30$CRCL / 21.8)^co30$e_crcl_cl
vc_i <- exp(co30$lvc) * (co30$BSA / 1.88)^co30$e_bsa_vc
k_i  <- cl_i / vc_i
rate <- 4000 / 0.5
cmax_closed <- rate / cl_i * (1 - exp(-k_i * 0.5)) / (1 - exp(-k_i * 8))
cmin_closed <- cmax_closed * exp(-k_i * 7.5)

cmp_struct <- sol30 |>
  dplyr::group_by(id) |>
  dplyr::summarise(cmax = max(Cc), cmin = min(Cc), .groups = "drop") |>
  dplyr::mutate(cmax_err = cmax / cmax_closed - 1,
                cmin_err = cmin / cmin_closed - 1)

sprintf("max |rel. error|: Cmax,ss %.2e, Cmin,ss %.2e",
        max(abs(cmp_struct$cmax_err)), max(abs(cmp_struct$cmin_err)))
#> [1] "max |rel. error|: Cmax,ss 2.66e-15, Cmin,ss 4.77e-15"
stopifnot(
  max(abs(cmp_struct$cmax_err)) < 1e-8,
  max(abs(cmp_struct$cmin_err)) < 1e-8
)

PKNCA validation

Steady-state NCA over one dosing interval, for the three non-dialysis eGFR groups on the 4.5 g q8h 30-min regimen. At steady state AUC(0-tau) must equal Dose / CL exactly, which makes the NCA an end-to-end audit of the clearance parameterisation including the eGFR covariate term.

nca_groups <- c(40, 30, 20)
cohorts <- lapply(nca_groups, function(g) make_cohort(n_sub, g))
names(cohorts) <- as.character(nca_groups)

nca_conc <- dplyr::bind_rows(lapply(nca_groups, function(g) {
  # A finer grid over the infusion phase: the linear-up trapezoid is exact on
  # the log-linear decay but not on the concave zero-order input ramp.
  co <- cohorts[[as.character(g)]]
  ev <- rxode2::et(amt = 4000, dur = 0.5, ii = 8, ss = 1, cmt = "central") |>
    rxode2::et(sort(unique(c(seq(0, 0.5, by = 0.01), seq(0, 8, by = 0.05)))),
               cmt = "central") |>
    rxode2::et(id = co$id)
  evdf <- as.data.frame(ev)
  evdf <- dplyr::left_join(evdf, co[, c("id", "CRCL", "BSA")], by = "id")
  evdf$RRT_HEMODIAL_ACTIVE <- 0
  suppressWarnings(
    rxode2::rxSolve(uiz, params = co[, par_cols], events = evdf,
                    returnType = "data.frame")
  ) |>
    dplyr::mutate(egfr_group = paste0("eGFR ", g, " mL/min"))
}))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

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

dose_df <- sim_nca |>
  dplyr::distinct(id, egfr_group) |>
  dplyr::mutate(time = 0, amt = 4000)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | egfr_group + id,
                             concu = "mg/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | egfr_group + id,
                             doseu = "mg")

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

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

auc_i <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD == "auclast") |>
  dplyr::select(egfr_group, id, auclast = PPORRES)

cl_all <- dplyr::bind_rows(lapply(nca_groups, function(g) {
  co <- cohorts[[as.character(g)]]
  data.frame(
    egfr_group = paste0("eGFR ", g, " mL/min"),
    id = co$id,
    cl = exp(co$lcl) * (co$CRCL / 21.8)^co$e_crcl_cl
  )
}))

ident <- auc_i |>
  dplyr::inner_join(cl_all, by = c("egfr_group", "id")) |>
  dplyr::mutate(auc_closed = 4000 / cl,
                rel_err = auclast / auc_closed - 1)

sprintf("AUC(0-tau) vs Dose/CL: max |rel. error| = %.2e across %d subjects",
        max(abs(ident$rel_err)), nrow(ident))
#> [1] "AUC(0-tau) vs Dose/CL: max |rel. error| = 6.21e-06 across 600 subjects"
stopifnot(max(abs(ident$rel_err)) < 1e-3)

The residual is pure trapezoidal discretisation on the zero-order input ramp; the identity itself holds exactly.

Comparison against reference NCA

Dohmann 2025 reports no NCA table, so the reference column below is the closed form implied by the paper’s own parameterisation - AUC(0-tau) = Dose / CL_ind and Cav = AUC(0-tau) / tau - evaluated on the same virtual cohort.

reference <- ident |>
  dplyr::group_by(egfr_group) |>
  dplyr::summarise(auclast = stats::median(auc_closed),
                   cav = stats::median(auc_closed) / 8,
                   .groups = "drop") |>
  as.data.frame()

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = reference,
  by = "egfr_group",
  params = c("auclast", "cav", "cmax", "cmin", "tmax"),
  units = c(auclast = "mg*h/L", cav = "mg/L", cmax = "mg/L",
            cmin = "mg/L", tmax = "h"),
  tolerance_pct = 20,
  label_first_column = "NCA parameter"
)

knitr::kable(
  cmp,
  caption = paste(
    "Simulated steady-state NCA (4.5 g q8h over 30 min) against the closed form",
    "implied by Dohmann 2025 Table 2 on the same cohort.",
    "* differs from reference by >20%. Cmax, Cmin and Tmax have no closed-form",
    "reference column and are reported for orientation only."
  ),
  digits = 3
)

Simulated steady-state NCA (4.5 g q8h over 30 min) against the closed form implied by Dohmann 2025 Table 2 on the same cohort. * differs from reference by >20%. Cmax, Cmin and Tmax have no closed-form reference column and are reported for orientation only.

NCA parameter	egfr_group	Reference	Simulated	% diff
AUClast (mg*h/L)	eGFR 20 mL/min	1180	1180	-0.0%
AUClast (mg*h/L)	eGFR 30 mL/min	937	937	-0.0%
AUClast (mg*h/L)	eGFR 40 mL/min	776	776	-0.0%
Cavg (mg/L)	eGFR 20 mL/min	148	148	-0.0%
Cavg (mg/L)	eGFR 30 mL/min	117	117	-0.0%
Cavg (mg/L)	eGFR 40 mL/min	97	97	-0.0%

attr(cmp, "footnote")
#> NULL

Orientation against the observed concentrations

Results 3.1 reports observed mean concentrations in the study cohort. These are not a like-for-like NCA comparison - the observed cohort had heterogeneous dosing (six group-1 patients received q12h and three group-2 patients received q8h), the “peak” was drawn within 30 min after the end of infusion rather than at the exact Cmax, and the mid-interval sample was taken at either 4 or 6 h. The comparison is nevertheless a useful check that the model is on the right absolute scale.

# Group-typical covariate values from Table 1, with eGFR de-normalized to the
# absolute mL/min scale the model covariate uses (x cohort BSA / 1.73).
obs_cmp <- dplyr::bind_rows(
  lapply(list(
    list(grp = "Group 1 (eGFR 20-40)", egfr = 29.2 * 1.90 / 1.73, tau = 8,
         peak_obs = 287.1, mid_obs = 74.6),
    list(grp = "Group 2 (eGFR <20)",   egfr = 14.6 * 1.90 / 1.73, tau = 12,
         peak_obs = 282.6, mid_obs = 97.0)
  ), function(g) {
    cl <- 3.52 * (g$egfr / 21.8)^0.54
    k  <- cl / 15.42
    r  <- 4000 / 0.5
    cmax <- r / cl * (1 - exp(-k * 0.5)) / (1 - exp(-k * g$tau))
    data.frame(
      Group = g$grp,
      `Peak observed (mg/L)` = g$peak_obs,
      `Peak typical-value predicted (mg/L)` = round(cmax, 1),
      `Mid-interval observed (mg/L)` = g$mid_obs,
      `Mid-interval predicted, 4-6 h post-infusion (mg/L)` =
        paste0(round(cmax * exp(-k * 6), 1), " to ", round(cmax * exp(-k * 4), 1)),
      check.names = FALSE
    )
  })
)
knitr::kable(obs_cmp, caption = paste(
  "Typical-value steady-state predictions against the observed group means of",
  "Dohmann 2025 Results 3.1."
))

Typical-value steady-state predictions against the observed group means of Dohmann 2025 Results 3.1.

Group	Peak observed (mg/L)	Peak typical-value predicted (mg/L)	Mid-interval observed (mg/L)	Mid-interval predicted, 4-6 h post-infusion (mg/L)
Group 1 (eGFR 20-40)	287.1	270.5	74.6	50.1 to 87.9
Group 2 (eGFR <20)	282.6	274.2	97.0	85.9 to 126.5

peak_err <- abs(as.numeric(obs_cmp[[3]]) / as.numeric(obs_cmp[[2]]) - 1)
sprintf("peak relative error: %s",
        paste(sprintf("%.1f%%", 100 * peak_err), collapse = ", "))
#> [1] "peak relative error: 5.8%, 3.0%"

# This check is fully deterministic - both sides are typical-value quantities
# with no random draw on either side (the published observed means are fixed
# numbers, and the prediction uses the Table 2 point estimates at the group
# covariate values). A tight bound is therefore appropriate here, unlike the
# cohort-level PTA gates below, and cannot be broken by reseeding.
stopifnot(all(peak_err < 0.08))

Both observed group-mean peaks are reproduced to within 6%. Observed troughs are not used as a check: in both groups the reported standard deviation is of the same size as the mean (27.8 +/- 35.9 and 58.2 +/- 57.9 mg/L), so the trough distribution is far too skewed to discriminate between candidate models.

Replicating Table 3 (conservative PK/PD target)

Table 3 is the paper’s headline answer key: PTA for free piperacillin above 16 mg/L for at least 60% of the dosing interval. This is the single strongest external validation available for the model, because it depends jointly on the clearance, the volume, both covariate terms, and both variance components.

regimens3 <- list(
  list(lab = "4.5 g q8 h over 30 min",  dose = 4000, tinf = 0.5, tau = 8),
  list(lab = "4.5 g q8 h over 3 h",     dose = 4000, tinf = 3.0, tau = 8),
  list(lab = "4.5 g q12 h over 30 min", dose = 4000, tinf = 0.5, tau = 12),
  list(lab = "4.5 g q12 h over 3 h",    dose = 4000, tinf = 3.0, tau = 12)
)
egfr3 <- c(40, 30, 20)

pta_grid <- function(bsa_scale) {
  dplyr::bind_rows(lapply(egfr3, function(g) {
    co <- cohorts[[as.character(g)]]
    dplyr::bind_rows(lapply(regimens3, function(r) {
      sol <- solve_ss(co, r$dose, r$tinf, r$tau, bsa_scale = bsa_scale)
      data.frame(egfr = g, regimen = r$lab,
                 pta = pta_from(sol, thresh_conservative, 0.60))
    }))
  }))
}

sim_188 <- pta_grid(1)             # reference 1.88 m^2 (as packaged)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
sim_166 <- pta_grid(1.88 / 1.66)   # reference 1.66 m^2 (as printed)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

published3 <- tibble::tribble(
  ~egfr, ~regimen,                     ~published,
  40, "4.5 g q8 h over 30 min",  91.6,
  40, "4.5 g q8 h over 3 h",     98.7,
  40, "4.5 g q12 h over 30 min", 61.2,
  40, "4.5 g q12 h over 3 h",    79.2,
  30, "4.5 g q8 h over 30 min",  95.8,
  30, "4.5 g q8 h over 3 h",     99.4,
  30, "4.5 g q12 h over 30 min", 73.8,
  30, "4.5 g q12 h over 3 h",    87.7,
  20, "4.5 g q8 h over 30 min",  99.2,
  20, "4.5 g q8 h over 3 h",     99.9,
  20, "4.5 g q12 h over 30 min", 90.0,
  20, "4.5 g q12 h over 3 h",    96.9
)

tab3 <- published3 |>
  dplyr::left_join(dplyr::rename(sim_188, sim188 = pta), by = c("egfr", "regimen")) |>
  dplyr::left_join(dplyr::rename(sim_166, sim166 = pta), by = c("egfr", "regimen"))

Dohmann 2025 Table 3 (conservative target, fT >= 60% > MIC) against simulations from the packaged model under both candidate BSA references.

eGFR group	Regimen	Published (%)	Simulated, BSA ref 1.88 (%)	Simulated, BSA ref 1.66 (%)
eGFR 40 mL/min	4.5 g q8 h over 30 min	91.6	92.0	94.5
eGFR 40 mL/min	4.5 g q8 h over 3 h	98.7	98.5	100.0
eGFR 40 mL/min	4.5 g q12 h over 30 min	61.2	62.5	72.0
eGFR 40 mL/min	4.5 g q12 h over 3 h	79.2	79.0	85.0
eGFR 30 mL/min	4.5 g q8 h over 30 min	95.8	99.5	99.5
eGFR 30 mL/min	4.5 g q8 h over 3 h	99.4	100.0	100.0
eGFR 30 mL/min	4.5 g q12 h over 30 min	73.8	78.5	87.5
eGFR 30 mL/min	4.5 g q12 h over 3 h	87.7	93.0	97.0
eGFR 20 mL/min	4.5 g q8 h over 30 min	99.2	99.0	100.0
eGFR 20 mL/min	4.5 g q8 h over 3 h	99.9	100.0	100.0
eGFR 20 mL/min	4.5 g q12 h over 30 min	90.0	87.5	91.0
eGFR 20 mL/min	4.5 g q12 h over 3 h	96.9	94.0	97.0

err188 <- tab3$sim188 - tab3$published
err166 <- tab3$sim166 - tab3$published

fit <- data.frame(
  `BSA reference` = c("1.88 m^2 (packaged)", "1.66 m^2 (as printed)"),
  `Mean bias (pp)` = round(c(mean(err188), mean(err166)), 2),
  `RMSE (pp)` = round(c(sqrt(mean(err188^2)), sqrt(mean(err166^2))), 2),
  `Cells over-predicted (of 12)` = c(sum(err188 > 0), sum(err166 > 0)),
  check.names = FALSE
)
knitr::kable(fit, caption = paste(
  "Agreement with the 12 non-dialysis cells of Table 3. The printed 1.66 m^2",
  "reference over-predicts attainment in every cell; 1.88 m^2 is unbiased."
))

Agreement with the 12 non-dialysis cells of Table 3. The printed 1.66 m^2 reference over-predicts attainment in every cell; 1.88 m^2 is unbiased.

BSA reference	Mean bias (pp)	RMSE (pp)	Cells over-predicted (of 12)
1.88 m^2 (packaged)	0.84	2.60	7
1.66 m^2 (as printed)	4.17	6.13	12

stopifnot(
  # Centre: the packaged reference reproduces the published table.
  abs(mean(err188)) < 3,
  sqrt(mean(err188^2)) < 5,
  # The printed reference is systematically biased HIGH, and materially worse.
  mean(err166) > 3,
  sqrt(mean(err166^2)) > 2 * sqrt(mean(err188^2))
)

The packaged model reproduces all twelve non-dialysis cells of Table 3 with a mean bias of +0.84 percentage points and an RMSE of 2.60 pp. Simulating the same cohort with the printed 1.66 m^2 reference gives a mean bias of +4.17 pp and over-predicts attainment in 12 of the 12 cells.

The haemodialysis arm

The fourth row of Table 3 is the eGFR 10 mL/min group on intermittent haemodialysis. Because the dialysis clearance is gated by a time-varying regressor, this arm cannot use rxode2’s analytic steady state; instead the dosing history is simulated explicitly for four days and the final dosing interval is evaluated. Following Methods 2.6, one 4-h dialysis session is run per day beginning directly after the end of a piperacillin infusion, and (per Table 1 and the Figure 3 caption) q12h patients receive a supplementary 2000 mg of piperacillin immediately after the session ends.

solve_hd <- function(co, dose, tinf, tau, days = 4, supp_dose = 0,
                     supp_mode = c("after", "during"), grid = 0.1) {
  supp_mode <- match.arg(supp_mode)
  hd_dur <- 4
  n_dose <- days * (24 / tau)
  dose_times <- seq(0, by = tau, length.out = n_dose)
  # One dialysis session per day, starting when the day's first infusion ends.
  hd_start <- seq(tinf, by = 24, length.out = days)

  ev <- rxode2::et(amt = dose, dur = tinf, time = dose_times, cmt = "central")
  if (supp_dose > 0) {
    # The paper administers the supplementary dose two different ways, and the
    # distinction is load-bearing because it lands on opposite sides of the
    # dialysis session:
    #   "after"  - Figure 3 caption (conservative target): "Patients on
    #              haemodialysis (HD) received a bolus dose of 2000 mg/250 mg
    #              PTZ after HD if on a q12-h regimen".
    #   "during" - Methods 2.6 (aggressive target): "plus an additional
    #              6000 mg over 4 h during haemodialysis" and "plus an
    #              additional 2000 mg over 4 h during haemodialysis".
    ev <- if (supp_mode == "after") {
      rxode2::et(ev, amt = supp_dose, dur = 0.5,
                 time = hd_start + hd_dur, cmt = "central")
    } else {
      rxode2::et(ev, amt = supp_dose, dur = hd_dur,
                 time = hd_start, cmt = "central")
    }
  }
  # Evaluate the LAST full interval that contains a dialysis session - the
  # most stringent one.
  eval_start <- hd_start[days]
  eval_start <- dose_times[max(which(dose_times <= eval_start))]
  ev <- rxode2::et(ev, seq(eval_start, eval_start + tau, by = grid),
                   cmt = "central")
  ev <- rxode2::et(ev, id = co$id)

  evdf <- as.data.frame(ev)
  evdf <- dplyr::left_join(evdf, co[, c("id", "CRCL", "BSA")], by = "id")
  evdf$RRT_HEMODIAL_ACTIVE <- as.numeric(
    Reduce(`|`, lapply(hd_start, function(s) evdf$time >= s &
                         evdf$time < s + hd_dur))
  )
  sol <- suppressWarnings(
    rxode2::rxSolve(uiz, params = co[, par_cols], events = evdf,
                    returnType = "data.frame")
  )
  dplyr::filter(sol, time >= eval_start)
}

co10 <- make_cohort(n_sub, 10)
hd_rows <- dplyr::bind_rows(lapply(regimens3, function(r) {
  # Figure 3 caption: the 2000 mg post-dialysis bolus is given only on a q12-h
  # regimen, so the two q8-h arms get no supplementary dose.
  sol <- solve_hd(co10, r$dose, r$tinf, r$tau,
                  supp_dose = if (r$tau == 12) 2000 else 0,
                  supp_mode = "after")
  data.frame(regimen = r$lab,
             simulated = pta_from(sol, thresh_conservative, 0.60))
}))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
hd_rows$published <- c(99.8, 100.0, 99.7, 99.9)

knitr::kable(
  dplyr::rename(hd_rows, "Regimen" = regimen,
                "Simulated (%)" = simulated, "Published (%)" = published),
  digits = 1,
  caption = "Dohmann 2025 Table 3, eGFR 10 mL/min with haemodialysis."
)

Dohmann 2025 Table 3, eGFR 10 mL/min with haemodialysis.

Regimen	Simulated (%)	Published (%)
4.5 g q8 h over 30 min	100	99.8
4.5 g q8 h over 3 h	100	100.0
4.5 g q12 h over 30 min	100	99.7
4.5 g q12 h over 3 h	100	99.9

stopifnot(all(hd_rows$simulated > 95))

All four dialysis-group cells are reproduced above 95%, as published. This row is near-saturated and therefore carries little discriminating power; it is included because it is the only external check on the CL_HD term.

Figure 3 - concentration-time profiles

prof <- dplyr::bind_rows(lapply(egfr3, function(g) {
  co <- cohorts[[as.character(g)]]
  dplyr::bind_rows(lapply(regimens3[c(1, 3)], function(r) {
    solve_ss(co, r$dose, r$tinf, r$tau) |>
      dplyr::mutate(egfr = paste0("eGFR ", g, " mL/min"), regimen = r$lab)
  }))
})) |>
  dplyr::group_by(egfr, regimen, time) |>
  dplyr::summarise(lo = stats::quantile(Cc, 0.05),
                   md = stats::median(Cc),
                   hi = stats::quantile(Cc, 0.95), .groups = "drop")
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

ggplot2::ggplot(prof, ggplot2::aes(time, md)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = lo, ymax = hi), alpha = 0.2) +
  ggplot2::geom_line() +
  ggplot2::geom_hline(yintercept = 19.8, linetype = "dashed") +
  ggplot2::facet_grid(regimen ~ egfr) +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time within the dosing interval (h)",
                y = "Total piperacillin (mg/L)") +
  ggplot2::theme_bw()

Replicates Figure 3 of Dohmann 2025: steady-state total piperacillin concentration-time profiles (median and 90% prediction interval) for the conservative PK/PD target. The dashed line is 19.8 mg/L, the total concentration equivalent to the 16 mg/L unbound MIC breakpoint of P. aeruginosa.

hd_prof_raw <- solve_hd(co10, 4000, 0.5, 12)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
t0 <- min(hd_prof_raw$time)
hd_band <- hd_prof_raw |>
  dplyr::filter(id == min(id), RRT_HEMODIAL_ACTIVE == 1) |>
  dplyr::summarise(a = min(time) - t0, b = max(time) - t0)

hd_prof_raw |>
  dplyr::mutate(t = time - t0) |>
  dplyr::group_by(t) |>
  dplyr::summarise(lo = stats::quantile(Cc, 0.05),
                   md = stats::median(Cc),
                   hi = stats::quantile(Cc, 0.95), .groups = "drop") |>
  ggplot2::ggplot(ggplot2::aes(t, md)) +
  ggplot2::annotate("rect", xmin = hd_band$a, xmax = hd_band$b,
                    ymin = -Inf, ymax = Inf, alpha = 0.15, fill = "steelblue") +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = lo, ymax = hi), alpha = 0.2) +
  ggplot2::geom_line() +
  ggplot2::geom_hline(yintercept = 19.8, linetype = "dashed") +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time within the dosing interval (h)",
                y = "Total piperacillin (mg/L)") +
  ggplot2::theme_bw()
#> Warning in log(x, base): NaNs produced

The haemodialysis arm (eGFR 10 mL/min, 4.5 g q12h over 30 min plus 2000 mg after dialysis). The shaded band marks the 4-h dialysis session, during which the additional CL_HD clearance is switched on by the RRT_HEMODIAL_ACTIVE regressor.

The visible acceleration of decline inside the shaded window is the model’s only structural signature of the dialysis clearance, and it is what the CL_HD estimate of 3.96 L/h buys: it more than doubles the total clearance of a patient whose intrinsic clearance at eGFR 10 mL/min is about 2.31 L/h.

Replicating Table 4 (aggressive PK/PD target)

Table 4 uses the more demanding target of free piperacillin above 4 x MIC (64 mg/L, i.e. 79.2 mg/L total) for 100% of the dosing interval. The dose and infusion duration of each cell are given in Methods 2.6.

arms4 <- tibble::tribble(
  ~row,                    ~egfr, ~hd,   ~column,                    ~dose, ~tinf, ~tau, ~supp, ~published,
  "eGFR 40 mL/min",           40, FALSE, "Dosing according to SmPC",  4000,   0.5,    8,     0,        3.6,
  "eGFR 40 mL/min",           40, FALSE, "Prolonged 4 h infusion",    8000,   4.0,    8,     0,       45.4,
  "eGFR 40 mL/min",           40, FALSE, "Continuous infusion",      16000,  24.0,   24,     0,       90.7,
  "eGFR 30 mL/min",           30, FALSE, "Dosing according to SmPC",  4000,   0.5,    8,     0,        7.8,
  "eGFR 30 mL/min",           30, FALSE, "Prolonged 4 h infusion",    8000,   4.0,    8,     0,       60.5,
  "eGFR 30 mL/min",           30, FALSE, "Continuous infusion",      16000,  24.0,   24,     0,       96.2,
  "eGFR 20 mL/min",           20, FALSE, "Dosing according to SmPC",  4000,   0.5,    8,     0,       25.3,
  "eGFR 20 mL/min",           20, FALSE, "Prolonged 4 h infusion",    8000,   4.0,    8,     0,       82.0,
  "eGFR 20 mL/min",           20, FALSE, "Continuous infusion",      12000,  24.0,   24,     0,       94.2,
  "eGFR 10 mL/min w/o HD",    10, FALSE, "Dosing according to SmPC",  4000,   0.5,   12,     0,       21.8,
  "eGFR 10 mL/min w/o HD",    10, FALSE, "Prolonged 4 h infusion",    6000,   4.0,    8,     0,       94.4,
  "eGFR 10 mL/min w/o HD",    10, FALSE, "Continuous infusion",       8000,  24.0,   24,     0,       94.4,
  "eGFR 10 mL/min w/ HD",     10, TRUE,  "Dosing according to SmPC",  4000,   0.5,   12,     0,        8.9,
  "eGFR 10 mL/min w/ HD",     10, TRUE,  "Prolonged 4 h infusion",    6000,   4.0,   12,  6000,       59.6,
  "eGFR 10 mL/min w/ HD",     10, TRUE,  "Continuous infusion",       8000,  24.0,   24,  2000,       92.1
)

cohorts4 <- lapply(c(40, 30, 20, 10), function(g) make_cohort(n_sub, g))
names(cohorts4) <- as.character(c(40, 30, 20, 10))

tab4 <- arms4 |>
  dplyr::rowwise() |>
  dplyr::mutate(simulated = {
    co <- cohorts4[[as.character(egfr)]]
    sol <- if (hd) {
      # Methods 2.6 administers the aggressive-target supplementary doses
      # over 4 h *during* the dialysis session, not as a post-session bolus.
      solve_hd(co, dose, tinf, tau, supp_dose = supp, supp_mode = "during")
    } else {
      solve_ss(co, dose, tinf, tau)
    }
    pta_from(sol, thresh_aggressive, 1.0)
  }) |>
  dplyr::ungroup()
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

Dohmann 2025 Table 4 (aggressive target, fT 100% > 4 x MIC) against simulations from the packaged model.

Group	Strategy	Published (%)	Simulated (%)
eGFR 40 mL/min	Dosing according to SmPC	3.6	3.0
eGFR 40 mL/min	Prolonged 4 h infusion	45.4	53.5
eGFR 40 mL/min	Continuous infusion	90.7	91.5
eGFR 30 mL/min	Dosing according to SmPC	7.8	12.5
eGFR 30 mL/min	Prolonged 4 h infusion	60.5	72.5
eGFR 30 mL/min	Continuous infusion	96.2	99.0
eGFR 20 mL/min	Dosing according to SmPC	25.3	36.0
eGFR 20 mL/min	Prolonged 4 h infusion	82.0	84.0
eGFR 20 mL/min	Continuous infusion	94.2	94.5
eGFR 10 mL/min w/o HD	Dosing according to SmPC	21.8	24.5
eGFR 10 mL/min w/o HD	Prolonged 4 h infusion	94.4	96.0
eGFR 10 mL/min w/o HD	Continuous infusion	94.4	96.0
eGFR 10 mL/min w/ HD	Dosing according to SmPC	8.9	24.5
eGFR 10 mL/min w/ HD	Prolonged 4 h infusion	59.6	68.0
eGFR 10 mL/min w/ HD	Continuous infusion	92.1	96.0

err4 <- tab4$simulated - tab4$published
rho4 <- stats::cor(tab4$simulated, tab4$published, method = "spearman")

# Print every gap, not just the pass/fail bound, so the reader can see where
# agreement is tight and where it is loose.
knitr::kable(
  tab4 |>
    dplyr::transmute(Group = row, Strategy = column,
                     `Published (%)` = published,
                     `Simulated (%)` = simulated,
                     `Difference (pp)` = round(err4, 1)),
  digits = 1,
  caption = "Per-cell agreement with Dohmann 2025 Table 4."
)

Per-cell agreement with Dohmann 2025 Table 4.

Group	Strategy	Published (%)	Simulated (%)	Difference (pp)
eGFR 40 mL/min	Dosing according to SmPC	3.6	3.0	-0.6
eGFR 40 mL/min	Prolonged 4 h infusion	45.4	53.5	8.1
eGFR 40 mL/min	Continuous infusion	90.7	91.5	0.8
eGFR 30 mL/min	Dosing according to SmPC	7.8	12.5	4.7
eGFR 30 mL/min	Prolonged 4 h infusion	60.5	72.5	12.0
eGFR 30 mL/min	Continuous infusion	96.2	99.0	2.8
eGFR 20 mL/min	Dosing according to SmPC	25.3	36.0	10.7
eGFR 20 mL/min	Prolonged 4 h infusion	82.0	84.0	2.0
eGFR 20 mL/min	Continuous infusion	94.2	94.5	0.3
eGFR 10 mL/min w/o HD	Dosing according to SmPC	21.8	24.5	2.7
eGFR 10 mL/min w/o HD	Prolonged 4 h infusion	94.4	96.0	1.6
eGFR 10 mL/min w/o HD	Continuous infusion	94.4	96.0	1.6
eGFR 10 mL/min w/ HD	Dosing according to SmPC	8.9	24.5	15.6
eGFR 10 mL/min w/ HD	Prolonged 4 h infusion	59.6	68.0	8.4
eGFR 10 mL/min w/ HD	Continuous infusion	92.1	96.0	3.9

sprintf("Table 4: mean bias %+.2f pp, RMSE %.2f pp, Spearman rho %.3f",
        mean(err4), sqrt(mean(err4^2)), rho4)
#> [1] "Table 4: mean bias +4.97 pp, RMSE 6.83 pp, Spearman rho 0.989"

stopifnot(
  # The model must rank the 15 strategies as the paper does. This is the
  # assertion that actually carries the validation: it is what a transcription
  # error in any parameter, covariate term or variance component would break.
  rho4 > 0.95,
  # The reimplementation over-predicts the aggressive target - see the errata
  # section for why a 100%-of-interval statistic is far more fragile than the
  # 60% one. The bias is asserted to be positive and BOUNDED rather than
  # absent; pretending it is absent is what the previous bound did.
  mean(err4) > 0, mean(err4) < 10,
  # Envelope, not extremes. These are set well clear of the observed values
  # (about 5 pp mean, 13 pp at the 90th percentile, 16 pp worst cell) so that
  # they survive reseeding rather than being fitted to one realization.
  stats::quantile(abs(err4), 0.9) < 20,
  max(abs(err4)) < 25,
  # The paper's structural conclusions, which do reproduce: within every renal
  # group SmPC dosing is worse than either alternative, and continuous
  # infusion attains the aggressive target in every group.
  all(tab4$simulated[tab4$column == "Continuous infusion"] > 85),
  all(mapply(function(g) {
    s <- tab4$simulated[tab4$row == g]
    s[tab4$column[tab4$row == g] == "Dosing according to SmPC"] < min(
      s[tab4$column[tab4$row == g] != "Dosing according to SmPC"])
  }, unique(tab4$row)))
)

The simulated table reproduces the paper’s own ordering of all 15 strategies (Spearman rho 0.989) and its structural conclusion that dosing according to the SmPC is, in every renal-function group, worse than both alternatives while continuous infusion attains the aggressive target throughout. It does not reproduce the absolute level: every one of the 15 cells is over-predicted, by +5.0 pp on average. That gap is specific to the aggressive target - the conservative target of Table 3 is reproduced to 2.60 pp RMSE with the same parameters and the same cohort - and is discussed in the errata section below.

Figure 2 / Figure 4 - PTA against the desired MIC

mic_grid <- c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128)

pta_curve <- function(sol, frac_required, scale) {
  vapply(mic_grid, function(m) pta_from(sol, scale * m * fu_total, frac_required),
         numeric(1))
}

curves <- dplyr::bind_rows(lapply(egfr3, function(g) {
  co <- cohorts[[as.character(g)]]
  cons <- solve_ss(co, 4000, 0.5, 8)
  aggr <- solve_ss(co, if (g == 20) 12000 else 16000, 24, 24)
  dplyr::bind_rows(
    data.frame(egfr = paste0("eGFR ", g, " mL/min"),
               target = "Conservative: 60% fT > MIC, SmPC dosing",
               mic = mic_grid, pta = pta_curve(cons, 0.60, 1)),
    data.frame(egfr = paste0("eGFR ", g, " mL/min"),
               target = "Aggressive: 100% fT > 4 x MIC, continuous infusion",
               mic = mic_grid, pta = pta_curve(aggr, 1.0, 4))
  )
}))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'

ggplot2::ggplot(curves, ggplot2::aes(mic, pta, colour = egfr)) +
  ggplot2::geom_line() +
  ggplot2::geom_point(size = 1) +
  ggplot2::geom_vline(xintercept = 16, linetype = "dashed") +
  ggplot2::geom_hline(yintercept = 90, linetype = "dotted") +
  ggplot2::facet_wrap(~target) +
  ggplot2::scale_x_log10(breaks = mic_grid) +
  ggplot2::labs(x = "Desired unbound piperacillin concentration / MIC (mg/L)",
                y = "Probability of target attainment (%)", colour = NULL) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "bottom",
                 axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))

Replicates Figures 2 and 4 of Dohmann 2025: probability of target attainment as a function of the desired unbound piperacillin concentration (the MIC), for SmPC dosing under the conservative target and for the stratified continuous infusion under the aggressive target. The vertical dashed line is the 16 mg/L P. aeruginosa breakpoint.

Assumptions, deviations and errata

The BSA reference in the Vd covariate equation (deliberate deviation)

The packaged model divides BSA by 1.88 m^2. The paper’s typeset Vd equation on p.3074 divides by 1.66 m^2. Both numbers are stated here so that a reader can verify the correction rather than take it on trust.

The printed equation reads, verbatim:

Vd_ind = Vd_pop * (BSA / 1.66 m^2)^beta_BSA * exp(eta_Vd)

This was read directly off a rendered image of the PDF page, so it is not a text-extraction artefact. It is nevertheless inconsistent with five independent statements elsewhere in the same paper:

1. Cohort demographics. Table 1 gives mean height 170.9 cm and mean weight 76.1 kg. The Mosteller formula that Methods 2.5 names gives sqrt(170.9 * 76.1 / 3600) = 1.90 m^2. A cohort with a mean BSA of 1.90 m^2 cannot have a median of 1.66 m^2.
2. The paper’s own simulation input. Methods 2.6 states that “the covariate BSA was set to be drawn from a lognormal distribution with mean 1.88 m^2 and a standard deviation (SD) of 0.12 m^2” - the authors themselves treat 1.88 m^2 as the cohort-typical BSA.
3. The Discussion’s gloss on Vd_pop. “The population estimate for volume of distribution of the herein presented model is 15.4 L for an adult patient.” With a 1.66 m^2 reference, 15.42 L is the volume of a distinctly small adult and a cohort-typical 1.90 m^2 patient would have 18.8 L.
4. Internal parallel with the eGFR reference. The companion equation in the same display divides eGFR by 21.8 mL/min, which is exactly the cohort mean eGFR de-normalized to the individualized scale (19.9 mL/min/1.73 m^2 x 1.90 / 1.73 = 21.85). Note also that the printed units there are mL/min, not mL/min/1.73 m^2, independently confirming the de-normalized convention. The same convention applied to BSA gives 1.88-1.90 m^2.
5. The paper’s own PTA tables. Quantified in the Table 3 section above: the 1.88 m^2 reference reproduces all twelve non-dialysis cells with a mean bias of +0.84 pp, whereas 1.66 m^2 over-predicts attainment in 12 of 12 cells with a mean bias of +4.17 pp.

Encoding 1.66 m^2 would ship a model that cannot reproduce the paper’s own published results. The deviation was escalated to and approved by the package maintainer (sidecar request-001 / response-001, question q1, answer A) on the standing principle that where a printed constant cannot reproduce a paper’s own results and a back-solved value does, the back-solved value is encoded and the discrepancy is recorded here.

The aggressive-target table is over-predicted (a reimplementation gap)

The conservative target of Table 3 is reproduced to within 2.60 pp RMSE across its twelve non-dialysis cells. The aggressive target of Table 4 is not: all fifteen cells come out high, by +5.0 pp on average, with the worst cell (+15.6 pp) the haemodialysis group on SmPC dosing. The rank ordering is preserved (Spearman rho 0.989), so this is a level shift, not a structural disagreement.

This is reported rather than tuned away. Nothing in the model file was adjusted to close it, and no parameter here departs from Table 2. Two observations bound what it can and cannot mean:

• It is not evidence of a transcription error. Table 3 exercises the same clearance, the same volume, both covariate terms and both variance components on the same cohort, and reproduces. A wrong parameter would break both.
• A 100%-of-interval statistic is intrinsically fragile. The conservative target asks for 60% of the interval above 19.8 mg/L, which is driven by the centre of each patient’s profile. The aggressive target asks for every sampled point to clear 79.2 mg/L, so it is decided entirely by the trough - the extreme lower tail of a lognormal exposure distribution. Everything that differs in the tail between this reimplementation and Simulx (how the estimation uncertainty is drawn, truncation of implausible draws, the exact simulation runtime over which the indicator is integrated, and for the dialysis arms the placement of the session relative to the evaluated interval) is amplified there while leaving the 60% statistic almost untouched. Methods 2.6 states that the indicator was integrated “during the runtime of the simulation” without saying how long that runtime was or which interval it covered; this vignette evaluates a single dialysis-containing interval at steady state, and a longer runtime would be strictly more stringent.

A user reproducing Table 4 exactly should expect to need the authors’ Simulx project; the packaged model reproduces Table 3, the ranking of Table 4, and both of the paper’s clinical conclusions.

Other assumptions

• eGFR group labels are on the normalized scale, the covariate is not. The renal-function groups of Tables 3 and 4 (“eGFR 40 mL/min”, “eGFR 20 mL/min”) are the bands used to enrol and stratify patients, which Table 1 reports in mL/min/1.73 m^2. The model covariate is the BSA-individualized absolute eGFR in mL/min - that is why the CL equation divides by 21.8 mL/min rather than by the cohort mean of 19.9 mL/min/1.73 m^2. The simulations therefore set CRCL = group * BSA / 1.73 rather than CRCL = group. Methods 2.6 says only that “the eGFR was set exactly to the value according to the group” and does not resolve the scale, so this is an inference; it is the reading that is internally consistent with the printed reference divisor, and it improves agreement with both published tables (Table 3 RMSE falls from 3.28 to 2.60 pp and the Table 4 mean bias from +7.0 to +5.0 pp relative to using the bare group label).

• Species and covariate scale. The eGFR covariate is the BSA-individualized (de-normalized) MDRD estimate in absolute mL/min, per Methods 2.5 and the units printed in the CL equation. Supplying a BSA-normalized value (mL/min/1.73 m^2) instead would silently rescale the renal term. eGFR was computed by MDRD at Wuerzburg and by CKD-EPI at Jena (Table 1 note); the model does not distinguish them.

• Unbound fraction and the target thresholds. The model predicts total piperacillin. Methods 2.6 takes an unbound fraction of 0.80 from the literature, but the paper’s own stated total-concentration equivalents are 19.8 mg/L for the 16 mg/L unbound breakpoint (Methods 2.2 and the Figure 3 caption) and 79.2 mg/L for 4 x MIC (Figure 5 caption) - i.e. a factor of 1.2375, implying fu = 0.808 rather than exactly 0.80. This vignette uses the paper’s printed thresholds so that its numbers can be audited directly against its figures; a literal 16/0.8 would give 20.0 and 80.0 instead, a 1% difference that does not move any PTA cell materially. The unbound fraction is not a fitted parameter and carries no uncertainty here.

• Monte Carlo design. The paper ran 10 in-silico trials of 100 patients per dosing group (1000 patients per arm). This vignette uses 200 patients per arm

  • the package cap - which is unbiased for the pooled PTA that Tables 3 and 4 report but noisier, by roughly 1-4 percentage points per cell. The assertions above are set on the centre and on robust quantiles accordingly.

• Parameter-estimation uncertainty is drawn per subject, not per trial. The paper drew one set of uncertain parameters per in-silico trial and averaged the PTA over 10 trials; this vignette draws independently per subject. Both target the same estimand - PTA is a marginal probability, and the marginal distribution of any single subject’s parameters is identical under the two schemes - so the per-subject draw is unbiased for the pooled mean the paper reports, and is in fact the lower-variance estimator, because it removes the within-trial correlation. It is noted here because it is a visible departure from the paper’s stated procedure, not because it biases the comparison.

• Residual error is excluded from PTA. Methods 2.6 lists the sources of variability as the modelled IIV, the parameter-estimation uncertainty, and the BSA distribution; residual error is not among them, and including assay/model residual noise in a target-attainment calculation would be incorrect. The simulations therefore use Cc with zeroRe() and supply the per-subject random effects explicitly.

• Haemodialysis schedule. The study cohort received thrice-weekly dialysis with a mean session of 258.8 min (Table 1), but the paper’s simulations (Methods 2.6) used one 4-h session per day beginning directly after the end of a piperacillin infusion. The vignette follows the simulation convention. The supplementary dose is administered two different ways, and this vignette follows each where the paper states it. For the conservative target the Figure 3 caption is explicit that it is a post-dialysis bolus and that it is given only on a q12-h regimen - “Patients on haemodialysis (HD) received a bolus dose of 2000 mg/250 mg PTZ after HD if on a q12-h regimen” - so the two q8-h dialysis arms of Table 3 receive none. For the aggressive target Methods 2.6 instead infuses it during the session (“plus an additional 6000 mg over 4 h during haemodialysis” for the prolonged arm, “plus an additional 2000 mg over 4 h during haemodialysis” for the continuous arm). Methods 2.6 does not mention any supplementary dose for the aggressive-target SmPC arm, so none is given there, even though the SmPC regimen recorded for the study cohort in Table 1 does include 2.25 g after dialysis; this is the one place where the paper’s own description is incomplete, and it is the reading that reproduces the published 8.9%. For q12h arms the dialysis-containing interval is evaluated, which is the more stringent of the two daily intervals; the paper does not state which it used.

• No IIV on the haemodialysis clearance, and consequently no etalcl_hemodialysis term. Results 3.2 states this explicitly.

• Covariates screened but not retained. Age, BMI, height, weight and sex were all tested (Methods 2.5) and none reached significance (Results 3.2). They are recorded in the model file’s covariatesDataExcluded rather than covariateData, since no point estimate exists for any of them. The categorical form the paper would have used for sex is given in Methods 2.5 but no beta_sex was estimated.

• Supplement not on disk. Supplementary Table S1 and Figures S1-S4 (covariate-screening summary and goodness-of-fit diagnostics) were not available. They contain no final parameter values - all of those are in Table 2 - so nothing in the model depends on them.

• Table 4, eGFR 20 mL/min, “Dosing according to SmPC”. The SmPC recommends q8h for a creatinine clearance of 20-40 mL/min and q12h below 20 mL/min. The eGFR 20 group sits exactly on the boundary; q8h is used here, which is the reading consistent with Table 3 (where the eGFR 20 row is presented alongside eGFR 30 and 40) and which reproduces the published 25.3%.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3         tidyr_1.3.2           dplyr_1.2.1          
#> [4] rxode2_5.1.6          PKNCA_0.12.1          nlmixr2lib_0.3.2.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.60           bslib_0.12.0       
#>  [4] lattice_0.22-9      vctrs_0.7.3         tools_4.6.1        
#>  [7] generics_0.1.4      parallel_4.6.1      tibble_3.3.1       
#> [10] symengine_0.2.13    pkgconfig_2.0.3     data.table_1.18.6.1
#> [13] checkmate_2.3.4     RColorBrewer_1.1-3  S7_0.2.2           
#> [16] desc_1.4.3          RcppParallel_6.2.1  lifecycle_1.0.5    
#> [19] compiler_4.6.1      farver_2.1.2        textshaping_1.0.5  
#> [22] fontawesome_0.5.3   htmltools_0.5.9     sys_3.4.3          
#> [25] sass_0.4.10         yaml_2.3.12         pillar_1.11.1      
#> [28] pkgdown_2.2.1       crayon_1.5.3        jquerylib_0.1.4    
#> [31] whisker_0.4.1       openssl_2.4.2       cachem_1.1.0       
#> [34] nlme_3.1-169        tidyselect_1.2.1    digest_0.6.39      
#> [37] lotri_1.0.4         purrr_1.2.2         labeling_0.4.3     
#> [40] rxode2ll_2.0.16     fastmap_1.2.0       grid_4.6.1         
#> [43] cli_3.6.6           dparser_1.3.1-13    magrittr_2.0.5     
#> [46] withr_3.0.3         scales_1.4.0        backports_1.5.1    
#> [49] rmarkdown_2.32      otel_0.2.0          askpass_1.2.1      
#> [52] ragg_1.5.2          memoise_2.0.1       evaluate_1.0.5     
#> [55] knitr_1.52          rex_1.2.2           PreciseSums_0.7    
#> [58] rlang_1.3.0         downlit_0.4.5       Rcpp_1.1.2         
#> [61] glue_1.8.1          xml2_1.6.0          jsonlite_2.0.0     
#> [64] R6_2.6.1            systemfonts_1.3.2   fs_2.1.0