---
title: "Ceftazidime-avibactam PK/PD in ICU and control patients (Zazo 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ceftazidime-avibactam PK/PD in ICU and control patients (Zazo 2024)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Cohort size per arm. 16 arms x N_PER_ARM subjects are simulated over the
# 168 h treatment week; the skill caps arms at 200 participants.
N_PER_ARM <- 100L
# Observation grid. T > MIC is read off this grid, so the step sets the
# resolution of that statistic (0.25 h -> ~0.15% of a week).
DT <- 0.25
set.seed(20240909)
```

# Model and source

Zazo, Aguazul and Lanao built a single semi-mechanistic PK/PD model of
ceftazidime-avibactam (CAZ-AVI) against *Pseudomonas aeruginosa* and applied it
to two distinct populations with their own pharmacokinetic parameter sets: a
non-ICU **control** population and a critically ill **ICU** population. Because
the two populations differ not only in typical values but also in their
variability (Table 5 reports a separate CV for every parameter in each
population), and nlmixr2 cannot express a covariate-dependent OMEGA, the
extraction carries them as two model files that share one vignette.

```{r models}
mod_control <- readModelDb("Zazo_2024_ceftazidime_avibactam_control")
mod_icu     <- readModelDb("Zazo_2024_ceftazidime_avibactam_icu")

# rxUi objects, needed to pipe model() edits in the static time-kill check.
ui_control <- rxode2::rxode(mod_control)
ui_icu     <- rxode2::rxode(mod_icu)
```

* Citation: `r ui_control$reference`
* Control model: `r ui_control$description`
* ICU model: `r ui_icu$description`
* Article: <https://doi.org/10.3390/antibiotics13090861>
* Supplement: <https://www.mdpi.com/article/10.3390/antibiotics13090861/s1>

# Population

This is a **simulation study**, not a fit to individual patient data. The
authors ran 1000 Monte Carlo replicates per subpopulation over the first week
(168 h) of treatment in GoldSim Pro v10.5. Two populations were simulated -
control (not admitted to the ICU) and ICU - each subdivided into six renal
strata by creatinine clearance (100, 60, 40, 20, 10 and 3 mL/min). Creatinine
clearance was drawn log-normally about each stratum mean with CV 20%, and body
weight log-normally with mean 75 kg and CV 13% (Section 4.1).

The pharmacokinetic parameters (Table 5) were derived by the authors from
published studies rather than estimated here: Stein 2019 (critically ill
CAZ-AVI), Merdjan 2017 (avibactam in renal impairment) and Welage 1984
(ceftazidime in renal insufficiency). The bacterial growth and kill parameters
(Table 6) were adapted from the Sy 2018 time-kill analysis of *P. aeruginosa*
strain 2154.

The MIC assumed throughout is the MIC50 for meropenem-non-susceptible
*P. aeruginosa*: 4 mg/L for ceftazidime and 1 mg/L for avibactam. The PK/PD
targets are T > MIC > 90% and Cmin/MIC > 1.3.

The same information is available programmatically via
`readModelDb("Zazo_2024_ceftazidime_avibactam_icu")()$population`.

# Source trace

| Equation / parameter | Value | Source location |
|---|---|---|
| `V = Vd * WT` | n/a | Eq 3, Section 4.2 |
| `Cl = Cli + CLs * CRCL` | n/a | Eq 4, Section 4.2 |
| `d/dt(bact_active)` (logistic growth, delayed kill, conversion loss) | n/a | Eq 5, Section 4.2 |
| `kdeath1` (sigmoidal Emax; AVI lowers the CAZ EC50 biexponentially) | n/a | Eq 6, Section 4.2 |
| `d/dt(bact_resting)` | n/a | Eq 7, Section 4.2 |
| `i1`, `i2` (delay functions) | n/a | Eq 8, Section 4.2 |
| `degrate` (bacterial degradation of CAZ, inhibited by AVI) | n/a | Eq 9, Section 4.2 |
| `lvd_caz` | 0.21 / 0.46 L/kg (control / ICU) | Table 5 |
| `lvd_avi` | 0.27 / 0.68 L/kg | Table 5 |
| `lcli_caz` | 0.39 / 1.15 L/h | Table 5 |
| `lcli_avi` | 0.53 / 0.89 L/h | Table 5 |
| `lcls_caz` | 0.06 / 0.04 (L/h)/(mL/min) | Table 5 |
| `lcls_avi` | 0.12 / 0.10 (L/h)/(mL/min) | Table 5 |
| `llog10nmax` | 9.89 log10 CFU/mL | Table 6 (`Nmax`) |
| `lkgrowth1` | 0.346 1/h | Table 6 |
| `kgrowth2` (derived) | `kgrowth1 / 1e7` | Table 6 |
| `lemax` | 0.240 1/h | Table 6 |
| `lec50_a`, `lec50_b` | 52.3, 12.6 mg/L | Table 6 (`A`, `B`) |
| `lalpha` | 2.38 L/mg | Table 6 |
| `lbeta` | 9.67e-2 L/mg | Table 6; **sign erratum**, see below |
| `lhill` | 2.60 | Table 6 (`gamma`) |
| `ldelta1` | 4.23e-2 1/h (fixed) | Table 6; **sign erratum** |
| `ldelta2` | 0.213 1/h | Table 6 |
| `lkconv` | 5.0e-3 1/h (fixed) | Table 6 (`k1-2`); **sign erratum** |
| `ldegmax` | 7.71e-2 1/h | Table 6; **sign erratum** |
| `lkm_deg` | 8.5 log10 CFU/mL (fixed) | Table 6 (`Km`) |
| `lhill_deg` | 1.46 | Table 6 (`phi`) |
| `lic50_deg` | 1.96 mg/L | Table 6 (`IC50`) |
| All IIV variances | `log(CV^2 + 1)` from the tabulated CV% | Tables 5 and 6 |
| Initial inocula | P1 = 10^6 CFU/mL, P2 = 1/10^7 CFU/mL | Section 4.2 |
| Dosage regimens | see Table 4 | Table 4 |
| MIC (CAZ 4, AVI 1 mg/L) | n/a | Section 4.3 |

# Validation 1: static time-kill (Table 1)

The strongest quantitative check available is Table 1, which reports the
model-predicted bacterial density for the three static time-kill conditions of
Sy 2018. "Static" means the drug concentrations are held constant, so the PK
and the Eq 9 degradation are bypassed. Here that is done by piping the two
concentration definitions to covariates - the packaged model is not modified.

```{r timekill}
# rxUi objects have REFERENCE semantics: piping model() edits mutates the
# object in place, so build the static variant from a fresh rxUi.
make_static <- function(modfun) {
  rxode2::rxode(modfun) |>
    rxode2::model(Cc_caz <- CAZ_STATIC) |>
    rxode2::model(Cc_avi <- AVI_STATIC) |>
    rxode2::zeroRe()
}
static_control <- make_static(mod_control)

