---
title: "Gentamicin (Abbasi 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gentamicin (Abbasi 2023)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r readModelDb("Abbasi_2023_gentamicin")()$reference`
* Description: `r readModelDb("Abbasi_2023_gentamicin")()$description`
* Article: <https://doi.org/10.2147/CPAA.S417298>

# Population

Abbasi et al. built a one-compartment intravenous gentamicin PK model by
systematically searching PubMed, EMBASE, SCOPUS, CINAHL, and EBSCO (MeSH terms
"gentamicin", "pharmacokinetics", "critically ill", "intensive care" up to
September 2022) and pooling the reported acute-phase PK parameters (Vd, CL, Ke)
and total body weight (TBW) from 21 previously published gentamicin PK studies in
critically ill adult ICU patients (Results; references 14-34 of the paper).
Among a total of 1215 pooled patients, 90.5% had confirmed severe infections.
Seven studies contributed acute-phase (first-dose, within 48-72 h) PK; two
contributed steady-state PK. The paper excluded pregnant women, patients on
extracorporeal membrane oxygenation, and patients on renal replacement therapy
(Methods, Pharmacokinetic Model Development).

The pooled acute-phase values (Results) were Vd = 0.33 +/- 0.20 L/kg, CL = 4.70
+/- 2.89 L/h, and Ke = 0.18 +/- 0.10 h^-1, with total body weight 70.8 +/- 19.9
kg. These pooled means and SDs were used as normal-distribution inputs to a
Monte Carlo Simulation (MCS) of 10,000 virtual patients in Crystal Ball, Oracle
(Methods, Monte Carlo Simulations). Total body weight was truncated at 40 kg to
infinity in the MCS. Doses of 5-10 mg/kg once daily as a 0.5-h IV infusion (24-h
dosing interval) were simulated over the first 72 h.

The pooled means and SDs are also accessible programmatically via the model's
`population` metadata:

```{r population, echo = TRUE}
mod_meta <- readModelDb("Abbasi_2023_gentamicin")()
mod_meta$population$n_subjects
mod_meta$population$n_studies
mod_meta$population$weight_range
mod_meta$population$disease_state
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Abbasi_2023_gentamicin.R`. This vignette
table collects them in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| Structural: 1-compartment IV, first-order elimination | n/a | Methods (Pharmacokinetic Model Development): "a one-compartment model with first-order elimination was developed". |
| Dose route: IV infusion, t_inf = 0.5 h, tau = 24 h | n/a | Methods (Pharmacokinetic Model Development): "t = infusion time (0.5 hour); T = dosing interval (24 hours)". |
| Dose = mg/kg * TBW | n/a | Methods (Pharmacokinetic Model Development): "The gentamicin dosage was calculated based on total body weight". |
| `lcl` (log CL) | log(4.70 L/h) | Results (paragraph 1): pooled acute-phase CL = 4.70 +/- 2.89 L/h across 21 studies. |
| `lvc` (log Vd per kg) | log(0.33 L/kg) | Results (paragraph 1): pooled acute-phase Vd = 0.33 +/- 0.20 L/kg across 21 studies. |
| `etalcl` (IIV on CL, log-scale variance) | log((2.89/4.70)^2 + 1) = 0.31989 | Results (paragraph 1): pooled CL SD 2.89 L/h; log-normal omega^2 = log(CV^2 + 1) conversion for nlmixr2 encoding. |
| `etalvc` (IIV on Vd, log-scale variance) | log((0.20/0.33)^2 + 1) = 0.31278 | Results (paragraph 1): pooled Vd SD 0.20 L/kg; log-normal omega^2 = log(CV^2 + 1) conversion. |
| `propSd` (proportional residual error) | fixed at 0 | Methods (no RUV; MCS is deterministic given sampled PK parameters). |
| TBW distribution | Mean 70.8 kg, SD 19.9 kg, truncated at 40 kg | Methods (Monte Carlo Simulations): "total body weight (kg) was set at the range 40 to infinity in MCS analysis"; Results: "Total body weight (70.8 +/- 19.9 kg) was used as dosing weight". |
| Dose range 5-10 mg/kg once daily | n/a | Methods (Monte Carlo Simulations): "various doses regimen ranging from 5-10 mg/kg once-daily gentamicin". |
| Efficacy targets AUC24h/MIC >= 110, Cmax/MIC ~ 8-10 | n/a | Methods (Probability of Target Attainment Prediction); Results Tables 1-2. |
| Nephrotoxicity targets AUC24h > 700 mg*h/L, Cmin > 2 mg/L | n/a | Methods (Probability of Target Attainment Prediction); Results Table 3. |

