Tigecycline (Broeker 2018) • nlmixr2lib

library(nlmixr2lib)
library(rxode2)
library(PKNCA)
library(dplyr)
library(ggplot2)

This vignette validates Broeker_2018_tigecycline, the two-compartment population PK model for intravenous tigecycline in critically ill adults on continuous renal replacement therapy (CRRT) published by Broeker and colleagues (doi:10.1186/s13054-018-2278-4).

The distinguishing feature of the analysis is that plasma and CRRT effluent concentrations were fitted simultaneously. Effluent is what separates the dialysis clearance from the physiological body clearance: without it, a model fitted to plasma alone can only report the total. Because the cohort spans two CRRT modalities with different solute-transport mechanisms – continuous venovenous hemodialysis (CVVHD, purely diffusive) and continuous venovenous hemodiafiltration (CVVHDF, diffusive plus convective) – the paper estimates a separate dialysis clearance for each while sharing one body-PK model.

The clinical conclusion is the negative one: despite tigecycline passing the membrane freely (saturation coefficients of 0.79 and 0.90), CRRT removes only about a tenth of the dose, because tigecycline’s volume of distribution is large enough that little drug is in plasma to be dialysed. No dose adjustment is needed.

mod <- readModelDb("Broeker_2018_tigecycline")()
zmod <- rxode2::zeroRe(mod)
mod
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>                     lcl                     lvc                     lvp 
#>               2.9069011               4.0724397               5.0369526 
#>                      lq  lcl_hemodialysis_cvvhd lcl_hemodialysis_cvvhdf 
#>               4.0324692               0.5247285               0.9969486 
#>              e_tbili_cl                  propSd        propSd_Ceffluent 
#>              -0.2900000               0.1690000               0.4060000 
#> 
#> Omega ($omega): 
#>                             etalcl   etalvc    etalq etalcl_hemodialysis_cvvhd
#> etalcl                    0.190096 0.000000 0.000000                  0.000000
#> etalvc                    0.000000 1.229881 0.000000                  0.000000
#> etalq                     0.000000 0.000000 0.174724                  0.000000
#> etalcl_hemodialysis_cvvhd 0.000000 0.000000 0.000000                  0.189225
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1          central
#> 2                  2      peripheral1
#>  ── Multiple Endpoint Model ($multipleEndpoint): ──  
#>        variable                      cmt                      dvid*
#> 1        Cc ~ …        cmt='Cc' or cmt=3        dvid='Cc' or dvid=1
#> 2 Ceffluent ~ … cmt='Ceffluent' or cmt=4 dvid='Ceffluent' or dvid=2
#>   * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc
#> 
#>  ── μ-referencing ($muRefTable): ──  
#>                    theta                       eta level
#> 1                    lcl                    etalcl    id
#> 2 lcl_hemodialysis_cvvhd etalcl_hemodialysis_cvvhd    id
#> 3                    lvc                    etalvc    id
#> 4                     lq                     etalq    id
#> 
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     compartmentData <- list(central = list(analyte = "tigecycline", 
#>         units = "mg", specimen = "plasma", verified = TRUE), 
#>         peripheral1 = list(analyte = "tigecycline", units = "mg", 
#>             specimen = "plasma", verified = TRUE))
#>     covariateData <- list(TBILI = list(description = "Total serum bilirubin", 
#>         units = "umol/L", type = "continuous", reference_category = NULL, 
#>         notes = "The ONLY covariate retained in the final model. Enters the body clearance as a median-normalised power term, Broeker 2018 Table 2 row heading: 'Clearance (L/h) = theta_1 x (bilirubin/2.3)^theta_2' with theta_1 = 18.3 L/h and theta_2 = -0.29. UNITS: the paper reports bilirubin in the US convention mg/dL (Table 1 column header 'Bilirubin (mg/dL)') and the reference value 2.3 is therefore 2.3 mg/dL, stated in the Results as 'normalized by the population median of bilirubin, 2.3 mg/dL' and confirmed by the Table 1 median row. This column carries the register's canonical SI umol/L, so model() converts back inline with TBILI / 17.1 before forming the ratio; equivalently the reference is 39.3 umol/L. The exponent is NEGATIVE, so higher bilirubin means LOWER clearance -- Results: 'lower bilirubin concentrations corresponded to higher clearances'. The effect is large because the cohort spans two orders of magnitude of bilirubin (0.7 to 43.3 mg/dL, Table 1) in a population with four liver-failure or cirrhosis patients and two liver transplants: the Results quantify it as 'Individual clearance values varied from 9.3 L/h (10th percentile) to 19.1 L/h (90th percentile) depending on the bilirubin concentration (24 mg/dL to 1.8 mg/dL)', which the encoded term reproduces as 18.3 * (24/2.3)^-0.29 = 9.3 and 18.3 * (1.8/2.3)^-0.29 = 19.7. Including it dropped the objective function by 5.71 (p = 0.017) and cut the unexplained interindividual variability on clearance from 58.6% to 43.6%, so the 43.6% in ini() is the POST-covariate value. The authors read bilirubin as a marker of hepatic function, tigecycline being predominantly biliary/faecally eliminated. Treated as a time-fixed baseline value (a single day-4 laboratory value per patient).", 
#>         source_name = "bilirubin"), RRT_CRRT_ACTIVE = list(description = "CRRT-active indicator (1 while the continuous renal replacement circuit is running, 0 otherwise)", 
#>         units = "(binary)", type = "binary", reference_category = "0 (circuit not running)", 
#>         notes = "Gates the dialysis clearance arm and, with it, the effluent observable. EVERY subject in Broeker 2018 was on CRRT continuously for the whole sampling period -- blood was drawn 'on day 4 of treatment with tigecycline after at least 24 h on CRRT' (Methods, Sampling and drug analysis) and the study has no off-CRRT subgroup -- so this column is identically 1 throughout the source analysis and the paper estimates no off-CRRT clearance. It is carried anyway, per the RRT_CRRT_EFFLUENT_FLOW register entry's instruction to pair the effluent flow with an on/off gate, so that a downstream user can simulate circuit interruption or discontinuation; setting it to 0 leaves the body clearance alone, which is the paper's CLbody and NOT a validated off-CRRT clearance for this population. The ACTIVE rather than the STATUS member of the RRT_<modality>_<kind> family is used because the quantity is physically time-varying (filter changes interrupt the circuit), matching ButraguenoLaiseca_2022_piperacillin.", 
#>         source_name = "(not a source column; implicitly 1 for all subjects)"), 
#>         RRT_CVVHDF_STATUS = list(description = "CRRT modality indicator (1 = continuous venovenous hemodiafiltration, 0 = continuous venovenous hemodialysis)", 
#>             units = "(binary)", type = "binary", reference_category = "0 (CVVHD, the reference modality carrying 8 of the 11 subjects)", 
#>             notes = "Selects which of the two published dialysis clearances applies. Broeker 2018 Table 2 estimates them as two separate parameters -- 'Dialysis clearance CVVHD (L/h) 1.69' and 'Dialysis clearance CVVHDF (L/h) 2.71' -- over a cohort that shares one body-PK model: 'Eleven patients ... receiving either continuous veno-venous hemodialysis (CVVHD, n = 8) or hemodiafiltration (CVVHDF, n = 3)' (Abstract, Methods). Per-subject modality assignment is recoverable from Table 1, whose footnotes read 'a Continuous veno-venous hemodialysis (CVVHD); b continuous veno-venous hemodiafiltration (CVVHDF)': patients 1, 2, 3, 6, 8, 9, 10 and 11 are CVVHD and patients 4, 5 and 7 are CVVHDF. MEANINGFUL ONLY WHEN RRT_CRRT_ACTIVE = 1; the pair (RRT_CRRT_ACTIVE, RRT_CVVHDF_STATUS) is what distinguishes off-circuit from CVVHD, which a modality indicator alone cannot do. The two modalities differ physically in that CVVHDF adds convective transport through a post-dilution ultrafiltrate stream (QFil = 1 L/h, Methods) on top of CVVHD's purely diffusive dialysate stream, which is why its dialysis clearance is the larger of the two and why its mean saturation coefficient is higher (0.90 versus 0.79, Results). Only the CVVHD arm carries IIV: Table 2 prints a dash for the CVVHDF arm and the Results explain that 'an IIV for this method was not supported by the data (IIV tended to zero during estimation)', unsurprising at n = 3.", 
#>             source_name = "Table 1 footnote markers a / b"), 
#>         RRT_CRRT_EFFLUENT_FLOW = list(description = "Total effluent flow rate leaving the CRRT circuit (dialysate + ultrafiltrate)", 
#>             units = "mL/h", type = "continuous", reference_category = NULL, 
#>             notes = "Enters ONLY the effluent observation equation, never the clearance model -- Broeker 2018 estimates the dialysis clearance directly rather than as a sieving coefficient times a flow. Broeker 2018 Eq. 3 defines CLDial,CVVHD = QDial * (Ceff/CPla) and Eq. 4 defines CLDial,CVVHDF = (QDial + QFil) * (Ceff/CPla), where 'QDial represents the dialysate flow rate, Ceff represents the concentration of tigecycline in the effluent, and CPla represents the pre-filter plasma concentration'. Rearranging for the observable gives Ceff = CLDial * CPla / Qeff with Qeff = QDial for CVVHD and QDial + QFil for CVVHDF -- exactly the form used by ButraguenoLaiseca_2022_piperacillin, so this column is the single Qeff denominator for both modalities. VALUES: the Methods prescribe the dialysate flow by weight band, 'Blood flow and dialysate flow were adjusted to body weight (< 90 kg/> 90 kg; 100/120 mL/min and 2000/2500 mL/h, respectively)', and fix the CVVHDF ultrafiltration rate at 'QFil was 1 L/h'. The weight-band rule gives 2000 mL/h for all eight CVVHD patients (all under 90 kg, Table 1) and 3000 mL/h for a CVVHDF patient under 90 kg. That reading is confirmed ARITHMETICALLY by the paper's own derived statistic: the published CVVHDF saturation coefficient is 0.90 (Results) and 2.71 / 3.000 = 0.903. The CVVHD counterpart, 1.69 / 2.000 = 0.845, sits inside the published mean of 0.79 plus or minus an SD of 0.36. The dialysate flow for the CVVHDF arm is not restated in the CVVHDF paragraph of the Methods and is assumed to follow the same weight band; see the vignette Assumptions section. Meaningful only when RRT_CRRT_ACTIVE = 1; converted to L/h inside model().", 
#>             source_name = "QDial, QFil"))
#>     covariatesDataExcluded <- list(AGE = list(description = "Age", 
#>         units = "years", type = "continuous", notes = "Screened on the body clearance and not retained (Broeker 2018 Methods, Pharmacometric analysis: 'Age, sex, serum creatinine, creatinine clearance (Cockcroft-Gault), and bilirubin were tested as covariates on the body clearance'; only bilirubin survived the likelihood-ratio criterion). Table 1: median 69 years, range 37 to 81. The protocol excluded patients over 85 or under 18."), 
#>         SEXF = list(description = "Female sex indicator", units = "(binary)", 
#>             type = "binary", notes = "Screened on the body clearance and not retained (Methods, Pharmacometric analysis). The cohort is 10 male and 1 female (Table 1), so the covariate is effectively unidentifiable here regardless of any true effect."), 
#>         SCR = list(description = "Serum creatinine", units = "mg/dL", 
#>             type = "continuous", notes = "Screened on the body clearance and not retained (Methods, Pharmacometric analysis). Table 1: median 1.2 mg/dL, range 0.5 to 2.4. Tigecycline is not renally eliminated to a meaningful extent, and in an anuric CRRT-dependent cohort serum creatinine reflects the dialysis prescription rather than native renal function."), 
#>         CRCL = list(description = "Cockcroft-Gault creatinine clearance", 
#>             units = "mL/min", type = "continuous", notes = "Screened on the body clearance and not retained (Methods, Pharmacometric analysis, which names the Cockcroft-Gault estimator explicitly). Not tabulated per patient; by convention Cockcroft-Gault is not interpretable in a CRRT-dependent patient, which the null result is consistent with."), 
#>         WT = list(description = "Total body weight", units = "kg", 
#>             type = "continuous", notes = "Tested as an allometric scalar on the structural parameters with both a fixed and a freely estimated exponent and NOT retained -- Methods: 'Allometric scaling models using total body weight with fixed and freely estimated scaling parameters were evaluated'; Results: 'Allometric scaling with a fixed exponent did not improve the model significantly and was not included.' Table 1: median 80 kg, range 68 to 104. Body weight nonetheless remains indirectly load-bearing through the CRRT prescription, because the Methods set the dialysate flow by a 90 kg weight band -- that pathway is carried by RRT_CRRT_EFFLUENT_FLOW, not by any covariate effect on a structural parameter."))
#>     description <- "Two-compartment population PK model for intravenous tigecycline in 11 critically ill adults with acute kidney injury receiving continuous renal replacement therapy, 8 on continuous venovenous hemodialysis (CVVHD) and 3 on continuous venovenous hemodiafiltration (CVVHDF) (Broeker 2018). Plasma and CRRT effluent concentrations were fitted SIMULTANEOUSLY, which is what separates the dialysis clearance from the physiological body clearance. Total elimination from the central compartment is the ADDITIVE sum of a body arm (CLbody, 18.3 L/h typical, encoded as lcl, carrying a median-normalised power effect of total bilirubin with a NEGATIVE exponent so that cholestatic patients clear tigecycline more slowly) and a dialysis arm whose typical value depends on the CRRT modality (CLdial 1.69 L/h for CVVHD, 2.71 L/h for CVVHDF). The dialysis arm is gated by RRT_CRRT_ACTIVE and selected by RRT_CVVHDF_STATUS, with CVVHD as the reference modality; only the CVVHD arm carries interindividual variability, an IIV on the CVVHDF arm having collapsed to zero during estimation. Besides the plasma concentration Cc the model returns the flow-normalised effluent concentration Ceffluent with its own proportional residual error, obtained by rearranging the paper's first-principle dialysis equations. Age, sex, serum creatinine and Cockcroft-Gault creatinine clearance were screened on the body clearance and not retained, and allometric scaling by total body weight did not improve the fit."
#>     population <- list(species = "human", n_subjects = 11L, n_studies = 1L, 
#>         n_observations = 217L, age_range = "37 to 81 years", 
#>         age_median = "69 years", weight_range = "68 to 104 kg", 
#>         weight_median = "80 kg", sex_female_pct = 9.1, race_ethnicity = "Not reported.", 
#>         disease_state = "Critically ill adults in a 40-bed anaesthesiological ICU of a tertiary care hospital who required renal replacement therapy for acute kidney injury and were treated with tigecycline. Ten of the eleven were treated for complicated intra-abdominal infection and one for an Acinetobacter baumannii infection. Relevant co-conditions were liver failure or cirrhosis (four patients), liver transplantation (two), and extracorporeal membrane oxygenation (one). Two patients died during follow-up. APACHE II median 29, range 15 to 45 (Table 1). Exclusion criteria were age over 85 or under 18 years, severe liver insufficiency (Child-Pugh C), acute pancreatitis, concomitant anticoagulation therapy, and a history of tigecycline allergy. EudraCT 2012-005617-39.", 
#>         dose_range = "Standard tigecycline dosing: a 100 mg intravenous loading dose followed by 50 mg twice daily (Methods, Setting and study population). Sampling took place on day 4 of treatment, after at least 24 h on CRRT, and is therefore at steady state. Samples were drawn immediately before the start of the infusion (time 0), at 1 h (the end of infusion), and at 1.25, 1.5, 1.75, 2, 4, 6, 8 and 12 h, with effluent collected from the circuit effluent port at the same time points. The 1 h infusion duration is implied by the 'after 1 h (i.e., the end of infusion)' sampling description rather than stated as a prescription.", 
#>         regions = "Single centre, Germany (University Hospital Tuebingen / University Hospital Regensburg collaboration).", 
#>         renal_function = "All patients had acute kidney injury requiring continuous renal replacement therapy: 8 on CVVHD and 3 on CVVHDF, all using the Fresenius MultiFiltrate system with an Ultraflux AV 1000 S polysulfone membrane. CVVHD used Ci-Ca Dialysate K2 with 4% sodium citrate regional anticoagulation at a median citrate flow of 176 mL/h (under 3% of blood flow), targeting a post-filter ionised calcium of 0.25 to 0.35 mmol/L. CVVHDF used multiBic fluid for both dialysis and post-filter (post-dilution) replacement with an ultrafiltration rate of 1 L/h and unfractionated heparin anticoagulation. Blood and dialysate flows were set by weight band: 100 mL/min and 2000 mL/h under 90 kg, 120 mL/min and 2500 mL/h above. Serum creatinine median 1.2 mg/dL (0.5 to 2.4). No predilution correction of the CRRT clearance was applied, the citrate flow being very low relative to blood flow.", 
#>         hepatic_function = "Deliberately broad, which is what identifies the bilirubin covariate. Total bilirubin median 2.3 mg/dL with a range of 0.7 to 43.3 mg/dL; albumin 2.1 to 3.1 g/dL; total protein 3.5 to 6.4 g/dL (Table 1). Four patients had liver failure or cirrhosis and two were liver transplant recipients, while Child-Pugh C liver insufficiency was an exclusion criterion.", 
#>         notes = "Baseline demographics from Broeker 2018 Table 1, which lists all eleven patients individually. A total of 109 blood and 108 effluent samples were used, after excluding two 12-h blood samples whose very high concentrations indicated the draw followed the start of the next infusion. Total and free tigecycline were measured by validated HPLC-UV, with a plasma limit of quantification of 0.05 mg/L and intra- and interassay imprecision under 6%; the corresponding effluent values were 0.025 mg/L and under 9%. Free concentrations were determined by ultrafiltration at 1, 2 and 12 h, giving a median unbound fraction of 61% (range 45 to 94%) -- the model is fitted to TOTAL concentrations and carries no protein-binding term. NONMEM 7.4 with FOCEI executed via PsN 4.5.16, ADVAN1 and ADVAN3 routines; model selection by likelihood-ratio test (dOFV > 3.84), AIC for non-nested models, goodness-of-fit plots and visual predictive checks (n = 1000); parameter uncertainty from a nonparametric bootstrap (n = 1000). Shrinkage of the individual parameters was at most 26%.")
#>     reference <- "Broeker A, Wicha SG, Dorn C, Kratzer A, Schleibinger M, Kees F, Heininger A, Kees MG, Haeberle H. Tigecycline in critically ill patients on continuous renal replacement therapy: a population pharmacokinetic study. Crit Care. 2018;22(1):341. doi:10.1186/s13054-018-2278-4"
#>     units <- list(time = "h", dosing = "mg", concentration = "mg/L")
#>     vignette <- "Broeker_2018_tigecycline"
#>     ini({
#>         lcl <- 2.90690105984738
#>         label("Typical body clearance CLbody at the median bilirubin of 2.3 mg/dL (L/h)")
#>         lvc <- 4.07243972683405
#>         label("Typical central volume of distribution V1 (L)")
#>         lvp <- 5.03695260241363
#>         label("Typical peripheral volume of distribution V2 (L)")
#>         lq <- 4.03246915850401
#>         label("Typical distribution clearance Q (L/h)")
#>         lcl_hemodialysis_cvvhd <- 0.524728528934982
#>         label("Typical dialysis clearance CLdial on CVVHD (L/h)")
#>         lcl_hemodialysis_cvvhdf <- 0.99694863489161
#>         label("Typical dialysis clearance CLdial on CVVHDF (L/h)")
#>         e_tbili_cl <- -0.29
#>         label("Power exponent of median-normalised total bilirubin on body clearance (unitless)")
#>         propSd <- c(0, 0.169)
#>         label("Proportional residual error for pre-filter plasma concentrations (fraction)")
#>         propSd_Ceffluent <- c(0, 0.406)
#>         label("Proportional residual error for effluent concentrations (fraction)")
#>         etalcl ~ 0.190096
#>         etalvc ~ 1.229881
#>         etalq ~ 0.174724
#>         etalcl_hemodialysis_cvvhd ~ 0.189225
#>     })
#>     model({
#>         tbili_mgdL <- TBILI/17.1
#>         cl_body <- exp(lcl + etalcl) * (tbili_mgdL/2.3)^e_tbili_cl
#>         cl_hemodialysis <- RRT_CRRT_ACTIVE * ((1 - RRT_CVVHDF_STATUS) * 
#>             exp(lcl_hemodialysis_cvvhd + etalcl_hemodialysis_cvvhd) + 
#>             RRT_CVVHDF_STATUS * exp(lcl_hemodialysis_cvvhdf))
#>         cl <- cl_body + cl_hemodialysis
#>         vc <- exp(lvc + etalvc)
#>         vp <- exp(lvp)
#>         q <- exp(lq + etalq)
#>         kel <- cl/vc
#>         k12 <- q/vc
#>         k21 <- q/vp
#>         d/dt(central) <- -(kel + k12) * central + k21 * peripheral1
#>         d/dt(peripheral1) <- k12 * central - k21 * peripheral1
#>         effluent_flow <- max(RRT_CRRT_EFFLUENT_FLOW/1000, 0.001)
#>         Cc <- central/vc
#>         Ceffluent <- cl_hemodialysis * Cc/effluent_flow
#>         Cc ~ prop(propSd)
#>         Ceffluent ~ prop(propSd_Ceffluent)
#>     })
#> }

Population

Eleven critically ill adults in a 40-bed anaesthesiological ICU of a German tertiary care hospital, all requiring renal replacement therapy for acute kidney injury and treated with tigecycline: eight on CVVHD and three on CVVHDF. Ten were treated for complicated intra-abdominal infection and one for an Acinetobacter baumannii infection. Median age was 69 years (range 37 to 81), median weight 80 kg (68 to 104), and the cohort was 10 male and 1 female. APACHE II was a median of 29 (15 to 45); two patients died during follow-up.

The covariate that matters is hepatic, not renal. Four patients had liver failure or cirrhosis and two were liver transplant recipients, so total bilirubin spans two orders of magnitude – 0.7 to 43.3 mg/dL against a median of 2.3 mg/dL. That spread is what makes the bilirubin effect on clearance estimable in eleven subjects.

A total of 109 blood and 108 effluent samples were used, drawn on day 4 of treatment after at least 24 h on CRRT and therefore at steady state. Two 12-h blood samples were excluded as having been drawn after the start of the next infusion.

pop <- mod$population
tibble::tibble(
  Field = c(
    "Species", "Subjects", "Age", "Weight", "Female", "Disease", "Dosing"
  ),
  Value = c(
    pop$species,
    as.character(pop$n_subjects),
    paste0(pop$age_median, " (", pop$age_range, ")"),
    paste0(pop$weight_median, " (", pop$weight_range, ")"),
    paste0(pop$sex_female_pct, "%"),
    "Critically ill adults with AKI on CVVHD or CVVHDF",
    "100 mg IV loading dose, then 50 mg q12h (1 h infusions)"
  )
) |>
  knitr::kable(caption = "Study population (Broeker 2018 Table 1 and Methods).")

Study population (Broeker 2018 Table 1 and Methods).
Field	Value
Species	human
Subjects	11
Age	69 years (37 to 81 years)
Weight	80 kg (68 to 104 kg)
Female	9.1%
Disease	Critically ill adults with AKI on CVVHD or CVVHDF
Dosing	100 mg IV loading dose, then 50 mg q12h (1 h infusions)

Source trace

Every value in ini() and every non-trivial equation in model(), with its location in the source.

tibble::tribble(
  ~Quantity, ~Value, ~`Source location`,
  "Structural model", "2-compartment, first-order disposition, IV", "Results, Pharmacometric analysis (dOFV -113.77 vs 1-compartment)",
  "lcl (CLbody)", "18.3 L/h", "Table 2, theta_1 (RSE 11.0%; 95% CI 13.2, 22.7)",
  "e_tbili_cl", "-0.29", "Table 2, theta_2 (RSE 33.1%; 95% CI -0.68, -0.10)",
  "Bilirubin reference", "2.3 mg/dL", "Table 2 row heading; Results ('normalized by the population median of bilirubin, 2.3 mg/dL')",
  "lvc (V1)", "58.7 L", "Table 2 (RSE 21.3%; 95% CI 29.3, 101.6)",
  "lvp (V2)", "154 L", "Table 2 (RSE 9.5%; 95% CI 124.3, 196.8)",
  "lq (Q)", "56.4 L/h", "Table 2 (RSE 15.3%; 95% CI 41.1, 76.6)",
  "lcl_hemodialysis_cvvhd", "1.69 L/h", "Table 2, 'Dialysis clearance CVVHD' (RSE 15.4%; 95% CI 1.26, 2.27)",
  "lcl_hemodialysis_cvvhdf", "2.71 L/h", "Table 2, 'Dialysis clearance CVVHDF' (RSE 8.9%; 95% CI 2.31, 3.16)",
  "etalcl", "43.6 %CV", "Table 2, IIV column (post-covariate; 58.6% in the base model)",
  "etalvc", "110.9 %CV", "Table 2, IIV column",
  "etalq", "41.8 %CV", "Table 2, IIV column",
  "etalcl_hemodialysis_cvvhd", "43.5 %CV", "Table 2, IIV column",
  "IIV on V2 and CVVHDF arm", "none", "Table 2 prints a dash; Results ('IIV tended to zero during estimation')",
  "propSd", "16.9 %CV", "Table 2, 'sigma proportional, pre-filter plasma'",
  "propSd_Ceffluent", "40.6 %CV", "Table 2, 'sigma proportional, effluent'",
  "IIV model (Eq. 1)", "P_k,i = theta_k * exp(eta_k,i)", "Methods, Pharmacometric analysis",
  "Residual model (Eq. 2)", "proportional only (additive dropped)", "Methods Eq. 2; Results ('additive error tended to zero')",
  "Effluent observable (Eq. 3/4)", "Ceff = CLdial * CPla / Qeff", "Methods, rearranged from Eq. 3 (CVVHD) and Eq. 4 (CVVHDF)",
  "Qeff (CVVHD)", "QDial = 2000 mL/h (< 90 kg)", "Methods, Continuous renal replacement therapy",
  "Qeff (CVVHDF)", "QDial + QFil = 2000 + 1000 mL/h", "Methods ('Ultrafiltration rate (QFil) was 1 L/h')",
  "Modality per subject", "a = CVVHD, b = CVVHDF", "Table 1 footnotes",
  "Dosing", "100 mg IV load, then 50 mg q12h", "Methods, Setting and study population",
  "Infusion duration", "1 h", "Methods, Sampling ('after 1 h (i.e., the end of infusion)')"
) |>
  knitr::kable(caption = "Source trace for Broeker 2018 tigecycline.")

Source trace for Broeker 2018 tigecycline.
Quantity	Value	Source location
Structural model	2-compartment, first-order disposition, IV	Results, Pharmacometric analysis (dOFV -113.77 vs 1-compartment)
lcl (CLbody)	18.3 L/h	Table 2, theta_1 (RSE 11.0%; 95% CI 13.2, 22.7)
e_tbili_cl	-0.29	Table 2, theta_2 (RSE 33.1%; 95% CI -0.68, -0.10)
Bilirubin reference	2.3 mg/dL	Table 2 row heading; Results (‘normalized by the population median of bilirubin, 2.3 mg/dL’)
lvc (V1)	58.7 L	Table 2 (RSE 21.3%; 95% CI 29.3, 101.6)
lvp (V2)	154 L	Table 2 (RSE 9.5%; 95% CI 124.3, 196.8)
lq (Q)	56.4 L/h	Table 2 (RSE 15.3%; 95% CI 41.1, 76.6)
lcl_hemodialysis_cvvhd	1.69 L/h	Table 2, ‘Dialysis clearance CVVHD’ (RSE 15.4%; 95% CI 1.26, 2.27)
lcl_hemodialysis_cvvhdf	2.71 L/h	Table 2, ‘Dialysis clearance CVVHDF’ (RSE 8.9%; 95% CI 2.31, 3.16)
etalcl	43.6 %CV	Table 2, IIV column (post-covariate; 58.6% in the base model)
etalvc	110.9 %CV	Table 2, IIV column
etalq	41.8 %CV	Table 2, IIV column
etalcl_hemodialysis_cvvhd	43.5 %CV	Table 2, IIV column
IIV on V2 and CVVHDF arm	none	Table 2 prints a dash; Results (‘IIV tended to zero during estimation’)
propSd	16.9 %CV	Table 2, ‘sigma proportional, pre-filter plasma’
propSd_Ceffluent	40.6 %CV	Table 2, ‘sigma proportional, effluent’
IIV model (Eq. 1)	P_k,i = theta_k * exp(eta_k,i)	Methods, Pharmacometric analysis
Residual model (Eq. 2)	proportional only (additive dropped)	Methods Eq. 2; Results (‘additive error tended to zero’)
Effluent observable (Eq. 3/4)	Ceff = CLdial * CPla / Qeff	Methods, rearranged from Eq. 3 (CVVHD) and Eq. 4 (CVVHDF)
Qeff (CVVHD)	QDial = 2000 mL/h (< 90 kg)	Methods, Continuous renal replacement therapy
Qeff (CVVHDF)	QDial + QFil = 2000 + 1000 mL/h	Methods (‘Ultrafiltration rate (QFil) was 1 L/h’)
Modality per subject	a = CVVHD, b = CVVHDF	Table 1 footnotes
Dosing	100 mg IV load, then 50 mg q12h	Methods, Setting and study population
Infusion duration	1 h	Methods, Sampling (‘after 1 h (i.e., the end of infusion)’)

Virtual cohort

The natural cohort here is the published one: Broeker 2018 Table 1 lists all eleven patients individually, including the bilirubin value that drives the only retained covariate and the footnote marker that gives each patient’s CRRT modality. No distribution needs to be invented.

The effluent flow is reconstructed from the Methods prescription rather than tabulated per patient: dialysate flow is 2000 mL/h below 90 kg and 2500 mL/h above, and CVVHDF adds a fixed 1 L/h of ultrafiltrate.

patients <- tibble::tribble(
  ~subject, ~sex, ~age, ~ht, ~WT, ~apache, ~scr, ~bili_mgdL, ~modality,
  1L, "M", 69, 176,  69, 21, 1.7,  2.3, "CVVHD",
  2L, "F", 47, 160,  70, 21, 1.2,  9.2, "CVVHD",
  3L, "M", 81, 172,  68, 15, 1.0,  2.2, "CVVHD",
  4L, "M", 52, 180,  80, 45, 1.5, 24.0, "CVVHDF",
  5L, "M", 78, 178,  70, 25, 0.5,  3.5, "CVVHDF",
  6L, "M", 73, 172,  86, 29, 2.4,  1.8, "CVVHD",
  7L, "M", 56, 164, 104, 31, 1.2, 11.1, "CVVHDF",
  8L, "M", 37, 182,  85, 21, 1.3, 43.3, "CVVHD",
  9L, "M", 60, 180,  80, 35, 0.8,  1.8, "CVVHD",
  10L, "M", 74, 170, 73, 30, 0.7,  2.2, "CVVHD",
  11L, "M", 75, 180, 80, 30, 0.7,  0.7, "CVVHD"
) |>
  mutate(
    RRT_CVVHDF_STATUS = as.integer(modality == "CVVHDF"),
    RRT_CRRT_ACTIVE = 1L,
    # SI umol/L is the canonical TBILI unit; the paper reports mg/dL.
    TBILI = bili_mgdL * 17.1,
    # Methods: dialysate 2000 mL/h below 90 kg, 2500 above; CVVHDF adds QFil = 1 L/h.
    RRT_CRRT_EFFLUENT_FLOW = ifelse(WT < 90, 2000, 2500) + 1000 * RRT_CVVHDF_STATUS
  )

patients |>
  select(subject, modality, WT, bili_mgdL, RRT_CRRT_EFFLUENT_FLOW) |>
  rename(
    "Patient" = subject, "Modality" = modality, "Weight (kg)" = WT,
    "Bilirubin (mg/dL)" = bili_mgdL, "Effluent flow (mL/h)" = RRT_CRRT_EFFLUENT_FLOW
  ) |>
  knitr::kable(caption = "The eleven study patients (Broeker 2018 Table 1), with the CRRT effluent flow reconstructed from the Methods prescription.")

The eleven study patients (Broeker 2018 Table 1), with the CRRT effluent flow reconstructed from the Methods prescription.
Patient	Modality	Weight (kg)	Bilirubin (mg/dL)	Effluent flow (mL/h)
1	CVVHD	69	2.3	2000
2	CVVHD	70	9.2	2000
3	CVVHD	68	2.2	2000
4	CVVHDF	80	24.0	3000
5	CVVHDF	70	3.5	3000
6	CVVHD	86	1.8	2000
7	CVVHDF	104	11.1	3500
8	CVVHD	85	43.3	2000
9	CVVHD	80	1.8	2000
10	CVVHD	73	2.2	2000
11	CVVHD	80	0.7	2000


stopifnot(
  nrow(patients) == 11L,
  sum(patients$RRT_CVVHDF_STATUS) == 3L,
  sum(patients$RRT_CVVHDF_STATUS == 0L) == 8L,
  # Table 1 median row.
  median(patients$bili_mgdL) == 2.3,
  median(patients$WT) == 80
)

Simulation

Dosing follows the study protocol: a 100 mg loading dose, then 50 mg every 12 h, each as a 1-hour infusion. The model is simulated out to day 10 and observed over the 24-hour window from 240 to 264 h, which is unambiguously at steady state.

The observation grid is log-spaced after each dose. With an interindividual variability of 110.9 %CV on the central volume, some simulated subjects have a very small V1 and therefore a very sharp post-infusion distribution phase; a uniform grid under-resolves their peak and would understate AUC by several percent, breaking the CL x AUC = dose identity for reasons that have nothing to do with the model.

tau <- 12      # dosing interval (h)
t_ss <- 240    # start of the steady-state observation window (h)
n_rep <- 16L   # replicates per patient in the stochastic cohort

# Log-spaced observation times within each of the two dosing intervals, plus
# exact anchors at the window start, the mid-dose and the window end (PKNCA
# needs a record exactly at each interval boundary).
obs_times <- sort(unique(c(
  t_ss,
  t_ss + exp(seq(log(0.01), log(tau), length.out = 60)),
  t_ss + tau,
  t_ss + tau + exp(seq(log(0.01), log(tau), length.out = 60)),
  t_ss + 2 * tau
)))

make_events <- function(cohort) {
  dose <- cohort |>
    tidyr::expand_grid(tibble::tibble(
      time = c(0, tau), amt = c(100, 50), rate = c(100, 50),
      ii = c(0, tau), addl = c(0L, 39L)
    )) |>
    mutate(evid = 1L, cmt = "central", dvid = NA_integer_)
  obs <- cohort |>
    tidyr::expand_grid(tibble::tibble(time = obs_times)) |>
    mutate(
      amt = NA_real_, rate = NA_real_, ii = 0, addl = 0L,
      evid = 0L, cmt = NA_character_, dvid = 1L
    )
  bind_rows(dose, obs) |>
    arrange(id, time, desc(evid)) |>
    as.data.frame()
}

Observation rows carry dvid = 1L and cmt = NA_character_: the model has two endpoints (Cc and Ceffluent), so rxode2 maps each endpoint to its own compartment slot and a cmt = "central" observation row cannot be matched to an endpoint. useLinCmt = FALSE is passed to every solve because rxode2’s automatic ODE-to-linCmt() conversion corrupts the endpoint mapping for multi-output models.

# Typical-value arm: the eleven published covariate patterns, no IIV.
cohort_typ <- patients |> mutate(id = subject)
sim_typ <- rxode2::rxSolve(
  zmod, make_events(cohort_typ),
  keep = c("modality", "bili_mgdL", "WT"),
  useLinCmt = FALSE, returnType = "data.frame"
) |>
  filter(!is.na(Cc))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq', 'etalcl_hemodialysis_cvvhd'
#> Warning: multi-subject simulation without without 'omega'

stopifnot(nrow(sim_typ) > 0, all(sim_typ$Cc >= 0), !anyNA(sim_typ$Cc))

# Stochastic arm: each published patient replicated n_rep times with IIV drawn.
rxode2::rxSetSeed(20181220)
cohort_sim <- patients |>
  tidyr::expand_grid(tibble::tibble(rep = seq_len(n_rep))) |>
  mutate(id = row_number())

sim_pop <- rxode2::rxSolve(
  mod, make_events(cohort_sim),
  keep = c("modality", "bili_mgdL", "WT"),
  useLinCmt = FALSE, returnType = "data.frame"
) |>
  filter(!is.na(Cc))

# 176 subjects total: 128 on CVVHD and 48 on CVVHDF, both well under the
# 200-per-arm cohort cap.
stopifnot(
  dplyr::n_distinct(sim_pop$id) == 11L * n_rep,
  all(sim_pop$Cc >= 0),
  all(sim_pop$Ceffluent >= 0)
)

Structural checks

These are deterministic identities evaluated on the typical-value arm. Each one fails loudly if a specific transcription error was made, and none of them depends on which cohort the random-number stream happened to draw.

trapz <- function(t, y) sum(diff(t) * (head(y, -1) + tail(y, -1)) / 2)

ss <- sim_typ |>
  group_by(id, modality) |>
  summarise(
    cl_body = first(cl_body),
    cl_dial = first(cl_hemodialysis),
    cl_total = first(cl),
    auc24 = trapz(time, Cc),
    sc = mean(Ceffluent / Cc),          # saturation coefficient Ceff / CPla
    qeff = first(effluent_flow),
    .groups = "drop"
  ) |>
  mutate(
    # At steady state the amount eliminated over 24 h equals the 24-h dose.
    dose_recovered = cl_total * auc24,
    frac_dialysed = cl_dial / cl_total
  )

ss |>
  mutate(across(where(is.numeric), \(x) signif(x, 4))) |>
  select(id, modality, cl_body, cl_dial, cl_total, auc24, dose_recovered, sc, frac_dialysed) |>
  rename(
    "Patient" = id, "Modality" = modality, "CLbody (L/h)" = cl_body,
    "CLdial (L/h)" = cl_dial, "CL total (L/h)" = cl_total,
    "AUC24,ss (mg*h/L)" = auc24, "CL x AUC (mg)" = dose_recovered,
    "Ceff/CPla" = sc, "Fraction dialysed" = frac_dialysed
  ) |>
  knitr::kable(caption = "Typical-value structural checks, one row per published patient.")

Typical-value structural checks, one row per published patient.
Patient	Modality	CLbody (L/h)	CLdial (L/h)	CL total (L/h)	AUC24,ss (mg*h/L)	CL x AUC (mg)	Ceff/CPla	Fraction dialysed
1	CVVHD	18.300	1.69	19.990	5.005	100.0	0.8450	0.08454
2	CVVHD	12.240	1.69	13.930	7.180	100.0	0.8450	0.12130
3	CVVHD	18.540	1.69	20.230	4.946	100.0	0.8450	0.08355
4	CVVHDF	9.270	2.71	11.980	8.349	100.0	0.9033	0.22620
5	CVVHDF	16.200	2.71	18.910	5.290	100.0	0.9033	0.14330
6	CVVHD	19.650	1.69	21.340	4.689	100.0	0.8450	0.07920
7	CVVHDF	11.590	2.71	14.300	6.993	100.0	0.7743	0.18950
8	CVVHD	7.812	1.69	9.502	10.530	100.0	0.8450	0.17790
9	CVVHD	19.650	1.69	21.340	4.689	100.0	0.8450	0.07920
10	CVVHD	18.540	1.69	20.230	4.946	100.0	0.8450	0.08355
11	CVVHD	25.840	1.69	27.530	3.635	100.1	0.8450	0.06139

Mass balance. At steady state the amount cleared over a 24-hour window must equal the 24-hour dose of 100 mg (two 50 mg doses). This is the gate that catches a dialysis arm that has been silently dropped from the elimination path, a mis-transcribed clearance, or a dose-unit error.

stopifnot(
  # Pure numerical-integration error: the two sides use the same drawn
  # parameters, so this is tight by construction.
  all(abs(ss$dose_recovered - 100) < 0.1)
)

The bilirubin covariate. The Results quantify it directly: “Individual clearance values varied from 9.3 L/h (10th percentile) to 19.1 L/h (90th percentile) depending on the bilirubin concentration (24 mg/dL to 1.8 mg/dL).”

cl_at_bili <- function(b) {
  s <- rxode2::rxSolve(
    zmod,
    make_events(tibble::tibble(
      id = 1L, TBILI = b * 17.1, RRT_CRRT_ACTIVE = 0L,
      RRT_CVVHDF_STATUS = 0L, RRT_CRRT_EFFLUENT_FLOW = 2000
    )),
    useLinCmt = FALSE, returnType = "data.frame"
  )
  s$cl_body[1]
}
cov_check <- tibble::tibble(
  `Bilirubin (mg/dL)` = c(24, 2.3, 1.8),
  `Published CLbody (L/h)` = c(9.3, 18.3, 19.1),
  `Model CLbody (L/h)` = signif(vapply(c(24, 2.3, 1.8), cl_at_bili, numeric(1)), 4)
) |>
  mutate(`% diff` = signif(100 * (`Model CLbody (L/h)` / `Published CLbody (L/h)` - 1), 3))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq', 'etalcl_hemodialysis_cvvhd'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq', 'etalcl_hemodialysis_cvvhd'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq', 'etalcl_hemodialysis_cvvhd'
knitr::kable(cov_check, caption = "Bilirubin effect on body clearance against the values quoted in the Results.")

Bilirubin effect on body clearance against the values quoted in the Results.
Bilirubin (mg/dL)	Published CLbody (L/h)	Model CLbody (L/h)	% diff
24.0	9.3	9.27	-0.323
2.3	18.3	18.30	0.000
1.8	19.1	19.65	2.880


stopifnot(
  # The reference value and the 10th-percentile anchor are exact; the
  # 90th-percentile anchor is quoted against a rounded bilirubin of 1.8 mg/dL
  # (1.98 mg/dL reproduces 19.1 L/h exactly), hence the wider bound.
  abs(cov_check$`% diff`[cov_check$`Bilirubin (mg/dL)` == 2.3]) < 0.1,
  abs(cov_check$`% diff`[cov_check$`Bilirubin (mg/dL)` == 24]) < 1,
  abs(cov_check$`% diff`[cov_check$`Bilirubin (mg/dL)` == 1.8]) < 5
)

The saturation coefficient. The effluent observable is Ceff = CLdial * CPla / Qeff, so the simulated Ceffluent / Cc ratio must reproduce the paper’s published mean saturation coefficients. This is the gate that confirms the effluent equation was rearranged correctly and that the effluent-flow reconstruction from the Methods prescription is right – a wrong Qeff moves the ratio proportionally.

# (1) Structural identity: the effluent observable must equal CLdial / Qeff at
#     every time point, with Qeff correctly converted from mL/h to L/h. This is
#     the gate on the equation itself and is exact.
stopifnot(
  all(abs(ss$sc - ss$cl_dial / (ss$qeff)) < 1e-9),
  all(abs(ss$qeff - patients$RRT_CRRT_EFFLUENT_FLOW[order(patients$subject)] / 1000) < 1e-9)
)

# (2) Comparison against the published mean saturation coefficients, for a
#     patient in the sub-90 kg dialysate band that eight of the eleven
#     patients (and all but one of the three CVVHDF patients) fall into.
sat <- ss |>
  left_join(patients |> select(id = subject, WT), by = "id") |>
  filter(WT < 90) |>
  group_by(modality) |>
  summarise(model_sc = mean(sc), .groups = "drop") |>
  mutate(
    published_sc = ifelse(modality == "CVVHD", 0.79, 0.90),
    published_sd = ifelse(modality == "CVVHD", 0.36, NA_real_)
  )
sat |>
  mutate(across(where(is.numeric), \(x) signif(x, 4))) |>
  rename(
    "Modality" = modality, "Model Ceff/CPla" = model_sc,
    "Published mean" = published_sc, "Published SD" = published_sd
  ) |>
  knitr::kable(caption = "Model-implied versus published mean saturation coefficients (Broeker 2018 Results), for the sub-90 kg dialysate band.")

Model-implied versus published mean saturation coefficients (Broeker 2018 Results), for the sub-90 kg dialysate band.
Modality	Model Ceff/CPla	Published mean	Published SD
CVVHD	0.8450	0.79	0.36
CVVHDF	0.9033	0.90	NA


stopifnot(
  # CVVHDF: 2.71 / 3.000 = 0.903 against a published 0.90. This near-exact
  # agreement is what identifies QFil = 1 L/h as additive to a 2000 mL/h
  # dialysate flow, which the CVVHDF paragraph of the Methods does not restate.
  abs(sat$model_sc[sat$modality == "CVVHDF"] - 0.90) < 0.01,
  # CVVHD: 1.69 / 2.000 = 0.845 against a published mean of 0.79 (SD 0.36).
  # The published figure is a mean of per-patient measured Ceff/CPla ratios,
  # not a model parameter, so agreement within the published SD is the
  # meaningful criterion.
  abs(sat$model_sc[sat$modality == "CVVHD"] - 0.79) < 0.36
)

The one CVVHDF patient above the 90 kg dialysate threshold (patient 7, 104 kg) gets Qeff = 3500 mL/h under the Methods’ weight-band rule and therefore a model-implied saturation coefficient of 0.77. Including that patient would move the CVVHDF mean to 0.86, which does not round to the published 0.90 – so the published coefficient is evidence that the higher dialysate band was not in fact applied to this patient. Since the paper tabulates neither the per-patient prescription nor the per-patient coefficient, the vignette follows the Methods as written and compares within the band; see Assumptions and deviations.

The dialysis arm is load-bearing. Turning RRT_CRRT_ACTIVE off must raise exposure by exactly the clearance ratio. Without this check a dialysis arm that never reached the ODE would leave every other gate above green.

off <- rxode2::rxSolve(
  zmod, make_events(cohort_typ |> mutate(id = subject, RRT_CRRT_ACTIVE = 0L)),
  keep = "modality", useLinCmt = FALSE, returnType = "data.frame"
) |>
  filter(!is.na(Cc)) |>
  group_by(id, modality) |>
  summarise(auc_off = trapz(time, Cc), ceff_off = max(Ceffluent), .groups = "drop")
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalq', 'etalcl_hemodialysis_cvvhd'
#> Warning: multi-subject simulation without without 'omega'

crrt <- ss |>
  select(id, modality, auc_on = auc24, cl_body, cl_total) |>
  left_join(off, by = c("id", "modality")) |>
  mutate(ratio = auc_off / auc_on, expected = cl_total / cl_body)

stopifnot(
  # Exposure rises by exactly CL_total / CL_body when the circuit stops.
  all(abs(crrt$ratio - crrt$expected) < 1e-3),
  # And the effluent observable goes to zero rather than dividing by zero.
  all(crrt$ceff_off == 0),
  # The CRRT contribution is real but moderate. The ratio is 1 + CLdial/CLbody,
  # so it is largest in the cholestatic patients whose body clearance is most
  # suppressed: 1.07 at the lowest bilirubin (patient 11, 0.7 mg/dL, CVVHD) up
  # to 1.29 at the highest on CVVHDF (patient 4, 24 mg/dL). A dialysis arm that
  # never reached the ODE would give exactly 1.
  all(crrt$ratio > 1.05), all(crrt$ratio < 1.35)
)
knitr::kable(
  crrt |>
    mutate(across(where(is.numeric), \(x) signif(x, 4))) |>
    select(id, modality, auc_on, auc_off, ratio, expected) |>
    rename(
      "Patient" = id, "Modality" = modality, "AUC24 on CRRT" = auc_on,
      "AUC24 off CRRT" = auc_off, "Ratio" = ratio, "CLtotal/CLbody" = expected
    ),
  caption = "Effect of stopping CRRT on steady-state exposure."
)

Effect of stopping CRRT on steady-state exposure.
Patient	Modality	AUC24 on CRRT	AUC24 off CRRT	Ratio	CLtotal/CLbody
1	CVVHD	5.005	5.467	1.092	1.092
2	CVVHD	7.180	8.171	1.138	1.138
3	CVVHD	4.946	5.397	1.091	1.091
4	CVVHDF	8.349	10.790	1.292	1.292
5	CVVHDF	5.290	6.174	1.167	1.167
6	CVVHD	4.689	5.092	1.086	1.086
7	CVVHDF	6.993	8.628	1.234	1.234
8	CVVHD	10.530	12.800	1.216	1.216
9	CVVHD	4.689	5.092	1.086	1.086
10	CVVHD	4.946	5.397	1.091	1.091
11	CVVHD	3.635	3.872	1.065	1.065

Replicating Figure 1

Broeker 2018 Figure 1 is a visual predictive check over a dosing interval, with plasma concentrations on the left and effluent concentrations on the right: observed and predicted medians with the 10th-to-90th percentile band. The panels below reproduce that layout from the stochastic cohort over the steady-state dosing interval.

vpc <- sim_pop |>
  filter(time >= t_ss, time <= t_ss + tau) |>
  mutate(tad = time - t_ss) |>
  tidyr::pivot_longer(c(Cc, Ceffluent), names_to = "matrix", values_to = "conc") |>
  mutate(matrix = factor(
    matrix, levels = c("Cc", "Ceffluent"),
    labels = c("Plasma (pre-filter)", "Effluent")
  )) |>
  group_by(matrix, modality, tad) |>
  summarise(
    p10 = quantile(conc, 0.1), med = median(conc), p90 = quantile(conc, 0.9),
    .groups = "drop"
  )

ggplot(vpc, aes(tad, med, colour = modality, fill = modality)) +
  geom_ribbon(aes(ymin = p10, ymax = p90), alpha = 0.18, colour = NA) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~matrix, scales = "free_y") +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)", y = "Tigecycline concentration (mg/L)",
    colour = "Modality", fill = "Modality",
    title = "Replicates Figure 1 of Broeker 2018",
    subtitle = "Median (line) and 10th-90th percentile band over the steady-state dosing interval"
  ) +
  theme_bw()