# Observation rows: ODE-state cmt + dvid, as in the main cohort below.
static_obs <- function(ids, times) {
  tidyr::crossing(id = ids, time = times) |>
    dplyr::mutate(evid = 0L, amt = NA_real_, rate = NA_real_,
                  cmt = "central_caz", dvid = 1L)
}

tk_cov <- data.frame(
  id = 1:3, CAZ_STATIC = c(2, 4, 8), AVI_STATIC = 4, WT = 75, CRCL = 100
)
tk_sim <- rxode2::rxSolve(
  static_control,
  dplyr::left_join(static_obs(1:3, seq(0, 168, by = DT)), tk_cov, by = "id"),
  omega = NA, sigma = NA, useLinCmt = FALSE
) |> as.data.frame()

tk_tab <- tk_sim |>
  dplyr::group_by(id) |>
  dplyr::slice_max(time, n = 1) |>
  dplyr::ungroup() |>
  dplyr::transmute(
    Condition = c("CAZ 2 + AVI 4 mg/L", "CAZ 4 + AVI 4 mg/L", "CAZ 8 + AVI 4 mg/L"),
    Simulated = round(bactTotal, 2),
    `Published (Table 1)` = c(9.23, 8.78, 6.57),
    `Observed (Sy 2018)` = c(9.1, 8.0, 4.7)
  ) |>
  dplyr::mutate(Difference = round(Simulated - `Published (Table 1)`, 2))

knitr::kable(
  tk_tab, align = c("l", "r", "r", "r", "r"),
  caption = "Bacterial density at 168 h under static concentrations (log10 CFU/mL). Replicates Table 1 of Zazo 2024."
)
```

All three reproduce the published predictions to within 0.4 log10 CFU/mL, and
the two higher-exposure conditions to within 0.15. This is the check that pins
down the two readings discussed under *Assumptions and deviations*: the
negative exponents in Table 6, and the exclusion of the Eq 9 degradation term
from the static arm.

The drug-free growth control should saturate at the carrying capacity
`log10(Nmax) = 9.89`:

```{r growth-control}
gc_sim <- rxode2::rxSolve(
  static_control,
  static_obs(1L, seq(0, 168, by = DT)) |>
    dplyr::mutate(CAZ_STATIC = 0, AVI_STATIC = 0, WT = 75, CRCL = 100),
  omega = NA, sigma = NA, useLinCmt = FALSE
) |> as.data.frame()