# Virtual cohort

A virtual cohort of 200 subjects per dose group (6 dose groups: 5, 6, 7, 8, 9,
10 mg/kg once daily) is generated with total body weight drawn from a truncated
normal distribution matching the paper's pooled TBW (mean 70.8 kg, SD 19.9 kg,
lower bound 40 kg to reproduce the MCS truncation described in Methods). Each
subject receives three consecutive daily 0.5-h IV infusions covering the 0-72 h
window in the paper.

```{r cohort}
set.seed(20230704L)

# Truncated normal (mean 70.8, SD 19.9, lower = 40 kg) via rejection sampling.
draw_tbw <- function(n, mean = 70.8, sd = 19.9, lower = 40) {
  out <- numeric(0)
  while (length(out) < n) {
    x <- rnorm(2L * n, mean = mean, sd = sd)
    out <- c(out, x[x >= lower])
  }
  out[seq_len(n)]
}

doses_mgkg <- c(5, 6, 7, 8, 9, 10)
n_per_arm  <- 200L

make_cohort <- function(dose_mgkg, n, id_offset) {
  tbw <- draw_tbw(n)
  amt_mg <- dose_mgkg * tbw

  # Three consecutive daily 0.5-h IV infusions to the central compartment
  # (rxode2 needs one row per event; rate = amt / 0.5 sets the 30-min infusion).
  dose_times <- c(0, 24, 48)
  dose_rows <- tibble(
    id       = rep(id_offset + seq_len(n), each = length(dose_times)),
    time     = rep(dose_times, times = n),
    amt      = rep(amt_mg, each = length(dose_times)),
    rate     = rep(amt_mg / 0.5, each = length(dose_times)),
    evid     = 1L,
    cmt      = "central",
    dose_mgkg = dose_mgkg,
    WT       = rep(tbw, each = length(dose_times))
  )

  # Observation grid: dense over the infusion (0-1 h), moderate elsewhere.
  # Repeat this pattern for each dosing interval so Cmax / Cmin are captured
  # on Day 1 and Day 3.
  obs_times_local <- c(seq(0, 1, by = 0.05), seq(1.5, 24, by = 0.5))
  obs_times <- unique(sort(c(
    obs_times_local,
    obs_times_local + 24,
    obs_times_local + 48,
    72
  )))

  obs_rows <- tibble(
    id       = rep(id_offset + seq_len(n), each = length(obs_times)),
    time     = rep(obs_times, times = n),
    amt      = NA_real_,
    rate     = NA_real_,
    evid     = 0L,
    cmt      = "central",
    dose_mgkg = dose_mgkg,
    WT       = rep(tbw, each = length(obs_times))
  )

  bind_rows(dose_rows, obs_rows) |>
    arrange(id, time, desc(evid))
}

events <- bind_rows(lapply(seq_along(doses_mgkg), function(i) {
  make_cohort(
    dose_mgkg = doses_mgkg[i],
    n         = n_per_arm,
    id_offset = (i - 1L) * n_per_arm
  )
})) |>
  mutate(treatment = paste0(dose_mgkg, " mg/kg/day"))

# Regression guard against duplicated (id, time, evid) rows across cohorts.
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Abbasi_2023_gentamicin")
sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("treatment", "dose_mgkg", "WT")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 1 (adapted) -- Concentration-time profiles over 72 h

Median and 5th-95th percentile envelope of simulated Cc over the 0-72 h window,
by once-daily dose. Compare qualitatively with the paper's Figure S1 (raw
time-concentration profiles) and Figure 1 (distributions of AUC24h/MIC ratios).

```{r figure-1-adapted}
sim |>
  group_by(treatment, time) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50, group = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), fill = "steelblue", alpha = 0.25) +
  geom_line(colour = "steelblue") +
  facet_wrap(~ treatment, ncol = 3) +
  labs(
    x = "Time (h)",
    y = "Gentamicin plasma concentration (mg/L)",
    title = "Simulated concentration-time profiles by once-daily dose",
    caption = "Median (line) and 5th-95th percentile envelope (band). n = 200 subjects per dose group. Compare to Figure S1 of Abbasi 2023."
  )
