---
title: "Beta-lactam and beta-lactamase-inhibitor PBP binding in Pseudomonas aeruginosa (Lopez-Arguello 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Beta-lactam and beta-lactamase-inhibitor PBP binding in Pseudomonas aeruginosa (Lopez-Arguello 2023)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.width = 9, fig.height = 6,
  warning = FALSE, message = FALSE
)
```

```{r libraries}
library(nlmixr2lib)
library(rxode2)
library(PKNCA)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Overview

Lopez-Arguello and colleagues measured the time course of covalent
penicillin-binding protein (PBP) inactivation by 15 beta-lactams and
beta-lactamase inhibitors (BLIs) in *Pseudomonas aeruginosa* PAO1, in two
parallel assays:

* **Intact (whole) cells** -- drug must first cross the outer membrane (OM) to
  reach the periplasmic PBPs, so the observed binding reflects the *net* rate of
  influx and PBP access (influx minus efflux minus beta-lactamase hydrolysis).
* **Lysed cells** (isolated PBP-containing membranes) -- the OM barrier is
  absent and drug is present in vast excess, so the observed binding reflects
  intrinsic PBP affinity alone.

A quantitative and systems pharmacology (QSP) model was fit **simultaneously**
to both assays. The difference between them identifies the rate of net influx
and PBP access, and the model implements a periplasmic **mass balance**: every
drug molecule that binds a PBP is consumed, so the six PBPs compete for the same
finite periplasmic pool. Because PBP5/6 accounts for 71% of all PBP molecules
per cell, it acts as a **decoy target (sink)** that soaks up slowly penetrating
drugs before they can reach the essential PBPs 1a, 1b, 2 and 3.

This vignette reproduces the published model for all 15 drugs and validates it
against the paper's own non-compartmental summary of PBP occupancy (Tables 2 and
3).

# Models in this family

The paper reports one common model structure with drug-specific parameters
(Table 1 and Table S1), so the extraction is one `.R` file per drug and one
vignette for the paper.

`MIC` and `Conc` (the studied extracellular concentration, `2x MIC` for the
beta-lactams and a fixed 4 mg/L for the BLIs) are columns 3 and 4 of Table 1.
`Conc` is the value supplied to each model's `CONC_<CODE>_MGL` covariate below.

```{r model-list}
pbp_drugs <- tibble::tribble(
  ~Drug,           ~Code, ~Class,          ~MIC,  ~Conc,
  "imipenem",      "IPM", "Carbapenem",       1,      2,
  "doripenem",     "DOR", "Carbapenem",       1,      2,
  "meropenem",     "MEM", "Carbapenem",     0.5,      1,
  "ertapenem",     "ETP", "Carbapenem",       4,      8,
  "ceftazidime",   "CAZ", "Cephalosporin",    1,      2,
  "cefepime",      "FEP", "Cephalosporin",    1,      2,
  "cefoxitin",     "FOX", "Cephalosporin", 1024,   2048,
  "aztreonam",     "ATM", "Monobactam",       4,      8,
  "piperacillin",  "PIP", "Penicillin",       4,      8,
  "carbenicillin", "CAR", "Penicillin",      48,     96,
  "ticarcillin",   "TIC", "Penicillin",      24,     48,
  "avibactam",     "AVI", "DBO-BLI",         NA,      4,
  "relebactam",    "REL", "DBO-BLI",         NA,      4,
  "sulbactam",     "SUL", "BLI",             NA,      4,
  "tazobactam",    "TZB", "BLI",             NA,      4
) |>
  mutate(
    Model    = paste0("LopezArguello_2023_", Drug, "_qsp"),
    ConcCol  = paste0("CONC_", Code, "_MGL")
  )

knitr::kable(
  pbp_drugs |> select(Model, Drug, Code, Class, `MIC (mg/L)` = MIC, `Conc (mg/L)` = Conc),
  caption = "The 15 model files contributed by this paper, with the MIC and studied extracellular concentration from Table 1."
)
```

# Experimental system

The `population` metadata of every model file records the same in vitro system.

```{r population}
pop <- rxode2::rxode2(readModelDb("LopezArguello_2023_imipenem_qsp"))$meta$population
tibble::tibble(Field = names(pop), Value = unlist(lapply(pop, paste, collapse = "; "))) |>
  knitr::kable(caption = "Experimental system (identical across the 15 model files except for the drug-specific concentration and estimation dataset).")
```