The effluent panel sits below the plasma panel by the saturation coefficient – a fixed multiple of plasma at every time point, since the model derives effluent from plasma algebraically rather than from a separate compartment. The CVVHDF band is the higher of the two effluent curves, as its saturation coefficient is 0.90 against 0.79.

PKNCA validation

Non-compartmental analysis of the simulated steady-state profiles, over the 24-hour window the paper reports its exposures on.

sim_nca <- sim_pop |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, modality)

# The q12h schedule is expanded via addl, so the two doses that fall inside the
# 240-264 h window are materialised explicitly for PKNCA.
dose_nca <- cohort_sim |>
  select(id, modality) |>
  tidyr::expand_grid(tibble::tibble(time = c(t_ss, t_ss + tau), amt = 50))

conc_obj <- PKNCA::PKNCAconc(
  as.data.frame(sim_nca), Cc ~ time | modality + id,
  concu = "mg/L", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  as.data.frame(dose_nca), amt ~ time | modality + id, doseu = "mg"
)

intervals <- data.frame(
  start = t_ss, end = t_ss + 2 * tau,
  cmax = TRUE, cmin = TRUE, tmax = TRUE, auclast = TRUE, cav = TRUE
)

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

stopifnot(
  nrow(nca_tbl) > 0,
  # A gate that cannot go red is worse than no gate: confirm the parameters we
  # are about to compare were actually computed.
  all(c("auclast", "cmax", "cmin") %in% nca_tbl$PPTESTCD),
  !anyNA(nca_tbl$PPORRES[nca_tbl$PPTESTCD == "auclast"])
)