```

# PKNCA validation

Per-subject Cmax, Cmin, and AUC over each 24-h dosing interval are computed
with PKNCA. Cmax is evaluated over the first infusion (0-24 h) and reported as
Day 1. AUC24h and Cmin are evaluated on Day 1 (0-24 h) and Day 3 (48-72 h) to
match the paper's Table 1 and Table 3 windows.

```{r pknca-setup}
# Concentration frame -- keep column named Cc per nlmixr2lib convention.
# Filter only on !is.na(Cc) (do NOT drop time > 0 or Cc > 0; both would strip
# the anchor row that PKNCA needs).
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# Ensure a time = 0 row per (id, treatment); pre-dose Cc is 0 for IV.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

dose_df <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

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

intervals <- data.frame(
  start   = c(0, 48),
  end     = c(24, 72),
  cmax    = c(TRUE, FALSE),
  cmin    = TRUE,
  auclast = TRUE
)

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

## Comparison against Abbasi 2023 Table 1, 2, 3 -- Probability of Target Attainment

The paper reports probability of target attainment (PTA, % of MCS subjects who
meet a given target) for AUC24h/MIC >= 110 (Table 1, Day 1 and Day 3), Cmax/MIC
>= 8 and >= 10 for four MIC values (Table 2, Day 1 only), and nephrotoxicity
risk given AUC24h > 700 mg*h/L and Cmin > 2 mg/L (Table 3, Day 1 and Day 3). The
following code recomputes PTA on the 200-subjects-per-arm validation cohort and
compares it side-by-side with the published values.

```{r pta-compute}
# Wide per-subject NCA: one row per (treatment, id) with Cmax_d1, Cmin_d1,
# Cmin_d3, AUC24h_d1, AUC24h_d3.
per_subj <- nca_tbl |>
  dplyr::filter(!is.na(PPORRES)) |>
  dplyr::mutate(
    day = ifelse(start == 0, "d1", "d3")
  ) |>
  tidyr::pivot_wider(
    id_cols     = c(treatment, id),
    names_from  = c(PPTESTCD, day),
    values_from = PPORRES
  )

# PTA by dose x MIC x day (AUC24h/MIC >= 110)
compute_pta_auc <- function(df, mic) {
  df |>
    dplyr::group_by(treatment) |>
    dplyr::summarise(
      d1 = mean(auclast_d1 / mic >= 110) * 100,
      d3 = mean(auclast_d3 / mic >= 110) * 100,
      .groups = "drop"
    ) |>
    dplyr::mutate(mic = mic)
}

mic_values <- c(0.5, 1, 2, 4)
pta_auc_sim <- do.call(rbind, lapply(mic_values, function(m) compute_pta_auc(per_subj, m)))
```

### Table 1 replication -- PTA for AUC24h/MIC >= 110