Late exponential-phase PAO1 cultures (7.6 log10 CFU/mL) in cation-adjusted
Mueller-Hinton broth at 37 C were exposed to each beta-lactam at 2x its MIC
(BLIs: a fixed 4 mg/L). Membranes were harvested at 0, 15, 30 and 60 min, and
the **unbound** PBPs were then labelled with 25 uM Bocillin FL and quantified by
SDS-PAGE band intensity. A normalised band intensity of 1.0 means no binding and
0.0 means complete inactivation.

# Source trace

Every model equation and every `ini()` entry, with its location in the source.

```{r source-trace}
tibble::tribble(
  ~Component, ~Source,
  "Rate_Influx/access = Rate_Influx/access,scaled x C_drug", "Main text, Materials and Methods, Eq 1 (image `aac.01603-22-m001`); Fig. S8 line 67",
  "MM_i = K_i x N_peri / (Km + N_peri)", "Fig. S8 lines 14-19 (authoritative); main-text Eq 2 (`m002`) writes the equivalent Vmax/(Km+N_peri) form",
  "d(N_PBPi)/dt = -MM_i x N_PBPi", "Fig. S8 lines 21-26; main-text Eq 3 (`m003`)",
  "d(N_peri)/dt = Rate_Influx/access - sum_i MM_i x N_PBPi", "Fig. S8 line 28; main-text Eq 4 (`m004`)",
  "Intact assay: influx active, N_peri(0) = 0", "Fig. S8 lines 66-68",
  "Lysed assay: influx = 0, N_peri(0) = 1e7 molecules", "Fig. S8 lines 69-72; supplement Results",
  "Acylation constants rescaled by 1/1000 (Table S1 is in 1e-3 /min)", "Fig. S8 lines 75-80",
  "Initial conditions N_PBPi(0) = N_i x (1 - Noise_i) x Fini_i", "Fig. S8 lines 84-89",
  "Output Y_i = N_PBPi + N_i x Noise_i x Fini_i", "Fig. S8 lines 102-107",
  "Negative-state guard max(N_PBPi, 0)", "Fig. S8 lines 94-100",
  "Residual variance V_i = (SDin_i + SDsl_i x Y_i)^2 -> combined1()", "Fig. S8 lines 111-116",
  "rate_influx_scaled (per drug)", "Table 1, column 'Rate_Influx/access,scaled (SE%)'",
  "CONC_<CODE>_MGL covariate (per drug)", "Table 1, column 'Studied extracellular drug conc. in PBP binding assays'",
  "CELLS_INTACT covariate", "Fig. S8 line 66 `IF (INTACT.EQ.1)`; canonical name ratified by operator sidecar (request-001 q1, answer A)",
  "k2_1a ... k2_5 (per drug); `fixed()` where footnote a applies", "Table S1 (footnote a = fixed after sensitivity analysis; footnote b = relebactam PBP2 shared with avibactam)",
  "km_pbp = 1000 drug molecules (fixed)", "Main text, Materials and Methods: \"The Km was fixed to 1,000 drug molecules\"",
  "n_pbp1a_0 = 153, 1b = 118, 2 = 50, 3 = 79, 4 = 99, 5/6 = 1232", "Main text, Materials and Methods, Mass balance equations (1,731 total PBP molecules per cell borrowed from E. coli; split by P. aeruginosa relative band intensity, Table S3)",
  "noise_*, fini_*, addSd_*, propSd_*", "Table S2, per estimation dataset"
) |>
  knitr::kable(caption = "Source trace for the model structure and parameters.")
```

## Units and dimensional analysis

```{r units-table}
tibble::tribble(
  ~Symbol, ~Units, ~Note,
  "N_PBPi", "molecules/cell", "State: unbound molecules of PBP i per bacterial cell",
  "N_peri", "molecules/cell", "State: free drug molecules in the periplasm per cell",
  "K_i (kacyl)", "1/min", "Table S1 value x 1e-3",
  "Km", "molecules/cell", "Fixed at 1000",
  "MM_i", "1/min", "K_i x (dimensionless N_peri/(Km+N_peri))",
  "Rate_Influx/access,scaled", "molecules/min per mg/L", "Table 1",
  "C_drug", "mg/L", "Static extracellular concentration",
  "Rate_Influx/access", "molecules/min", "scaled rate x C_drug",
  "d(N_PBPi)/dt", "molecules/cell/min", "MM_i [1/min] x N_PBPi [molecules/cell]",
  "d(N_peri)/dt", "molecules/cell/min", "Influx [molecules/min] minus sum of binding fluxes [molecules/min]"
) |>
  knitr::kable(caption = "Units of every symbol in the ODE system. Each ODE right-hand side resolves to molecules per cell per minute.")
```