nca_tbl |>
  filter(PPTESTCD %in% c("auclast", "cmax", "cmin", "cav")) |>
  group_by(modality, PPTESTCD) |>
  summarise(
    median = median(PPORRES), p10 = quantile(PPORRES, 0.1),
    p90 = quantile(PPORRES, 0.9), .groups = "drop"
  ) |>
  mutate(
    Parameter = nlmixr2lib::ncaParamLabel(PPTESTCD),
    across(c(median, p10, p90), \(x) signif(x, 3))
  ) |>
  select(Parameter, modality, median, p10, p90) |>
  rename("Modality" = modality, "Median" = median, "10th pct" = p10, "90th pct" = p90) |>
  knitr::kable(caption = "Simulated steady-state NCA over the 24-hour window (mg/L; AUC in mg*h/L).")

Simulated steady-state NCA over the 24-hour window (mg/L; AUC in mg*h/L).
Parameter	Modality	Median	10th pct	90th pct
AUClast	CVVHD	4.750	2.9300	9.500
Cavg	CVVHD	0.198	0.1220	0.396
Cmax	CVVHD	0.550	0.3220	0.954
Cmin	CVVHD	0.109	0.0439	0.265
AUClast	CVVHDF	7.670	4.1600	11.900
Cavg	CVVHDF	0.319	0.1730	0.494
Cmax	CVVHDF	0.571	0.3360	1.170
Cmin	CVVHDF	0.209	0.0859	0.367