cat(sprintf(
  "Drug-free control: t=0 %.2f, t=24 %.2f, t=168 %.2f log10 CFU/mL (Nmax = 9.89)\n",
  gc_sim$bactTotal[1],
  gc_sim$bactTotal[which.min(abs(gc_sim$time - 24))],
  dplyr::last(gc_sim$bactTotal)
))
```

# Validation 2: derived PK parameters (Table 2)

Eqs 3 and 4 evaluated at the mean body weight of 75 kg should reproduce the
volumes and clearances tabulated in Table 2.

```{r table2}
pk_coef <- tibble::tribble(
  ~Population, ~Drug,          ~vd,   ~cli,  ~cls,
  "Control",   "Ceftazidime",  0.21,  0.39,  0.06,
  "Control",   "Avibactam",    0.27,  0.53,  0.12,
  "ICU",       "Ceftazidime",  0.46,  1.15,  0.04,
  "ICU",       "Avibactam",    0.68,  0.89,  0.10
)

published_t2 <- tibble::tribble(
  ~Population, ~Drug,         ~CRCL, ~cl_pub, ~v_pub,
  "Control", "Ceftazidime", 100, 6.08, 15.7, "Control", "Ceftazidime", 60, 3.82, 15.7,
  "Control", "Ceftazidime",  40, 2.26, 15.7, "Control", "Ceftazidime", 20, 1.59, 15.7,
  "Control", "Ceftazidime",  10, 0.93, 15.7, "Control", "Ceftazidime",  3, 0.48, 15.7,
  "Control", "Avibactam",   100, 12.2, 20.2, "Control", "Avibactam",   60, 7.57, 20.2,
  "Control", "Avibactam",    40, 5.19, 20.2, "Control", "Avibactam",   20, 2.99, 20.2,
  "Control", "Avibactam",    10, 1.65, 20.2, "Control", "Avibactam",    3, 0.73, 20.2,
  "ICU", "Ceftazidime",     100, 5.44, 34.5, "ICU", "Ceftazidime",     60, 3.73, 34.5,
  "ICU", "Ceftazidime",      40, 2.86, 34.5, "ICU", "Ceftazidime",     20, 2.05, 34.5,
  "ICU", "Ceftazidime",      10, 1.56, 34.5, "ICU", "Ceftazidime",      3, 1.22, 34.5,
  "ICU", "Avibactam",       100, 10.9, 51.0, "ICU", "Avibactam",       60, 6.90, 51.0,
  "ICU", "Avibactam",        40, 4.86, 51.0, "ICU", "Avibactam",       20, 2.99, 51.0,
  "ICU", "Avibactam",        10, 1.85, 51.0, "ICU", "Avibactam",        3, 1.06, 51.0
)

t2 <- published_t2 |>
  dplyr::left_join(pk_coef, by = c("Population", "Drug")) |>
  dplyr::mutate(
    v_model  = vd * 75,
    cl_model = cli + cls * CRCL,
    `V % diff`  = round(100 * (v_model - v_pub) / v_pub, 1),
    `Cl % diff` = round(100 * (cl_model - cl_pub) / cl_pub, 1)
  ) |>
  dplyr::transmute(
    Population, Drug, `Clcr (mL/min)` = CRCL,
    `V model (L)` = round(v_model, 2), `V Table 2 (L)` = v_pub, `V % diff`,
    `Cl model (L/h)` = round(cl_model, 2), `Cl Table 2 (L/h)` = cl_pub, `Cl % diff`
  )

knitr::kable(t2, caption = "Eqs 3-4 at 75 kg vs the Table 2 Monte Carlo summaries.")
```

The volumes match exactly. The clearances agree to within about 10% at most
strata: Table 2 reports the **geometric mean of the simulated distribution**,
whereas `Cli + CLs * Clcr` is the typical value, and the geometric mean of a
sum of log-normal variates is not the sum of their geometric means.

# Virtual cohort

```{r cohort}
regimens <- tibble::tribble(
  ~source,     ~pop,      ~crcl, ~caz_g, ~avi_g, ~tau,
  "SmPC",      "Control", 100,   2.00,   0.50,    8,
  "SmPC",      "Control",  60,   2.00,   0.50,    8,
  "SmPC",      "Control",  40,   1.00,   0.25,    8,
  "SmPC",      "Control",  20,   0.75,   0.19,   12,
  "SmPC",      "Control",  10,   0.75,   0.19,   24,
  "SmPC",      "Control",   3,   0.75,   0.19,   48,
  "SmPC",      "ICU",     100,   2.00,   0.50,    8,
  "SmPC",      "ICU",      60,   2.00,   0.50,    8,
  "SmPC",      "ICU",      40,   1.00,   0.25,    8,
  "SmPC",      "ICU",      20,   0.75,   0.19,   12,
  "SmPC",      "ICU",      10,   0.75,   0.19,   24,
  "SmPC",      "ICU",       3,   0.75,   0.19,   48,
  "Suggested", "Control",  10,   0.75,   0.19,   12,
  "Suggested", "Control",   3,   0.75,   0.19,   24,
  "Suggested", "ICU",      10,   0.75,   0.19,   12,
  "Suggested", "ICU",       3,   0.75,   0.19,   24
) |>
  dplyr::mutate(
    arm = sprintf("%s | %s | Clcr %g | (%.2f/%.2f)/%gh",
                  source, pop, crcl, caz_g, avi_g, tau)
  )