Both ODE lines balance: `MM_i` carries `1/min` (the Michaelis-Menten factor
`N_peri/(Km + N_peri)` is dimensionless), so `MM_i * N_PBPi` is
`molecules/cell/min`, matching `d(N_PBPi)/dt`; the same product appears with
opposite sign in the periplasmic balance alongside the influx rate.

# Simulation

Each model is solved twice: once with `intact = 1` (whole cells, drug enters at
`Rate_Influx/access`) and once with `intact = 0` (lysed cells, `N_peri(0)` =
10^7 molecules). These are deterministic typical-value solves -- the paper
reports population means only, no between-replicate variance components (see
Assumptions below).

```{r simulate}
obs_times <- c(0, 15, 30, 60)                 # the paper's sampling schedule
dense_times <- seq(0, 60, by = 0.5)           # for smooth trajectory plots

solve_one <- function(i, times, intact_value) {
  mod <- rxode2::rxode2(readModelDb(pbp_drugs$Model[i]))
  # The two covariates ride on the event table: CELLS_INTACT selects the assay
  # arm, and the drug-specific CONC_<CODE>_MGL column supplies the studied
  # extracellular concentration.
  # A multi-endpoint model needs an endpoint on each observation record; PBP1a
  # is used here and rxSolve returns all six endpoints as columns.
  ev <- data.frame(
    id = 1L, time = times, amt = 0, evid = 0L, cmt = "PBP1a",
    CELLS_INTACT = intact_value
  )
  ev[[pbp_drugs$ConcCol[i]]] <- pbp_drugs$Conc[i]
  as.data.frame(rxode2::rxSolve(mod, events = ev))
}

simulate_all <- function(times) {
  per_drug <- lapply(seq_len(nrow(pbp_drugs)), function(i) {
    lapply(c(Intact = 1, Lysed = 0), function(sc) {
      solve_one(i, times, sc) |>
        select(time, PBP1a, PBP1b, PBP2, PBP3, PBP4, PBP56) |>
        mutate(Drug = pbp_drugs$Drug[i], Code = pbp_drugs$Code[i], Cells = ifelse(sc == 1, "Intact", "Lysed"))
    }) |> bind_rows()
  })
  bind_rows(per_drug) |>
    pivot_longer(
      c(PBP1a, PBP1b, PBP2, PBP3, PBP4, PBP56),
      names_to = "PBP", values_to = "n_unbound"
    ) |>
    mutate(PBP = recode(PBP, PBP56 = "PBP5/6")) |>
    group_by(Drug, Code, Cells, PBP) |>
    # The paper normalises every band to its own 0 h intensity (1.0 = no binding).
    mutate(frac_unbound = n_unbound / n_unbound[time == 0]) |>
    ungroup()
}

sim_obs <- simulate_all(obs_times)
sim_dense <- simulate_all(dense_times)

drug_levels <- pbp_drugs$Drug
sim_obs <- mutate(sim_obs, Drug = factor(Drug, levels = drug_levels))
sim_dense <- mutate(sim_dense, Drug = factor(Drug, levels = drug_levels))

nrow(sim_dense)
```

# Published extent of PBP binding

Tables 2 (imipenem through aztreonam) and 3 (piperacillin through tazobactam)
report the average fraction of unbound PBPs over the first 60 min -- the linear
trapezoidal AUC divided by 60. A value of 1.00 means no binding and 0.00 means
complete inactivation. These 180 values are the paper's quantitative summary of
the whole experiment and are used below both to annotate Figure 4 and as the
reference for the NCA comparison.