Comparison against the published exposures

Broeker 2018 does not publish an observed NCA table. What it does publish is the steady-state AUC24h distribution from its own clinical-trial simulation: a median of 6.15 mgh/L with a 10th-to-90th percentile range of 3.39 to 11.13 mgh/L in the CVVHD patients.

published <- tibble::tibble(modality = "CVVHD", auclast = 6.15)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_tbl |> filter(modality == "CVVHD"),
  reference = published,
  by = "modality",
  params = "auclast",
  units = c(auclast = "mg*h/L"),
  tolerance_pct = 20
)
knitr::kable(
  cmp,
  caption = "Simulated versus published steady-state AUC24h. * differs from the reference by more than 20%.",
  align = c("l", "l", "r", "r", "r")
)

Simulated versus published steady-state AUC24h. * differs from the reference by more than 20%.
NCA parameter	modality	Reference	Simulated	% diff
AUClast (mg*h/L)	CVVHD	6.15	4.75	-22.7%*

attr(cmp, "footnote")
#> [1] "* differs from reference by more than ±20%."

This row is a documented deviation, not a passing check, and it is discussed in the next section. The structural gates above – which are exact – all pass, so the disagreement is not a transcription error in ini().

# The paper's two simulation-derived statistics imply a body clearance, and they
# imply the SAME one -- which is not the 18.3 L/h of its own Table 2.
implied <- tibble::tibble(
  `Published statistic` = c(
    "Median steady-state AUC24h = 6.15 mg*h/L (CVVHD)",
    "Median fraction eliminated by CVVHD = 11.2%"
  ),
  `Implied CLtotal (L/h)` = signif(c(100 / 6.15, 1.69 / 0.112), 4),
  `Implied CLbody (L/h)` = signif(c(100 / 6.15 - 1.69, 1.69 / 0.112 - 1.69), 4)
) |>
  mutate(`Table 2 CLbody (L/h)` = 18.3)