# Log-normal draws matched to the paper's ARITHMETIC mean and CV (Section 4.1).
rlnorm_mean_cv <- function(n, mean, cv) {
  s2 <- log(cv^2 + 1)
  stats::rlnorm(n, meanlog = log(mean) - s2 / 2, sdlog = sqrt(s2))
}

# Observation rows sit on an ODE STATE (never on an algebraic observable) and
# carry dvid = 1L; this model declares three endpoints, so without dvid rxode2
# cannot map the observation to a prediction. See
# known-vignette-failure-patterns.md patterns 2 and 5b.
make_arm <- function(row, id_offset) {
  ids <- id_offset + seq_len(N_PER_ARM)
  cov <- tibble::tibble(
    id = ids,
    WT = rlnorm_mean_cv(N_PER_ARM, 75, 0.13),
    CRCL = rlnorm_mean_cv(N_PER_ARM, row$crcl, 0.20),
    arm = row$arm, source = row$source, pop = row$pop,
    crcl_stratum = row$crcl, tau = row$tau
  )
  dose_times <- seq(0, 168 - row$tau, by = row$tau)
  one_dose <- function(mg, cmt) {
    tidyr::crossing(id = ids, time = dose_times) |>
      dplyr::mutate(evid = 1L, amt = mg, rate = mg / 2,
                    cmt = cmt, dvid = NA_integer_)
  }
  ev <- dplyr::bind_rows(
    one_dose(row$caz_g * 1000, "central_caz"),
    one_dose(row$avi_g * 1000, "central_avi"),
    tidyr::crossing(id = ids, time = seq(0, 168, by = DT)) |>
      dplyr::mutate(evid = 0L, amt = NA_real_, rate = NA_real_,
                    cmt = "central_caz", dvid = 1L)
  ) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
  dplyr::left_join(ev, cov, by = "id")
}

events <- dplyr::bind_rows(lapply(
  seq_len(nrow(regimens)),
  function(i) make_arm(regimens[i, ], id_offset = (i - 1L) * N_PER_ARM)
))

# Disjoint-ID guard: duplicate ids across arms would silently merge subjects.
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
cat(sprintf("%d arms x %d subjects = %d subjects, %s event rows\n",
            nrow(regimens), N_PER_ARM, dplyr::n_distinct(events$id),
            format(nrow(events), big.mark = ",")))
```

# Simulation

The control and ICU arms are solved with their respective model files and
recombined.

```{r simulate}
solve_pop <- function(mod, ev) {
  if (nrow(ev) == 0) return(NULL)
  rxode2::rxSolve(
    mod, events = ev,
    keep = c("arm", "source", "pop", "crcl_stratum", "tau", "WT", "CRCL"),
    useLinCmt = FALSE
  ) |>
    as.data.frame() |>
    # rxSolve returns every model variable; keep only what the figures, the
    # Table 3 replication and the NCA need, so peak memory stays modest.
    dplyr::select(id, time, Cc_caz, Cc_avi, bactTotal,
                  arm, source, pop, crcl_stratum, tau, WT, CRCL)
}

# rxSolve returns observation rows only, so there is no evid column to filter.
sim <- dplyr::bind_rows(
  solve_pop(mod_control, dplyr::filter(events, pop == "Control")),
  solve_pop(mod_icu,     dplyr::filter(events, pop == "ICU"))
)