```{r pta-table1-compare}
pta_auc_sim_wide <- pta_auc_sim |>
  tidyr::pivot_wider(
    id_cols     = treatment,
    names_from  = mic,
    values_from = c(d1, d3),
    names_glue  = "{.value}_MIC{mic}"
  ) |>
  dplyr::select(
    treatment,
    d1_MIC0.5, d1_MIC1, d1_MIC2, d1_MIC4,
    d3_MIC0.5, d3_MIC1, d3_MIC2, d3_MIC4
  )

# Paper Table 1 (Abbasi 2023): PTA for AUC24h/MIC >= 110 (columns are Day 1 for
# each MIC, then Day 3 for each MIC).
pta_auc_pub <- tibble::tribble(
  ~treatment,        ~d1_MIC0.5, ~d1_MIC1, ~d1_MIC2, ~d1_MIC4, ~d3_MIC0.5, ~d3_MIC1, ~d3_MIC2, ~d3_MIC4,
  "5 mg/kg/day",     77.3,       34.7,     6.3,      0.4,      80.2,       38.2,     7.4,      0.4,
  "6 mg/kg/day",     85.4,       46.5,     11.1,     0.9,      87.7,       50.0,     12.8,     1.0,
  "7 mg/kg/day",     90.9,       55.8,     16.4,     1.8,      92.6,       59.5,     18.5,     2.1,
  "8 mg/kg/day",     94.3,       64.8,     22.0,     3.0,      95.5,       68.3,     24.9,     3.7,
  "9 mg/kg/day",     96.8,       71.6,     28.7,     4.7,      97.9,       75.0,     31.8,     5.4,
  "10 mg/kg/day",    98.0,       77.7,     34.8,     6.3,      98.6,       80.7,     38.1,     7.6
)

# Merge side-by-side per dose x MIC x day.
tab1_long_sim <- pta_auc_sim_wide |>
  tidyr::pivot_longer(-treatment, names_to = "cell", values_to = "sim_PTA")
tab1_long_pub <- pta_auc_pub |>
  tidyr::pivot_longer(-treatment, names_to = "cell", values_to = "pub_PTA")

tab1_compare <- dplyr::inner_join(tab1_long_pub, tab1_long_sim,
                                  by = c("treatment", "cell")) |>
  dplyr::mutate(
    diff_pct = sim_PTA - pub_PTA,
    day  = ifelse(grepl("^d1_", cell), "Day 1", "Day 3"),
    MIC  = sub("^d[13]_MIC", "", cell)
  ) |>
  dplyr::select(treatment, day, MIC, pub_PTA, sim_PTA, diff_pct) |>
  dplyr::arrange(treatment, day, as.numeric(MIC))

knitr::kable(
  tab1_compare |>
    dplyr::rename(
      "Dose"           = treatment,
      "Day"            = day,
      "MIC (mg/L)"     = MIC,
      "Published PTA (%)" = pub_PTA,
      "Simulated PTA (%)" = sim_PTA,
      "Difference (percentage points)" = diff_pct
    ),
  digits = 1,
  caption = "Table 1 (Abbasi 2023): Probability of target attainment for AUC24h/MIC >= 110. Simulated values from a 200-subjects-per-dose validation cohort; the paper's MCS used 10,000 per dose."
)
```

### Table 2 replication -- PTA for Cmax/MIC >= 8 and >= 10 (Day 1)

```{r pta-table2-compare}
compute_pta_cmax <- function(df) {
  out <- lapply(mic_values, function(m) {
    df |>
      dplyr::group_by(treatment) |>
      dplyr::summarise(
        `>=8`  = mean(cmax_d1 / m >= 8) * 100,
        `>=10` = mean(cmax_d1 / m >= 10) * 100,
        .groups = "drop"
      ) |>
      dplyr::mutate(mic = m)
  })
  do.call(rbind, out)
}

pta_cmax_sim <- compute_pta_cmax(per_subj) |>
  tidyr::pivot_wider(
    id_cols     = treatment,
    names_from  = mic,
    values_from = c(`>=8`, `>=10`),
    names_glue  = "{.value}_MIC{mic}"
  )

pta_cmax_pub <- tibble::tribble(
  ~treatment,       ~`>=8_MIC0.5`, ~`>=10_MIC0.5`, ~`>=8_MIC1`, ~`>=10_MIC1`, ~`>=8_MIC2`, ~`>=10_MIC2`, ~`>=8_MIC4`, ~`>=10_MIC4`,
  "5 mg/kg/day",    95.2,          91.2,           76.4,        66.3,         41.2,        30.6,         13.2,        8.2,
  "6 mg/kg/day",    97.5,          95.0,           83.1,        75.6,         51.5,        39.9,         19.4,        12.5,
  "7 mg/kg/day",    98.7,          96.9,           88.0,        80.5,         59.3,        47.5,         25.5,        17.6,
  "8 mg/kg/day",    99.1,          97.8,           91.2,        85.5,         66.4,        55.3,         31.3,        21.5,
  "9 mg/kg/day",    99.6,          98.7,           93.4,        88.4,         72.2,        61.5,         37.6,        27.0,
  "10 mg/kg/day",   99.7,          99.1,           95.0,        91.2,         76.2,        66.9,         42.5,        31.8
)

tab2_long_sim <- pta_cmax_sim |>
  tidyr::pivot_longer(-treatment, names_to = "cell", values_to = "sim_PTA")
tab2_long_pub <- pta_cmax_pub |>
  tidyr::pivot_longer(-treatment, names_to = "cell", values_to = "pub_PTA")

tab2_compare <- dplyr::inner_join(tab2_long_pub, tab2_long_sim,
                                  by = c("treatment", "cell")) |>
  dplyr::mutate(
    diff_pct = sim_PTA - pub_PTA,
    target   = ifelse(grepl("^>=8_", cell), "Cmax/MIC >= 8", "Cmax/MIC >= 10"),
    MIC      = sub("^>=[0-9]+_MIC", "", cell)
  ) |>
  dplyr::select(treatment, target, MIC, pub_PTA, sim_PTA, diff_pct) |>
  dplyr::arrange(treatment, target, as.numeric(MIC))

knitr::kable(
  tab2_compare |>
    dplyr::rename(
      "Dose"           = treatment,
      "Target"         = target,
      "MIC (mg/L)"     = MIC,
      "Published PTA (%)" = pub_PTA,
      "Simulated PTA (%)" = sim_PTA,
      "Difference (percentage points)" = diff_pct
    ),
  digits = 1,
  caption = "Table 2 (Abbasi 2023): Probability of target attainment for Cmax/MIC >= 8 and >= 10 on Day 1 of therapy. Simulated values from a 200-subjects-per-dose validation cohort."
)
```