```{r published-values}
published_wide <- tibble::tribble(
  ~PBP,     ~Cells,   ~IPM, ~DOR, ~MEM, ~ETP, ~CAZ, ~FEP, ~FOX, ~ATM, ~PIP, ~CAR, ~TIC, ~AVI, ~REL, ~SUL, ~TZB,
  "PBP1a",  "Intact", 0.29, 0.67, 0.68, 0.87, 0.61, 0.61, 0.24, 1.02, 0.79, 0.41, 0.48, 0.93, 0.87, 0.88, 0.98,
  "PBP1b",  "Intact", 0.30, 0.33, 0.37, 0.79, 0.87, 0.94, 0.38, 0.94, 0.89, 0.36, 0.43, 0.78, 0.84, 0.75, 1.05,
  "PBP2",   "Intact", 0.38, 0.30, 0.31, 0.63, 0.85, 0.96, 0.45, 1.02, 0.91, 0.82, 0.77, 0.80, 0.81, 0.97, 1.04,
  "PBP3",   "Intact", 0.35, 0.38, 0.44, 0.47, 0.51, 0.60, 0.47, 0.39, 0.48, 0.39, 0.48, 0.89, 0.81, 0.98, 1.10,
  "PBP4",   "Intact", 0.44, 0.22, 0.38, 0.33, 0.80, 0.69, 0.27, 1.01, 0.69, 0.49, 0.59, 0.75, 0.84, 0.74, 1.06,
  "PBP5/6", "Intact", 0.27, 0.36, 0.56, 0.93, 0.83, 0.98, 0.18, 0.97, 1.11, 0.80, 1.08, 0.90, 0.90, 0.84, 1.14,
  "PBP1a",  "Lysed",  0.19, 0.26, 0.29, 0.31, 0.17, 0.15, 0.14, 0.26, 0.14, 0.13, 0.10, 1.07, 0.95, 0.68, 0.79,
  "PBP1b",  "Lysed",  0.18, 0.26, 0.21, 0.27, 0.36, 0.34, 0.16, 0.24, 0.14, 0.15, 0.14, 0.64, 0.98, 0.48, 0.73,
  "PBP2",   "Lysed",  0.37, 0.34, 0.37, 0.35, 0.57, 0.47, 0.37, 0.64, 0.18, 0.20, 0.32, 0.79, 0.88, 0.79, 1.10,
  "PBP3",   "Lysed",  0.52, 0.47, 0.57, 0.35, 0.21, 0.14, 0.43, 0.23, 0.23, 0.17, 0.17, 1.07, 1.05, 1.09, 1.05,
  "PBP4",   "Lysed",  0.21, 0.23, 0.24, 0.20, 0.35, 0.17, 0.19, 0.36, 0.14, 0.17, 0.15, 0.61, 0.82, 0.70, 0.97,
  "PBP5/6", "Lysed",  0.19, 0.43, 0.54, 0.40, 0.82, 1.02, 0.13, 1.01, 0.82, 0.58, 0.67, 0.49, 0.86, 0.46, 0.90
)

published <- published_wide |>
  pivot_longer(-c(PBP, Cells), names_to = "Code", values_to = "auclast") |>
  left_join(select(pbp_drugs, Drug, Code), by = "Code") |>
  select(Drug, Cells, PBP, auclast) |>
  as.data.frame()

nrow(published)
```

# Replicate published figures

## Figure 3 / Figure S7 -- time course of PBP binding

Replicates Figure 3 (and its full-y-axis companion Figure S7) of
Lopez-Arguello 2023: model-predicted fraction of unbound PBPs in intact (solid)
and lysed (dashed) cells over 60 min.

```{r fig3, fig.height = 8, fig.width = 10}
ggplot(sim_dense, aes(time, frac_unbound, colour = PBP, linetype = Cells)) +
  geom_line(linewidth = 0.6) +
  facet_wrap(~Drug, ncol = 5) +
  scale_linetype_manual(values = c(Intact = "solid", Lysed = "22")) +
  coord_cartesian(ylim = c(0, 1.05)) +
  labs(
    x = "Time (min)", y = "Fraction of unbound PBPs",
    title = "Model-predicted PBP binding in intact and lysed P. aeruginosa PAO1"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

Imipenem and cefoxitin inactivate every PBP almost as fast in intact cells as in
lysed cells; ertapenem, aztreonam and piperacillin show a large intact-vs-lysed
gap, which is the signature of slow target-site penetration.

## Figure 2 -- difference between intact and lysed cells

Replicates Figure 2: `delta PBP` = (fraction unbound, intact) - (fraction
unbound, lysed) at each sampling time. A positive value means *less* binding in
intact cells, i.e. a penetration penalty.

```{r fig2, fig.height = 5, fig.width = 10}
delta_pbp <- sim_obs |>
  filter(time > 0) |>
  select(Drug, Code, PBP, time, Cells, frac_unbound) |>
  pivot_wider(names_from = Cells, values_from = frac_unbound) |>
  mutate(delta = Intact - Lysed)