stopifnot(
  dplyr::n_distinct(sim$id) == nrow(regimens) * N_PER_ARM,
  !anyNA(sim$Cc_caz), !anyNA(sim$bactTotal)
)
```

# Replicate published figures

## Figures 1 and 3 - severe renal impairment, SmPC vs shortened interval

```{r figures-1-3, fig.height = 6.5}
fig_data <- sim |>
  dplyr::filter(crcl_stratum == 10) |>
  dplyr::select(time, pop, source, tau, Cc_caz, Cc_avi, bactTotal) |>
  tidyr::pivot_longer(c(Cc_caz, Cc_avi, bactTotal),
                      names_to = "panel", values_to = "value") |>
  dplyr::group_by(pop, source, tau, panel, time) |>
  dplyr::summarise(
    mean = mean(value),
    lo = quantile(value, 0.05), hi = quantile(value, 0.95),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    panel = factor(panel, c("Cc_caz", "Cc_avi", "bactTotal"),
                   c("Ceftazidime (mg/L)", "Avibactam (mg/L)",
                     "Bacterial density (log10 CFU/mL)")),
    regimen = sprintf("%s (q%gh)", source, tau)
  )

mic_lines <- data.frame(
  panel = factor(c("Ceftazidime (mg/L)", "Avibactam (mg/L)"),
                 levels = levels(fig_data$panel)),
  mic = c(4, 1)
)

ggplot(fig_data, aes(time, mean)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "darkgreen", alpha = 0.20) +
  geom_line(colour = "red") +
  geom_hline(data = mic_lines, aes(yintercept = mic),
             colour = "blue", linetype = "dashed") +
  facet_grid(panel ~ pop + regimen, scales = "free_y") +
  labs(
    x = "Time (h)", y = NULL,
    title = "Clcr = 10 mL/min: SmPC vs suggested regimen",
    caption = paste("Replicates Figures 1 (control) and 3 (ICU) of Zazo 2024.",
                    "Red: mean. Green band: 5th-95th percentile.",
                    "Blue dashed: MIC (4 mg/L CAZ, 1 mg/L AVI).")
  ) +
  theme_bw(base_size = 9) +
  theme(legend.position = "none")
```

As in the paper, the SmPC q24h regimen lets the bacterial density oscillate and
regrow between doses in both populations, while halving the dosing interval to
q12h holds it near the limit of detection.

## Figure 5 - probability of exceeding a bacterial density

```{r figure-5}
final_density <- sim |>
  dplyr::group_by(id, arm, pop, source, crcl_stratum, tau) |>
  dplyr::slice_max(time, n = 1) |>
  dplyr::ungroup()

fig5 <- final_density |>
  dplyr::filter(pop == "ICU", crcl_stratum %in% c(3, 10)) |>
  dplyr::mutate(regimen = sprintf("Clcr %g, q%gh (%s)", crcl_stratum, tau, source))

grid_d <- seq(0, 10, by = 0.1)
fig5_curves <- fig5 |>
  dplyr::group_by(regimen) |>
  dplyr::reframe(
    density = grid_d,
    prob = vapply(grid_d, function(d) mean(bactTotal > d), numeric(1))
  )

ggplot(fig5_curves, aes(density, prob, colour = regimen)) +
  geom_step() +
  labs(x = "Bacterial density at 168 h (log10 CFU/mL)",
       y = "Probability of exceeding", colour = NULL,
       title = "Probability of exceeding a bacterial density, ICU patients",
       caption = "Replicates Figure 5 of Zazo 2024.") +
  theme_bw(base_size = 9) +
  theme(legend.position = "bottom")
```

# Replicate Table 3 - PK/PD indices

`T > MIC` is the fraction of the treatment week with ceftazidime above 4 mg/L;
`Cmin/MIC` uses the trough of the final dosing interval; the bacterial density
change is expressed relative to the 6 log10 CFU/mL initial inoculum, as in the
Table 3 footnote.

```{r table3}
per_subject <- sim |>
  dplyr::group_by(id, arm, source, pop, crcl_stratum, tau) |>
  dplyr::summarise(
    t_gt_mic = 100 * mean(Cc_caz > 4),
    cmin_mic = min(Cc_caz[time >= max(time) - dplyr::first(tau)]) / 4,
    dens_chg = 100 * (dplyr::last(bactTotal[order(time)]) - 6) / 6,
    .groups = "drop"
  )

fmt <- function(x, lo, hi, d = 1) {
  sprintf("%.*f (%.*f/%.*f)", d, x, d, lo, d, hi)
}