knitr::kable(implied, caption = "The paper's two simulation outputs are mutually consistent but disagree with its own Table 2.")

The paper’s two simulation outputs are mutually consistent but disagree with its own Table 2.
Published statistic	Implied CLtotal (L/h)	Implied CLbody (L/h)	Table 2 CLbody (L/h)
Median steady-state AUC24h = 6.15 mg*h/L (CVVHD)	16.26	14.57	18.3
Median fraction eliminated by CVVHD = 11.2%	15.09	13.40	18.3


# Patient 1 sits exactly at the cohort median bilirubin of 2.3 mg/dL and is on
# CVVHD, so this is the typical CVVHD patient of Table 2.
model_frac <- ss |>
  filter(id == 1L) |>
  pull(frac_dialysed)
stopifnot(
  # The two published statistics agree with each other to within 10%...
  abs((100 / 6.15) / (1.69 / 0.112) - 1) < 0.1,
  # ...and the model's typical-value dialysed fraction is exactly the ratio of
  # the Table 2 clearances, i.e. 8.5%, not the published median of 11.2%.
  abs(model_frac - 1.69 / (18.3 + 1.69)) < 1e-9
)

Assumptions and deviations

Scale of the %CV interindividual-variability column. Broeker 2018 Table 2 reports IIV as a %CV without stating the conversion. The model reads it as omega = CV / 100, i.e. 100 times the log-scale standard deviation. Two things support that reading: the Table 2 caption states that RSE is “reported on standard deviation scale for variability parameters”, and the same column reports the proportional residual errors (16.9% and 40.6%), where a proportional sigma’s %CV is unambiguously 100 times its standard deviation – reading one column two ways would be inconsistent. The alternative lognormal convention omega^2 = log(1 + CV^2) would give variances of 0.174, 0.802, 0.161 and 0.173 instead of 0.190, 1.230, 0.175 and 0.189. The difference is about 5% in omega for the three moderate etas and larger for the central volume. No published statistic in the paper discriminates between the two readings; every gate in this vignette is either deterministic (evaluated with zeroRe(), so unaffected) or wide enough to admit either.