ggplot(delta_pbp, aes(factor(time), delta, fill = PBP)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  facet_wrap(~Drug, ncol = 5) +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  labs(
    x = "Time (min)", y = "delta PBP (intact - lysed fraction unbound)",
    title = "Penetration penalty: model-predicted difference between assays"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

## Figure 4 -- PBP5/6 as a decoy target

Replicates Figure 4: the extent of PBP5/6 binding (mean fraction unbound over
0-60 min) against the rate of net influx and PBP access at 2x MIC. The paper
restricts this plot to the beta-lactams, because the BLIs were studied at a
fixed 4 mg/L rather than at 2x MIC.

```{r fig4, fig.height = 4.5, fig.width = 8}
# Rate_Influx/access at the studied concentration = scaled rate x C_drug
# (Table 1, last column). The scaled rate is read back out of each model's
# ini() block; the concentration is the CONC_<CODE>_MGL covariate value.
rate_influx <- lapply(seq_len(nrow(pbp_drugs)), function(i) {
  ini_tbl <- rxode2::rxode2(readModelDb(pbp_drugs$Model[i]))$iniDf
  scaled <- ini_tbl$est[match("rate_influx_scaled", ini_tbl$name)]
  tibble::tibble(
    Drug = pbp_drugs$Drug[i],
    rate_scaled = scaled,
    cdrug = pbp_drugs$Conc[i],
    rate_at_conc = scaled * pbp_drugs$Conc[i]
  )
}) |> bind_rows()

pbp56_extent <- sim_obs |>
  filter(PBP == "PBP5/6") |>
  group_by(Drug, Cells) |>
  # Linear trapezoidal AUC over 0-60 min divided by 60 = mean fraction unbound,
  # the same summary the paper computed in Phoenix WinNonlin.
  summarise(
    extent = sum(diff(time) * (head(frac_unbound, -1) + tail(frac_unbound, -1)) / 2) / 60,
    .groups = "drop"
  ) |>
  left_join(rate_influx, by = "Drug") |>
  left_join(select(pbp_drugs, Drug, Code, Class), by = "Drug") |>
  left_join(
    published |> filter(PBP == "PBP5/6") |> select(Drug, Cells, observed = auclast),
    by = c("Drug", "Cells")
  )

beta_lactams <- filter(pbp56_extent, !Class %in% c("BLI", "DBO-BLI"))

ggplot(beta_lactams, aes(rate_at_conc, extent)) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.4, colour = "grey40") +
  geom_point(aes(y = observed, shape = "Published (Tables 2-3)"), size = 2.6) +
  geom_point(aes(shape = "Model-predicted"), size = 2.6) +
  geom_text(aes(label = Code), size = 3, vjust = -0.9) +
  scale_shape_manual(values = c("Model-predicted" = 16, "Published (Tables 2-3)" = 1)) +
  facet_wrap(~Cells) +
  scale_x_log10() +
  labs(
    x = "Rate of net influx and PBP access at 2x MIC (molecules/min, log scale)",
    y = "Extent of PBP5/6 binding\n(mean fraction unbound, 0-60 min)",
    shape = NULL,
    title = "PBP5/6 acts as a decoy target for rapidly penetrating beta-lactams"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

```{r fig4-correlation}
beta_lactams |>
  group_by(Cells) |>
  summarise(
    `r2 (model-predicted)` = cor(log10(rate_at_conc), extent)^2,
    `r2 (published extents)` = cor(log10(rate_at_conc), observed)^2,
    .groups = "drop"
  ) |>
  knitr::kable(digits = 3, caption = "Correlation between log10(Rate_Influx/access at 2x MIC) and the extent of PBP5/6 binding, for the 11 beta-lactams. The paper reports r2 = 0.96 for this relationship (abstract; Fig. 4).")
```

Recomputing the correlation directly from the paper's own Table 2 / Table 3
values gives r2 = 0.94 in intact cells, close to the reported 0.96 (the paper
does not state the exact regression it fitted, so the small difference is not
diagnostic). The reproduced model gives a slightly weaker but qualitatively
identical relationship. In the lysed-cell panel ertapenem sits well off the line
-- the paper explicitly identifies it as the lysed-assay outlier, because only
the intact-cell assay sees ertapenem's poor penetration.

The drugs with the fastest net influx (imipenem, cefoxitin, doripenem) are the
ones that inactivate PBP5/6, exactly the inverse relationship the paper
describes: a beta-lactam that binds the highly expressed PBP5/6 must penetrate
rapidly to still reach the essential PBPs.

# PKNCA validation

The paper summarised PBP occupancy by the linear-trapezoidal AUC of the unbound
fraction over the first 60 min, divided by 60 (Tables 2 and 3, "Average +/- SD of
unbound PBPs over the first 60 min"). The same summary is computed here with
PKNCA on the four published sampling times, so the trapezoidal approximation
matches the paper's.

```{r pknca}
nca_input <- sim_obs |>
  filter(!is.na(frac_unbound)) |>
  transmute(
    Drug = as.character(Drug), Cells, PBP,
    time = time, frac = frac_unbound
  ) |>
  as.data.frame()

conc_obj <- PKNCA::PKNCAconc(nca_input, frac ~ time | Drug + Cells + PBP)
nca_res <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(
    conc_obj,
    intervals = data.frame(start = 0, end = 60, auclast = TRUE)
  )
)

nca_tbl <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD == "auclast") |>
  # AUC0-60 / 60 min = the paper's mean fraction of unbound PBPs
  mutate(PPORRES = PPORRES / 60)