### Table 3 replication -- Nephrotoxicity risk (Day 1 and Day 3)

```{r pta-table3-compare}
neph_sim <- per_subj |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    auc700_d1 = mean(auclast_d1 > 700) * 100,
    auc700_d3 = mean(auclast_d3 > 700) * 100,
    cmin2_d1  = mean(cmin_d1  > 2)   * 100,
    cmin2_d3  = mean(cmin_d3  > 2)   * 100,
    .groups = "drop"
  )

neph_pub <- tibble::tribble(
  ~treatment,       ~auc700_d1, ~auc700_d3, ~cmin2_d1, ~cmin2_d3,
  "5 mg/kg/day",    0.10,       0.11,       3.9,       5.5,
  "6 mg/kg/day",    0.10,       0.13,       5.7,       8.0,
  "7 mg/kg/day",    0.21,       0.25,       8.0,       10.4,
  "8 mg/kg/day",    0.43,       0.46,       10.8,      13.2,
  "9 mg/kg/day",    0.72,       0.78,       13.0,      15.5,
  "10 mg/kg/day",   1.10,       1.29,       15.9,      18.4
)

tab3_compare <- dplyr::inner_join(
  neph_pub |> tidyr::pivot_longer(-treatment, names_to = "endpoint", values_to = "pub_pct"),
  neph_sim |> tidyr::pivot_longer(-treatment, names_to = "endpoint", values_to = "sim_pct"),
  by = c("treatment", "endpoint")
) |>
  dplyr::mutate(
    diff_pct = sim_pct - pub_pct,
    day      = ifelse(grepl("_d1$", endpoint), "Day 1", "Day 3"),
    metric   = ifelse(grepl("^auc700", endpoint),
                      "AUC24h > 700 mg*h/L (%)",
                      "Cmin > 2 mg/L (%)")
  ) |>
  dplyr::select(treatment, metric, day, pub_pct, sim_pct, diff_pct) |>
  dplyr::arrange(treatment, metric, day)

knitr::kable(
  tab3_compare |>
    dplyr::rename(
      "Dose"           = treatment,
      "Nephrotoxicity metric" = metric,
      "Day"            = day,
      "Published (%)"  = pub_pct,
      "Simulated (%)"  = sim_pct,
      "Difference (percentage points)" = diff_pct
    ),
  digits = 2,
  caption = "Table 3 (Abbasi 2023): Probability of developing nephrotoxicity by dose and duration of therapy. Simulated values from a 200-subjects-per-dose validation cohort."
)
```

# Assumptions and deviations