Dialysate flow on CVVHDF. The Methods give the weight-banded dialysate flow (2000 mL/h below 90 kg, 2500 above) in the CVVHD paragraph and state the CVVHDF ultrafiltration rate separately as 1 L/h, but do not restate the dialysate flow for CVVHDF. This vignette assumes the same weight band applies, giving an effluent flow of 3000 mL/h for a CVVHDF patient below 90 kg. That assumption is confirmed arithmetically by the paper’s own derived statistic rather than being merely plausible: 2.71 / 3.000 = 0.903 against a published CVVHDF saturation coefficient of 0.90.

The same arithmetic raises a question the paper cannot settle. Patient 7 is the only CVVHDF patient above the 90 kg threshold (104 kg), so the weight-band rule gives him Qeff = 3500 mL/h and a saturation coefficient of 0.77; the mean over all three CVVHDF patients would then be 0.86, which does not round to the published 0.90. All three at 3000 mL/h gives 0.903, which does. The published coefficient is therefore evidence that the higher dialysate band was not applied to this patient – but the paper tabulates neither the per-patient CRRT prescription nor the per-patient coefficient, so this cannot be confirmed. The cohort in this vignette follows the Methods as written, and the saturation- coefficient comparison is made within the sub-90 kg band so that the published value is compared against like. A user simulating a heavy CVVHDF patient should be aware that RRT_CRRT_EFFLUENT_FLOW is a free input and that the dialysis clearance itself (2.71 L/h) is not rescaled by it – the paper estimates the clearance directly, so the effluent flow affects only the effluent observable.