head(nca_tbl)
```

# Comparison against the published NCA

Reference values are the Table 2 / Table 3 extents transcribed above.

```{r nca-comparison}
comparison <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_tbl,
  reference = published,
  by = c("Drug", "Cells", "PBP"),
  params = "auclast",
  units = c(auclast = "fraction unbound, mean over 0-60 min"),
  label_first_column = "NCA parameter"
)

knitr::kable(
  comparison, digits = 3,
  caption = "Model-predicted vs published extent of PBP occupancy for all 15 drugs, 6 PBPs and 2 assays (Tables 2 and 3)."
)
attr(comparison, "footnote")
```

`ncaComparisonTable()` returns formatted character columns, so the numeric
agreement is summarised from the underlying values.

```{r nca-comparison-summary}
agreement <- nca_tbl |>
  select(Drug, Cells, PBP, Simulated = PPORRES) |>
  inner_join(rename(published, Reference = auclast), by = c("Drug", "Cells", "PBP")) |>
  mutate(abs_diff = abs(Simulated - Reference))

agreement |>
  summarise(
    n = dplyr::n(),
    `median |difference|` = median(abs_diff),
    `90th pctile |difference|` = quantile(abs_diff, 0.9),
    `max |difference|` = max(abs_diff),
    `n within 0.10` = sum(abs_diff <= 0.10),
    `n within 0.20` = sum(abs_diff <= 0.20)
  ) |>
  knitr::kable(digits = 3, caption = "Agreement between the reproduced model and the published extent-of-binding values, on the fraction-unbound scale (range 0-1).")
```

The ten largest absolute deviations:

```{r nca-worst}
agreement |>
  arrange(desc(abs_diff)) |>
  head(10) |>
  knitr::kable(digits = 3, caption = "Largest deviations between the reproduced model and Tables 2-3.")
```

Across all 180 drug x PBP x assay combinations the median absolute deviation
from the published extents is about 0.05 on the 0-1 fraction scale, and no
combination deviates by more than about 0.25. That is the expected level of
agreement when a fitted model is regenerated from published population-mean
point estimates and compared against the *observed* means it was fit to: the
published extents carry replicate-to-replicate noise (the SDs in Tables 2 and 3
reach 0.3 for some entries) that a typical-value simulation cannot reproduce.
The largest deviations are concentrated in PBP4 and PBP5/6 for the slowly
penetrating penicillins and the BLIs, where the observed curves are flat near
1.0 and small differences in the fitted acylation constants move the predicted
plateau.

# Comparison against the published rate of net influx

Table 1 also reports each drug's rate of net influx and PBP access relative to
imipenem. That column is a direct function of the `rate_influx_scaled` values in
the model files and is reproduced exactly.

```{r rate-comparison}
rate_influx |>
  left_join(select(pbp_drugs, Drug, Code, Class), by = "Drug") |>
  mutate(`Relative to IPM` = rate_scaled[Drug == "imipenem"] / rate_scaled) |>
  select(Drug, Code, Class, `Scaled rate (molecules/min per mg/L)` = rate_scaled,
         `C_drug (mg/L)` = cdrug,
         `Rate at studied conc. (molecules/min)` = rate_at_conc,
         `Relative to IPM`) |>
  knitr::kable(digits = 3, caption = "Reproduces Table 1: scaled rate of net influx and PBP access, the rate at the studied concentration, and the fold-difference relative to imipenem.")