- **Cohort size**: 200 subjects per dose group (6 doses = 1200 subjects total),
  compared with the paper's 10,000 subjects per dose. The Monte Carlo
  discretisation error at n = 200 is on the order of +/- 3-4 percentage points
  for a 50% PTA and smaller at the tails (near 0% or 100%).
- **Log-normal random effects vs. the paper's normal distributions**: the
  paper's MCS sampled Vd/kg and CL from normal distributions with the reported
  means and SDs directly in Crystal Ball (truncation at zero was not stated).
  The nlmixr2lib canonical is log-normal random effects on structural PK
  parameters; this extraction transforms the reported %CV (SD/mean = 61.5% for
  CL, 60.6% for Vd) to log-scale variance via `omega^2 = log(CV^2 + 1)` and
  uses `lcl <- log(4.70)` / `lvc <- log(0.33)` so the log-median matches each
  reported mean. This is the standard nlmixr2lib encoding but is NOT identical
  to the paper's distributional assumption. At the reported CVs the log-normal
  has a positive skew, so the arithmetic mean of the simulated CL distribution
  (`4.70 * exp(0.32/2)` = 5.5 L/h) is ~17% higher than the paper's reported
  arithmetic mean of 4.70 L/h; the simulated distribution has fewer very-low-CL
  patients (in whom AUC24h is largest) and correspondingly lower PTA at
  intermediate MICs. The Table 1-3 comparison above shows this systematic
  offset: simulated AUC24h/MIC and Cmax/MIC PTAs run ~5-15 percentage points
  below the paper at MIC 1 and 2 mg/L, and simulated nephrotoxicity
  proportions run correspondingly lower. Users who want to reproduce the
  paper's PTA numbers exactly should re-sample Vd/kg and CL from the reported
  normal distributions externally rather than relying on the model's built-in
  log-normal IIV.
- **Redundant paper parameterisation**: Abbasi 2023 reports Vd, CL, and Ke as
  three independently-pooled parameters. Because Ke = CL / (Vd * TBW), only two
  of the three are structurally independent. This extraction uses CL and Vd as
  the primary parameters and derives Ke inside `model()`, giving an implied
  typical Ke of 4.70 / (0.33 * 70.8) = 0.201 h^-1 at the reference weight,
  slightly higher than the paper's independently pooled Ke of 0.18 h^-1. The
  ~12% inconsistency is a known artefact of the paper's independent parameter
  pooling and cannot be reconciled without departing from at least one of the
  reported values.
- **CL not scaled by body weight**: the paper reports pooled CL as total L/h
  (4.70 +/- 2.89 L/h) rather than L/h/kg. This extraction encodes CL as
  WT-independent to reproduce the paper's parameterisation exactly, even though
  a biologically motivated model would allometrically scale CL with weight.
  Downstream users who intend to simulate outside the pooled TBW range should
  be aware that CL will not scale with weight in this model.
- **Residual error set to 0**: Abbasi 2023 does not fit an OMEGA / SIGMA to
  observed concentrations; the MCS is deterministic given each virtual patient's
  sampled PK parameters and body weight. `propSd` is therefore `fixed(0)` to
  reproduce the paper's noise-free simulation.
- **TBW distribution**: the paper's MCS truncated TBW at 40 kg (Methods, Monte
  Carlo Simulations). This extraction reproduces that truncation via rejection
  sampling from N(70.8, 19.9^2) restricted to TBW >= 40 kg.
- **Table S2 (`Additional parameters used in the model`) not on disk**: the
  paper references a Table S2 in the Results section. All parameter values
  needed to reproduce the model appear in the main text (Results paragraph 1
  and Methods); the on-disk main PDF was used as the sole extraction source.
  Table S2 is inferred from context (Methods) to consolidate MCS input choices
  (dose range 5-10 mg/kg, MIC 0.5-4 mg/L, TBW >= 40 kg, dosing interval 24 h,
  infusion time 0.5 h) that are already stated in the main text.
- **Age / sex / race**: not reported at the pooled level in the main text.
  The 21 contributing studies (references 14-34) enrolled critically ill
  adult ICU patients; pregnant women, patients on ECMO, and patients on renal
  replacement therapy were excluded.

# Session info

```{r sessioninfo}
sessionInfo()
```