published_t3 <- tibble::tribble(
  ~source, ~pop, ~crcl_stratum, ~dens_pub, ~tmic_pub, ~cmin_pub,
  "SmPC", "Control", 100, -63.7,  95.4, 0.81,  "SmPC", "Control", 60, -78.3, 100, 4.42,
  "SmPC", "Control",  40, -76.0, 100.0, 4.25,  "SmPC", "Control", 20, -69.5, 100, 3.40,
  "SmPC", "Control",  10, -51.2, 100.0, 2.04,  "SmPC", "Control",  3, -44.7, 100, 2.08,
  "SmPC", "ICU",     100, -77.5, 100.0, 4.47,  "SmPC", "ICU",      60, -84.8, 100, 9.76,
  "SmPC", "ICU",      40, -80.0, 100.0, 6.43,  "SmPC", "ICU",      20, -65.7, 100, 3.50,
  "SmPC", "ICU",      10, -30.2, 100.0, 1.18,  "SmPC", "ICU",       3,  13.0, 59.6, 0.35,
  "Suggested", "Control", 10, -83.7, 100.0, 10.10,
  "Suggested", "Control",  3, -80.8, 100.0,  9.37,
  "Suggested", "ICU",     10, -79.5, 100.0,  6.00,
  "Suggested", "ICU",      3, -56.5, 100.0,  2.63
)

t3 <- per_subject |>
  dplyr::group_by(source, pop, crcl_stratum, tau) |>
  dplyr::summarise(
    dens = fmt(mean(dens_chg), quantile(dens_chg, .05), quantile(dens_chg, .95)),
    dens_m = mean(dens_chg),
    tmic = fmt(mean(t_gt_mic), quantile(t_gt_mic, .05), quantile(t_gt_mic, .95)),
    tmic_m = mean(t_gt_mic),
    cmin = fmt(mean(cmin_mic), quantile(cmin_mic, .05), quantile(cmin_mic, .95), 2),
    cmin_m = mean(cmin_mic),
    .groups = "drop"
  ) |>
  dplyr::left_join(published_t3, by = c("source", "pop", "crcl_stratum")) |>
  dplyr::arrange(source, pop, dplyr::desc(crcl_stratum))

t3 |>
  dplyr::transmute(
    Regimen = source, Population = pop, `Clcr (mL/min)` = crcl_stratum,
    `Interval (h)` = tau,
    `Density change % (sim)` = dens, `Density change % (pub)` = round(dens_pub, 1),
    `T>MIC % (sim)` = tmic, `T>MIC % (pub)` = round(tmic_pub, 1),
    `Cmin/MIC (sim)` = cmin, `Cmin/MIC (pub)` = round(cmin_pub, 2)
  ) |>
  knitr::kable(
    caption = paste("Replicates Table 3 of Zazo 2024. Simulated entries are",
                    "mean (5th/95th percentile).")
  )
```

```{r table3-agreement}
cat(sprintf(
  "Qualitative target agreement (T>MIC > 90%% and Cmin/MIC > 1.3): %d of %d arms.\n",
  sum((t3$tmic_m > 90) == (t3$tmic_pub > 90) &
        (t3$cmin_m > 1.3) == (t3$cmin_pub > 1.3)),
  nrow(t3)
))
cat(sprintf(
  "Bacterial-density change, simulated vs published: Pearson r = %.3f over %d arms.\n",
  stats::cor(t3$dens_m, t3$dens_pub), nrow(t3)
))
```

The two purely pharmacokinetic indices reproduce well: `T > MIC` lands within a
few percentage points of Table 3 in every arm, and `Cmin/MIC` preserves both
the ordering and the pass/fail verdict against the 1.3 target. The paper's
central conclusion also survives - the SmPC regimen for Clcr = 3 mL/min in ICU
patients is the worst arm on every index, it is the only arm the paper scores
as failing `T > MIC > 90%` (simulated 68.6% vs published 59.6%), and halving
its dosing interval rescues it.

The **bacterial-density change is systematically understated in magnitude**
(simulated means roughly -35% to +35%, against -85% to +13% published), though
it remains well correlated with the published values. This is a property of the
published parameter set rather than of the encoding, and is discussed as
deviation 10 below. Note also that the simulation puts the *control* Clcr =
3 mL/min arm on the wrong side of zero (a net increase, where the paper reports
-44.7%).

# PKNCA validation

Steady-state NCA is computed over the final dosing interval of the Clcr =
100 mL/min arms, where Table 1 reports predicted steady-state peaks.

```{r pknca}
nca_arms <- sim |>
  dplyr::filter(crcl_stratum == 100) |>
  dplyr::mutate(regimen = paste(pop, "Clcr 100"))

start_ss <- 160
end_ss   <- 168