```

Because the reference values are fractions bounded in `[0, 1]`, the percentage
difference reported by `ncaComparisonTable()` is inflated wherever the reference
value is small; the absolute difference on the fraction scale is the meaningful
metric and is summarised above. Rows flagged with `*` are dominated by such
small-denominator cases.

# Mechanistic checks

## Mass balance of PBP molecules

Every PBP molecule is either unbound (a state) or acylated. The number bound
must equal the initial number minus the current number, and can never exceed the
initial number or fall below zero.

```{r mass-balance-pbp}
mb <- sim_dense |>
  group_by(Drug, Cells, PBP) |>
  summarise(
    start = first(n_unbound[time == 0]),
    min_unbound = min(n_unbound),
    max_unbound = max(n_unbound),
    .groups = "drop"
  )

stopifnot(
  all(mb$min_unbound >= -1e-8),               # no negative PBP counts
  all(mb$max_unbound <= mb$start + 1e-8)      # monotone non-increasing
)
tibble::tibble(
  Check = c("No negative unbound-PBP counts", "Unbound PBPs never exceed the initial condition"),
  Result = c("pass", "pass")
) |> knitr::kable()
```

## Mass balance of drug molecules in periplasm

In the intact-cell assay the cumulative drug that entered the periplasm must
equal the free periplasmic pool plus the drug consumed by acylation. Because one
drug molecule is consumed per PBP acylated, the consumed amount equals the total
number of PBP molecules inactivated.

```{r mass-balance-drug}
drug_balance <- lapply(seq_len(nrow(pbp_drugs)), function(i) {
  s <- solve_one(i, dense_times, 1)
  pbp_cols <- c("npbp1a", "npbp1b", "npbp2", "npbp3", "npbp4", "npbp56")
  bound <- sum(s[s$time == 0, pbp_cols]) - rowSums(s[, pbp_cols])
  tibble::tibble(
    Drug = pbp_drugs$Drug[i],
    time = s$time,
    entered = s$rate_influx * s$time,   # constant influx rate
    accounted = s$nperi + bound
  )
}) |> bind_rows()

worst <- drug_balance |>
  filter(time > 0) |>
  mutate(rel_err = abs(entered - accounted) / pmax(entered, 1e-12)) |>
  summarise(`max relative mass-balance error` = max(rel_err))

stopifnot(worst[[1]] < 1e-4)
knitr::kable(worst, digits = 8, caption = "Periplasmic drug mass balance in the intact-cell assay: cumulative influx equals free periplasmic drug plus drug consumed by PBP acylation.")
```

## Saturation regime

The paper fixed `Km` at 1,000 drug molecules to keep periplasmic binding in the
linear (non-saturated) regime in intact cells, while the 10^7 molecules used in
the lysed assay saturate the reaction 10,000-fold above `Km`.

```{r saturation}
sim_peri <- lapply(seq_len(nrow(pbp_drugs)), function(i) {
  s <- solve_one(i, dense_times, 1)
  tibble::tibble(Drug = pbp_drugs$Drug[i], time = s$time, nperi = s$nperi)
}) |> bind_rows()

# While PBP binding is still ongoing (>= 5% of PBP5/6 still unbound), N_peri
# should stay well below Km = 1000; once the PBPs are exhausted N_peri simply
# accumulates because the model has no periplasmic efflux term.
binding_window <- sim_dense |>
  filter(PBP == "PBP5/6", Cells == "Intact", frac_unbound > 0.05) |>
  group_by(Drug) |>
  summarise(t_end = max(time), .groups = "drop")

sim_peri |>
  left_join(binding_window, by = "Drug") |>
  filter(!is.na(t_end), time <= t_end) |>
  group_by(Drug) |>
  summarise(`max N_peri while PBP5/6 binding` = max(nperi), .groups = "drop") |>
  mutate(`below Km = 1000` = `max N_peri while PBP5/6 binding` < 1000) |>
  knitr::kable(digits = 1, caption = "Periplasmic drug molecules per cell during the PBP-binding window in intact cells.")