Effluent as a derived observable. The paper estimates the dialysis clearance directly and reports the saturation coefficients as descriptive Ceff/CPla ratios. The model therefore derives the effluent concentration from plasma algebraically, Ceffluent = cl_hemodialysis * Cc / Qeff, rearranged from Eq. 3 and Eq. 4, rather than carrying an effluent compartment. A consequence is that the effluent-to-plasma ratio is constant in time, so the effluent VPC panel is a scaled copy of the plasma panel. The paper notes one patient in whom effluent concentrations lagged, probably from adsorption to the membrane after a filter change; that time-dependence is not in the published model and is not reproduced here.

RRT_CRRT_ACTIVE is not estimated. Every subject was on CRRT throughout, so the column is identically 1 in the source analysis. It is carried so a downstream user can simulate circuit interruption, but setting it to 0 leaves the body clearance alone, which is the paper’s CLbody and not a validated off-CRRT clearance for this population.

Published steady-state AUC24h is not reproduced (deviation). The paper’s simulated median AUC24h of 6.15 mgh/L in CVVHD patients implies a total clearance of 16.3 L/h, whereas Table 2’s typical values give 18.3 + 1.69 = 19.99 L/h and a corresponding AUC24h of 5.00 mgh/L – about 23% lower than published. The disagreement is not a transcription error: the mass-balance identity CL x AUC24,ss = 100 mg holds to better than one part in a million in the table above, and the bilirubin covariate reproduces the Results’ own quoted clearances. It is also not isolated. The paper’s mass-balance analysis reports a median 11.2% of the dose eliminated by CVVHD, which implies a total clearance of 15.1 L/h, whereas the Table 2 typical values give 1.69 / 19.99 = 8.5%. The two simulation-derived statistics agree with each other – both imply a median body clearance near 13 to 15 L/h – and both disagree with Table 2 in the same direction, which points to the clinical-trial simulation having drawn bilirubin from a distribution whose central tendency is nearer the cohort mean of 7.9 mg/dL than its median of 2.3 mg/dL. At a bilirubin of 7.9 mg/dL the model gives CLbody = 12.8 L/h and AUC24h = 6.9 mg*h/L, which brackets the published figure. The simulation’s covariate-resampling scheme is not described in enough detail to reproduce, so the published exposure percentiles are recorded here as a deviation rather than used as a gate.