nca_for <- function(conc_col, concu) {
  sim_nca <- nca_arms |>
    dplyr::filter(!is.na(.data[[conc_col]])) |>
    dplyr::transmute(id, time, Cc = .data[[conc_col]], regimen)

  # Defensive time-zero row (see pknca-recipes.md); existing rows win.
  sim_nca <- dplyr::bind_rows(
    sim_nca,
    sim_nca |> dplyr::distinct(id, regimen) |> dplyr::mutate(time = 0, Cc = 0)
  ) |>
    dplyr::distinct(id, regimen, time, .keep_all = TRUE) |>
    dplyr::arrange(id, regimen, time)

  dose_df <- events |>
    dplyr::filter(evid != 0, crcl_stratum == 100,
                  cmt == if (conc_col == "Cc_caz") "central_caz" else "central_avi") |>
    dplyr::transmute(id, time, amt, regimen = paste(pop, "Clcr 100")) |>
    dplyr::distinct()

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

  intervals <- data.frame(
    start = start_ss, end = end_ss,
    cmax = TRUE, tmax = TRUE, cmin = TRUE, auclast = TRUE, cav = TRUE
  )
  PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
}

nca_caz <- nca_for("Cc_caz", "mg/L")
nca_avi <- nca_for("Cc_avi", "mg/L")
```

## Comparison against the published steady-state peaks

```{r nca-compare-caz}
published_caz <- tibble::tribble(
  ~regimen,           ~cmax,
  "Control Clcr 100",  95.9,
  "ICU Clcr 100",      67.8
)

cmp_caz <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_caz, reference = published_caz, by = "regimen",
  units = c(cmax = "mg/L"), tolerance_pct = 20
)
knitr::kable(
  cmp_caz, align = c("l", "l", "r", "r", "r"),
  caption = "Ceftazidime steady-state Cmax vs Table 1 of Zazo 2024. * differs by >20%."
)
```

```{r nca-compare-avi}
published_avi <- tibble::tribble(
  ~regimen,           ~cmax,
  "Control Clcr 100",  15.9
)

cmp_avi <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_avi, reference = published_avi, by = "regimen",
  units = c(cmax = "mg/L"), tolerance_pct = 20
)
knitr::kable(
  cmp_avi, align = c("l", "l", "r", "r", "r"),
  caption = "Avibactam steady-state Cmax vs Table 1 of Zazo 2024. * differs by >20%."
)
```

```{r nca-full}
as.data.frame(summary(nca_caz)) |>
  knitr::kable(caption = "Full steady-state NCA summary, ceftazidime (interval 160-168 h).")