```

# Assumptions and deviations

* **Main-text Eq 2/3/4 vs Fig. S8 code.** The published equations and the
  published estimation code express the Michaelis-Menten term differently. Main
  text Eq 2 defines `MM_1a = Vmax_1a / (Km + N_peri)` and Eq 3 then multiplies by
  `N_peri` (`d(N_PBP1a)/dt = -MM_1a x N_peri x N_PBP1a`), whereas Eq 4 omits that
  `N_peri` factor from the periplasmic balance. Fig. S8 instead folds `N_peri`
  into `MM` (`MM1a = K1a*Nperi/(KM_1a+Nperi)`) and uses the same `MM_i x N_PBPi`
  product in both the PBP ODEs and the periplasmic balance. The two forms give
  identical PBP ODEs, but only the Fig. S8 form conserves drug mass (each
  acylation consumes exactly one drug molecule). **The Fig. S8 code form is the
  one implemented here**, and the mass-balance check above confirms it.
* **No between-replicate variability.** The supplement states that all
  parameters were described by log-normal (or logistic, for the noise terms)
  distributions, but Tables S1 and S2 report population means only -- no variance
  components are published. The model files therefore contain no `eta` terms and
  simulate typical values. Residual error *is* reported and is implemented.
* **Residual error form.** Fig. S8 lines 111-116 give
  `V = (SDin + SDsl * Y)^2`, i.e. the additive and proportional standard
  deviations add linearly. That is nlmixr2's `combined1()` form, not the default
  `combined2()` root-sum-of-squares.
* **Choice of estimation dataset for the nuisance parameters.** The 15 drugs
  were fit in five datasets of four drugs each, with imipenem as the shared
  backbone (Table S2). Each drug's `noise_*`, `fini_*`, `addSd_*` and `propSd_*`
  values are taken from its own dataset. Cefoxitin appeared in datasets 2 and 5;
  the supplement states the final cefoxitin estimates come from dataset 5, which
  is what `LopezArguello_2023_cefoxitin_qsp` uses. **Imipenem is the one
  ambiguous case**: it appears in all five datasets and Table S1 reports a single
  set of imipenem acylation rate constants without attributing them to a dataset.
  `LopezArguello_2023_imipenem_qsp` uses dataset 1 (the carbapenem dataset), as
  directed by the operator (sidecar `request-001` q2, answer A). Only the
  background-noise, initial-condition and residual-error parameters are affected;
  the drug-specific acylation constants and net-influx rate are unambiguous.
* **`CELLS_INTACT` and `CONC_<CODE>_MGL` as covariates.** The assay arm and the
  studied extracellular concentration are experimental design variables, not
  estimated quantities, and are supplied on the event table rather than as
  `ini()` entries. `CELLS_INTACT` is a new canonical covariate ratified by the
  operator for this extraction (sidecar `request-001` q1, answer A); it is
  general-scope because the intact-versus-lysed contrast is a standard
  target-site-penetration design rather than something specific to this paper.
  It is structural rather than a coefficient multiplier: it gates both the influx
  term and the periplasmic initial condition, so it has no `e_<cov>_<param>`
  coefficient. The concentration covariates are well-formed members of the
  existing `CONC_<DRUG>_MGL` family (`CONC_IPM_MGL` was already registered for
  Landersdorfer 2018); the other fourteen were added alongside this extraction.
* **Shared parameters for relebactam.** Table 1 footnote c and Table S1 footnote
  b state that relebactam's rate of net influx and its PBP2 acylation constant
  were estimated as shared parameters with avibactam because relebactam bound
  PBPs only to a limited extent in intact bacteria. These are *shared estimates*,
  not fixed constants, and are therefore not wrapped in `fixed()`; the sharing is
  recorded in the in-file source-trace comments. The relebactam estimates carry
  extra uncertainty.
* **Total PBP molecules per cell.** No published count exists for
  *P. aeruginosa*, so the authors borrowed the *E. coli* figure of 1,731
  molecules per cell and split it using *P. aeruginosa* relative band intensities
  (Table S3). This is the authors' documented assumption, carried through
  unchanged.
* **Ceftazidime PBP5/6 at 60 min.** The paper notes that ceftazidime showed
  considerable PBP5/6 binding at 60 min but not earlier, "which may have been an
  'outlier' that could not be explained by QSP modeling" (Results, Time-course
  PBP binding in whole cells). The reproduced model shows the same behaviour --
  little PBP5/6 binding -- so this row differs from the observed Table 2 value by
  design, not by extraction error.
* **No dosing events.** The extracellular concentration is held static
  throughout the assay, so the model has no dose records; drug enters the
  periplasm at a constant rate while the intact-cell switch is on. The model has
  no periplasmic efflux or degradation term, so once the PBPs are exhausted
  `N_peri` accumulates linearly. That is a faithful consequence of the published
  structure and does not affect the 0-60 min binding predictions.
* **Negative-state guard.** Fig. S8 lines 94-100 clamp any negative state to
  zero before computing the output. This is reproduced as `max(state, 0)` in the
  output equations. With the published parameter values the states remain
  non-negative, so the guard never activates.

# Session information

```{r session}
sessionInfo()
```