Not reproduced. The probability-of-target-attainment analysis of Figure 2 compares this model against two external published models (Van Wart 2013 for cSSSI/cIAI patients, Korth-Bradley 2011 for healthy volunteers) that are not part of this extraction, so Figure 2 is out of scope. The concentration-dependent unbound fraction characterised in Additional file 1 is descriptive: the model is fitted to total concentrations and carries no protein-binding term.

Screened but not retained. Age, sex, serum creatinine and Cockcroft-Gault creatinine clearance were tested on body clearance and not retained; allometric scaling by total body weight was tested with both fixed and estimated exponents and rejected. These are recorded in the model file’s covariatesDataExcluded rather than covariateData.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 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         dplyr_1.2.1           PKNCA_0.12.1         
#> [4] rxode2_5.1.8          nlmixr2lib_0.3.2.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.61           bslib_0.12.0       
#>  [4] rxode2lincmt_0.1.0  lattice_0.22-9      vctrs_0.7.3        
#>  [7] tools_4.6.1         generics_0.1.4      parallel_4.6.1     
#> [10] tibble_3.3.1        symengine_0.2.13    pkgconfig_2.0.3    
#> [13] data.table_1.18.6.1 checkmate_2.3.4     RColorBrewer_1.1-3 
#> [16] S7_0.2.2            desc_1.4.3          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       tidyr_1.3.2         openssl_2.4.2      
#> [34] cachem_1.1.0        nlme_3.1-169        tidyselect_1.2.1   
#> [37] digest_0.6.39       lotri_1.0.5         purrr_1.2.2        
#> [40] labeling_0.4.3      rxode2ll_2.0.18     fastmap_1.2.0      
#> [43] grid_4.6.1          cli_3.6.6           dparser_1.3.1-13   
#> [46] magrittr_2.0.5      withr_3.0.3         scales_1.4.0       
#> [49] backports_1.5.1     rmarkdown_2.32      otel_0.2.0         
#> [52] askpass_1.2.1       ragg_1.5.2          memoise_2.0.1      
#> [55] evaluate_1.0.5      knitr_1.52          rex_1.2.2          
#> [58] PreciseSums_0.7     rlang_1.3.0         downlit_0.4.5      
#> [61] Rcpp_1.1.2          glue_1.8.1          xml2_1.6.0         
#> [64] jsonlite_2.0.0      R6_2.6.1            systemfonts_1.3.2  
#> [67] fs_2.1.0