```

# Assumptions and deviations

**1. Table 6 exponent signs (erratum).** The published Table 6 prints four
values in E-notation with a superscript but **no minus sign**: `beta` as
"9.67 E2", `delta1` as "4.23 E2", `k1-2` as "5.0 E3" and `Degmax` as
"7.71 E2". All four are encoded here with **negative** exponents
(9.67e-2, 4.23e-2, 5.0e-3, 7.71e-2). The evidence:

* *Arithmetic.* Simulating the three static time-kill conditions with the
  positive exponents gives a bacterial density of essentially 0 log10 CFU/mL
  for all three, against the 9.23 / 8.78 / 6.57 printed in Table 1 - a sum of
  squared errors of 205, versus 0.16 for the negative reading (reproduced in
  *Validation 1* above).
* *Dimensional analysis.* `k1-2` is added to `kdeath,1` inside the bracket of
  Eq 5, so it must carry units of 1/h; 5000/h would eradicate the inoculum
  within seconds. `Degmax` = 771/h would give ceftazidime a 3-second half-life.
* *Plausibility.* `delta1` = 423/h would make the "delay" function reach 1
  within milliseconds, i.e. no delay at all, contradicting its stated role.

**2. `Nmax` is a log10 density, not a linear count.** Table 6 gives
`Nmax = 9.89` with units "CFU/mL", but Eq 5 divides by `log10(Nmax)`. Taking
the units literally gives `log10(9.89) = 0.995` and a carrying capacity near
0.7 log10 CFU/mL, which cannot produce the 6.57-9.23 log10 CFU/mL of Table 1.
The value 9.89 is therefore `log10(Nmax)`, i.e. `Nmax = 10^9.89 CFU/mL`.
The same applies to `Km` = 8.5, which the text under Eq 9 explicitly calls
"the log10 transformed CFU number density".

**3. Eq 7 is encoded literally, on the linear scale - and the data cannot
settle it.** Eq 5 is written as `d(log10(P1))/dt` while Eq 7 is written as
`d(P2)/dt` with linear `(P2)` and `(P1)`. Table 6 pulls the other way,
describing `kgrowth,2` as "associated with the log10 of the active population
P2", which would imply a log10-scale Eq 7. The printed equation is followed
here, per the standing "trust the equation over the prose" rule, and because
the stated P2 initial condition - "1/10^7 CFU/mL" - is expressed as a linear
concentration.

Readers should know that **the published numbers do not settle the question**.
Both readings were simulated and compared against the paper's own results:

* *Table 1 (static time-kill)* does not discriminate them at all - total
  densities are nearly identical (sum of squared errors 0.160 linear vs 0.152
  log10 at 168 h; 0.156 vs 0.152 at 72 h).
* *Table 3 (bacterial-density change)* discriminates weakly and inconsistently.
  Over a six-arm probe the printed linear reading fits the published means
  better (SSE 13,274 vs 21,088) while the log10 reading better reproduces the
  published 5th percentiles, which reach -99.3% - near-complete eradication -
  and which the linear reading cannot reach at all, because its resting
  population accumulates to about 10^6 CFU/mL and floors the total.

The printed equation is therefore kept. The readings differ mainly in the
*composition* of the total: under the linear reading the resting population
becomes a substantial fraction of the total within hours, whereas under the
log10 reading it stays negligible all week. Any downstream use that depends on
the size of the resting reservoir - rather than on the total density the paper
reports - should treat this as an open structural ambiguity. The resulting
model carries `bact_active` on the log10 scale and `bact_resting` on the linear
scale; that asymmetry is the paper's, not the extraction's.

**4. The Table 1 static time-kill predictions exclude Eq 9.** Holding the
concentrations constant - the definition of a *static* time-kill - reproduces
Table 1 (SSE 0.16); letting ceftazidime degrade under Eq 9 lets all three
conditions regrow to > 9.7 log10 CFU/mL (SSE 11.1). Unlike the Eq 7 question
above, this one *is* decisive. Eq 9 is therefore retained in the packaged
model, where it acts on the dynamic patient-simulation concentrations, but is
bypassed in *Validation 1* above.

**5. Estimation error encoded as inter-individual variability.** Section 4.2
states that the Table 5 CVs on `Cli` and `CLs` are *estimation errors*
(parameter uncertainty) while the CV on `Vd` is inter-individual variability.
The paper's Monte Carlo drew all of them the same way - one log-normal draw per
virtual patient - so all are encoded as `eta` terms with
`omega^2 = log(CV^2 + 1)`. The same applies to the Table 6 CVs, which the
paper likewise sampled per virtual patient. A user who wants only structural
variability can zero the PD etas with `rxode2::zeroRe()` or by editing `ini()`.

**6. No residual-error model.** The source is a forward Monte Carlo simulation
in GoldSim, not an estimation run, so there is no `$SIGMA` to transcribe. The
three residual SDs are `fixed(0)` rather than invented.

**7. `CRCL` is raw, not BSA-normalised.** The Eq 4 slope is reported per
mL/min and the strata are defined by raw creatinine clearance. The paper does
not state which creatinine-clearance formula was used.

**8. The Zenodo dataset is unavailable.** The Data Availability Statement cites
DOI 10.5281/zenodo.13731674, which is unregistered (HTTP 404 as of the date of
this extraction). Every parameter and equation above comes from the article
itself.

**9. Quantities the paper does not fully specify.** Table 3 does not state the
averaging window for `T > MIC` or which dosing interval `Cmin` is taken from;
the whole 168 h week and the final dosing interval are used here. Body weight
and creatinine clearance are described as log-normal with a stated *mean* and
CV, which is the parameterisation used in the virtual cohort.

The Table 3 footnote likewise does not say whether the bacterial-density change
is a percentage of the **log10** density or of the linear CFU/mL count. It is
the log10 scale: computing it on the linear scale gives values in the thousands
of percent (up to +188,000% in the control Clcr = 3 arm, where a handful of
regrowing subjects dominate the mean) against the -85% to +13% range Table 3
prints, and the correlation with the published values falls from 0.88 to 0.63.

**10. Bacterial-density change is understated in magnitude.** Across the 16
arms the simulated density change correlates well with Table 3 (`r` reported
above), but the simulated means are consistently less negative than published
(roughly -35% to +35% simulated against -85% to +13% published). The cause is
visible in the published parameter set itself, and is not an encoding choice:
`Emax` (0.240 1/h, the maximum ceftazidime kill rate) is **smaller** than
`kgrowth,1` (0.346 1/h, the maximum growth rate). Once the population is driven
low, the logistic bracket `1 - log10(P1 + P2)/log10(Nmax)` approaches 1, so
growth outruns the maximum achievable kill and the population rebounds to the
density at which the two balance - about 2.9 log10 CFU/mL, i.e. a change of
about -52%. Eradication is only reachable for virtual patients whose `Emax` and
`kgrowth,1` draws happen to invert that inequality, which is the mechanism
behind the strongly bimodal 5th/95th percentiles Table 3 reports. Switching Eq
7 to the log10 reading (deviation 3) does not remove the gap; it shifts it. The
PK-driven indices `T > MIC` and `Cmin/MIC`, which do not depend on the
bacterial sub-model, agree with Table 3 closely.
