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

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Liang_2023_polymyxinB"))$reference`
* Description: `r rxode2::rxode(readModelDb("Liang_2023_polymyxinB"))$description`
* Article: <https://doi.org/10.3389/fphar.2023.1122310>
* PubMed Central: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10090446/>

Liang et al. fitted a two-compartment intravenous model to steady-state
polymyxin B concentrations from critically ill adults using the nonparametric
adaptive grid (NPAG) algorithm in Pmetrics 1.9.7. Serum albumin was retained as
a power covariate on clearance and age as a power covariate on the central
volume of distribution; the resulting model (Table 2, "model 6") was then used
in Monte Carlo simulations to recommend a dosing regimen.

# Population

Twenty-two critically ill adults (age 18 years or older) treated with
intravenous polymyxin B at a single center in Guangzhou, China between November
2021 and September 2022 contributed 64 plasma concentrations (Table 1, Results
3.1). The cohort was 77.3% male, with median age 68 years (range 31-94), median
weight 60 kg (range 50-80), median APACHE-II score 21.5 points (range 15-46),
median baseline serum creatinine 75 umol/L (range 33-125), median baseline
Cockcroft-Gault creatinine clearance 68.29 mL/min (range 34.05-192.59), and
median baseline serum albumin 31.45 g/L (range 23.1-41.5).

Half the cohort had pulmonary infection alone and a further 45.5% had pulmonary
infection with sepsis or septicemia; *Acinetobacter baumannii* was the most
common pathogen (50.0%). Patients receiving continuous renal replacement therapy
(n = 10) or extracorporeal membrane oxygenation (n = 2) at the time of blood
collection were excluded, so the model must not be extrapolated to those
settings (Discussion, study limitations).

Dosing was a loading dose of 100-150 mg followed by a maintenance dose of
50-75 mg q12h, infused over 1-2 h. Sampling began at least 48 h after the start
of treatment: 24 trough samples, 23 peak samples, and 17 samples 6-8 h after the
end of infusion. Total (not unbound) polymyxin B was measured as the sum of
polymyxin B1, B2, and B1-Ile by LC-MS/MS.

The same information is available programmatically via the model's `population`
metadata:

```{r population}
str(readModelDb("Liang_2023_polymyxinB")()$population)
```

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `d/dt(central)` | n/a | Equation 1: `dX1/dt = RateIV + (Q/Vp)*X2 - (CL/Vc + Q/Vc)*X1` |
| `d/dt(peripheral1)` | n/a | Equation 2: `dX2/dt = (Q/Vc)*X1 - (Q/Vp)*X2` |
| `lcl` | `log(1.22)` | Table 3, CL median (mean 1.24, SD 0.38 L/h) |
| `lvc` | `log(13.52)` | Table 3, Vc median (mean 16.64, SD 12.74 L) |
| `lq` | `log(2.34)` | Table 3, Q median (mean 3.04, SD 2.27 L/h) |
| `lvp` | `log(76.53)` | Table 3, Vp median (mean 66.20, SD 36.25 L) |
| `e_alb_cl` | `fixed(-0.95)` | Table 2, model 6: `CL0*(ALB/31.45)^(-0.95)` |
| `e_age_vc` | `fixed(0.95)` | Table 2, model 6: `V0*(age/68)^0.95` |
| ALB reference 31.45 g/L | n/a | Table 1, median baseline ALB |
| AGE reference 68 years | n/a | Table 1, median age |
| `etalcl` | `0.088837` | Table 3, CL between-subject variability 30.48 %CV |
| `etalvc` | `0.461209` | Table 3, Vc between-subject variability 76.55 %CV |
| `etalq` | `0.444176` | Table 3, Q between-subject variability 74.78 %CV |
| `etalvp` | `0.262261` | Table 3, Vp between-subject variability 54.76 %CV |
| `addSd` | `fixed(0.072)` | Methods 2.4 (`C0 = 0.1`) x Results 3.2 (final gamma 0.72) |
| `propSd` | `fixed(0.108)` | Methods 2.4 (`C1 = 0.15`) x Results 3.2 (final gamma 0.72) |

Two transcription steps deserve comment.

**Between-subject variability.** Pmetrics NPAG returns a discrete nonparametric
joint distribution rather than a parametric OMEGA. Table 3 summarises it by mean,
median, SD, and "%CV", and the %CV column is exactly `SD / mean` for all four
parameters (for example CL: `0.38 / 1.24 = 30.6%`). Encoding this in nlmixr2
requires a parametric approximation; a log-normal is used, with the Table 3
**median** mapped onto `exp(l<param>)` (the log-normal median) and the variance
obtained from `omega^2 = log(CV^2 + 1)`.

**Residual error.** Pmetrics weights each observation by the reciprocal of
`SD x gamma`, where `SD = C0 + C1*Cobs + C2*Cobs^2 + C3*Cobs^3` with
`C0 = 0.1`, `C1 = 0.15`, `C2 = C3 = 0` (Methods 2.4) and the final-cycle
`gamma` was 0.72 (Results 3.2). The total residual standard deviation is
therefore `0.72 * (0.1 + 0.15*C) = 0.072 + 0.108*C` mg/L, which is a combined
additive plus proportional error.

# What drives exposure in this model

Before simulating, it is worth writing down what the final model implies. At
steady state, the area under the concentration-time curve over 24 h is fixed by
clearance alone:

$$
\mathrm{AUC}_{ss,24h} = \frac{\text{total daily dose}}{CL}
= \frac{\text{daily dose}}{CL_0 \cdot (ALB/31.45)^{-0.95}}
$$

Age enters the model only through `Vc`, and `Vc` does not appear in this
expression. **Under the published model, age changes the shape of the
concentration-time curve (peak and trough) but cannot change AUC, and therefore
cannot change any AUC-based probability of target attainment.** This is
consistent with the source's own Figure 3, in which the three panels of each
row (same albumin stratum, ages 34 / 68 / 93) are visually indistinguishable,
but it is *not* consistent with the paper's Conclusion that "the possibility of
obtaining an AUC/MIC >= 50 decreased with older age". See the Assumptions and
deviations section.

The age-invariance is checked explicitly below.

# Virtual cohort

Original observed data are not publicly available. The cohorts below are virtual
populations whose covariate distributions follow the three albumin strata and the
five dosing regimens defined in Methods 2.5.

```{r cohort}
set.seed(20230329)

n_per_arm <- 100L

alb_strata <- tibble::tribble(
  ~alb_group,        ~alb_lo, ~alb_hi,
  "Ultra-low ALB",      23.1,    24.9,
  "Low ALB",            25.0,    34.9,
  "Normal ALB",         35.0,    41.5
)

regimens <- tibble::tibble(
  regimen    = c("100 mg q12h", "75 mg q12h", "60 mg q12h",
                 "50 mg q12h", "40 mg q12h"),
  dose_mg    = c(100, 75, 60, 50, 40)
) |>
  mutate(daily_mg = dose_mg * 2)

arms <- tidyr::crossing(alb_strata, regimens) |>
  mutate(arm = paste(alb_group, regimen, sep = " | "))

# 30 days of q12h dosing takes the model well past steady state (terminal
# half-life is about 71 h at the typical parameter values), so the final 24 h
# window is a true steady-state interval.
n_doses  <- 60L
tau      <- 12
obs_from <- (n_doses - 2L) * tau   # start of the final 24 h window
obs_to   <- obs_from + 24

# `n` subjects with consecutive IDs starting at `id_offset + 1`. IDs must be
# consecutive across the whole event table: rxSolve re-indexes non-consecutive
# IDs in its output, which silently breaks the id-based join PKNCA makes
# between the concentration and dose objects.
make_arm <- function(n, dose_mg, alb_lo, alb_hi, arm, alb_group, regimen,
                     daily_mg, id_offset = 0L, age = 68) {
  subj <- tibble::tibble(
    id        = id_offset + seq_len(n),
    ALB       = stats::runif(n, alb_lo, alb_hi),
    AGE       = age,
    arm       = arm,
    alb_group = alb_group,
    regimen   = regimen,
    daily_mg  = daily_mg
  )

  doses <- subj |>
    mutate(
      time = 0, amt = dose_mg, evid = 1L, dur = 1,
      ii = tau, addl = n_doses - 1L, cmt = "central"
    )

  obs <- subj |>
    tidyr::crossing(time = seq(obs_from, obs_to, by = 0.25)) |>
    mutate(
      amt = NA_real_, evid = 0L, dur = NA_real_,
      ii = NA_real_, addl = NA_integer_, cmt = "central"
    )

  bind_rows(doses, obs) |> arrange(id, time, desc(evid))
}

events <- do.call(
  bind_rows,
  lapply(seq_len(nrow(arms)), function(i) {
    a <- arms[i, ]
    make_arm(
      n         = n_per_arm,
      dose_mg   = a$dose_mg,
      alb_lo    = a$alb_lo,
      alb_hi    = a$alb_hi,
      arm       = a$arm,
      alb_group = a$alb_group,
      regimen   = a$regimen,
      daily_mg  = a$daily_mg,
      id_offset = (i - 1L) * n_per_arm
    )
  })
)

# Disjoint IDs across arms are mandatory; duplicate IDs silently merge subjects.
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

# ... and they must also be consecutive from 1. rxSolve re-indexes
# non-consecutive IDs in its output, which leaves the simulated concentrations
# carrying different IDs than the event table the PKNCA dose object is built
# from -- the join then fails silently and produces scrambled NCA values.
assert_consecutive_ids <- function(x) {
  stopifnot(identical(sort(unique(x$id)), seq_len(dplyr::n_distinct(x$id))))
  invisible(x)
}
assert_consecutive_ids(events)

events <- events |>
  mutate(
    alb_group = factor(alb_group, levels = alb_strata$alb_group),
    regimen   = factor(regimen, levels = regimens$regimen)
  )

nrow(events)
dplyr::n_distinct(events$id)
```

# Simulation

```{r simulate}
mod <- readModelDb("Liang_2023_polymyxinB")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("arm", "alb_group", "regimen", "daily_mg", "ALB", "AGE")
) |>
  as.data.frame()

sim <- sim |> mutate(time_in_window = time - obs_from)

summary(sim$Cc)
```

# Steady-state concentration-time profiles

```{r profiles}
sim |>
  filter(!is.na(Cc)) |>
  group_by(alb_group, regimen, time_in_window) |>
  summarise(
    Q05 = quantile(Cc, 0.05),
    Q50 = quantile(Cc, 0.50),
    Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  ggplot(aes(time_in_window, Q50, colour = regimen, fill = regimen)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.15, colour = NA) +
  geom_line() +
  facet_wrap(~alb_group) +
  labs(
    x = "Time within the steady-state 24 h window (h)",
    y = "Polymyxin B concentration (mg/L)",
    colour = "Regimen", fill = "Regimen",
    title = "Simulated steady-state profiles by albumin stratum",
    caption = paste(
      "Median with 5th-95th percentile band,", n_per_arm,
      "subjects per arm. Model: Liang 2023 Table 2 model 6."
    )
  ) +
  theme(legend.position = "bottom")
```

Lower albumin gives higher clearance (exponent -0.95), so the ultra-low albumin
stratum has the lowest exposure for a given regimen -- the mechanism behind the
paper's recommendation that low-albumin patients need higher doses.

# Age-invariance check

The model's structure implies that age shifts `Vc` but leaves AUC untouched.
A typical-value simulation at the three ages used in the paper's Figure 3
confirms this.

```{r age-check}
mod_typ <- mod |> rxode2::zeroRe()

# Every typical-value solve below passes `omega = NA`. `zeroRe()` alone is NOT
# sufficient once a stochastic `rxSolve()` has already run in the same session:
# rxode2 keeps the previous solve's omega in its solve options, so a later call
# re-samples etas even though `mod_typ$omega` is a zero matrix -- silently
# turning a "typical value" run back into a one-subject-per-arm random draw.
# `omega = NA` is the explicit "no random effects" sentinel. The stopifnot()
# checks below fail loudly if this ever regresses.

ages <- c(34, 68, 93)

# One typical subject per age, consecutive IDs.
age_events <- do.call(
  bind_rows,
  lapply(seq_along(ages), function(i) {
    make_arm(
      n         = 1L,
      dose_mg   = 75,
      alb_lo    = 31.45, alb_hi = 31.45,
      arm       = paste0("age ", ages[i]),
      alb_group = "Cohort median ALB",
      regimen   = "75 mg q12h",
      daily_mg  = 150,
      id_offset = i - 1L,
      age       = ages[i]
    )
  })
)

assert_consecutive_ids(age_events)

age_sim <- rxode2::rxSolve(
  mod_typ, events = age_events, keep = c("arm", "AGE", "daily_mg"),
  omega = NA
) |>
  as.data.frame() |>
  filter(!is.na(Cc))

# CL must not vary at all in a typical-value run (ALB is identical here).
stopifnot(dplyr::n_distinct(round(age_sim$cl, 8)) == 1L)

age_summary <- age_sim |>
  group_by(AGE) |>
  summarise(
    `Cmax (mg/L)`  = max(Cc),
    `Cmin (mg/L)`  = min(Cc),
    `AUC0-24 (mg*h/L)` =
      sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups = "drop"
  ) |>
  rename("Age (years)" = AGE)

# The structural claim under test: AUC is age-invariant. Anything above
# trapezoidal-integration noise (< 0.1%) means the covariate model is wrong.
stopifnot(
  diff(range(age_summary$`AUC0-24 (mg*h/L)`)) <
    0.001 * mean(age_summary$`AUC0-24 (mg*h/L)`)
)

knitr::kable(
  age_summary,
  digits = 2,
  caption = paste(
    "Typical-subject steady-state exposure at ALB = 31.45 g/L,",
    "75 mg q12h. Cmax and Cmin change with age; AUC0-24 does not."
  )
)
```

`Cmax` and `Cmin` move with age (larger `Vc` in older patients flattens the
profile), while `AUC0-24` is identical to numerical precision across all three
ages, exactly as the closed-form expression predicts.

# PKNCA validation

```{r pknca}
sim_nca <- sim |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, arm)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | arm + id,
  concu = "mg/L", timeu = "h"
)

dose_df <- events |>
  filter(evid == 1) |>
  transmute(id, time = obs_from, amt = amt * 2, arm)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id, doseu = "mg")

intervals <- data.frame(
  start   = obs_from,
  end     = obs_to,
  cmax    = TRUE,
  cmin    = TRUE,
  auclast = TRUE,
  cav     = TRUE
)

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

The dose record is placed at the start of the steady-state window and carries
the **daily** dose so that `auclast` over the 24 h interval is directly
comparable with the paper's `AUCss,24h`.

## Comparison against the model-implied steady-state exposure

Liang et al. report no NCA table, so there is no published Cmax / AUC to compare
against. What can be checked exactly is whether the packaged model reproduces the
closed-form steady-state exposure implied by its own published parameters:
`AUCss,24h = daily dose / CL` with `CL = 1.22 * (ALB/31.45)^(-0.95)` (Table 3
median CL, Table 2 model 6 covariate equation). A typical-value simulation at the
midpoint albumin of each stratum isolates the structural and covariate model from
between-subject variability.

```{r nca-typical}
alb_mid <- alb_strata |>
  mutate(ALB_mid = (alb_lo + alb_hi) / 2)

typ_arms <- tidyr::crossing(alb_mid, regimens) |>
  mutate(arm = paste(alb_group, regimen, sep = " | "))

# A typical-value run needs only one subject per arm; IDs are consecutive.
typ_events <- do.call(
  bind_rows,
  lapply(seq_len(nrow(typ_arms)), function(i) {
    a <- typ_arms[i, ]
    make_arm(
      n         = 1L,
      dose_mg   = a$dose_mg,
      alb_lo    = a$ALB_mid, alb_hi = a$ALB_mid,
      arm       = a$arm,
      alb_group = a$alb_group,
      regimen   = a$regimen,
      daily_mg  = a$daily_mg,
      id_offset = i - 1L
    )
  })
)

assert_consecutive_ids(typ_events)

typ_sim <- rxode2::rxSolve(
  mod_typ, events = typ_events,
  keep = c("arm", "alb_group", "regimen", "daily_mg", "ALB"),
  omega = NA
) |>
  as.data.frame() |>
  filter(!is.na(Cc))

# In a typical-value run CL is a deterministic function of ALB alone, so each
# distinct ALB must map to exactly one CL.
stopifnot(
  typ_sim |>
    group_by(ALB) |>
    summarise(n_cl = dplyr::n_distinct(round(cl, 8)), .groups = "drop") |>
    pull(n_cl) |>
    max() == 1L
)

typ_nca <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(
    PKNCA::PKNCAconc(
      typ_sim |> select(id, time, Cc, arm),
      Cc ~ time | arm + id, concu = "mg/L", timeu = "h"
    ),
    PKNCA::PKNCAdose(
      typ_events |> filter(evid == 1) |>
        transmute(id, time = obs_from, amt = amt * 2, arm),
      amt ~ time | arm + id, doseu = "mg"
    ),
    intervals = data.frame(
      start = obs_from, end = obs_to, auclast = TRUE, cav = TRUE
    )
  )
)

published_implied <- typ_arms |>
  mutate(
    cl_typ  = 1.22 * (ALB_mid / 31.45)^(-0.95),
    auclast = daily_mg / cl_typ,
    cav     = auclast / 24
  ) |>
  select(arm, auclast, cav)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = typ_nca,
  reference = published_implied,
  by        = "arm",
  units     = c(auclast = "mg*h/L", cav = "mg/L"),
  # This is a deterministic identity check, not a noisy comparison against a
  # published NCA table, so the tolerance is much tighter than the usual 20%.
  tolerance_pct = 5
)

knitr::kable(
  cmp,
  digits  = 2,
  caption = paste(
    "Simulated typical-value steady-state exposure vs. the closed-form",
    "value implied by the published CL and albumin covariate equation.",
    "* differs from reference by >5%."
  ),
  align = c("l", "l", "r", "r", "r")
)

attr(cmp, "footnote")
```

Every row agrees to better than 0.2% (the residual difference is trapezoidal
integration error on the 0.25 h observation grid), confirming that the ODE
system, the albumin power covariate, and the unit conventions in the packaged
model reproduce the published parameterisation exactly.

# Replicate Figure 3 -- PTA of AUC/MIC >= 50

Methods 2.5 sets the efficacy target as `AUC/MIC >= 50` (derived from a
consensus `fAUC/MIC` target of about 20 and an unbound fraction of 0.42) over
the MIC range 0.25-4 mg/L.

```{r figure-3}
mic_grid <- c(0.25, 0.5, 1, 2, 4)

auc_subject <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD == "auclast") |>
  select(arm, id, auc24 = PPORRES) |>
  left_join(
    arms |> select(arm, alb_group, regimen),
    by = "arm"
  ) |>
  mutate(
    alb_group = factor(alb_group, levels = alb_strata$alb_group),
    regimen   = factor(regimen, levels = regimens$regimen)
  )

pta_efficacy <- tidyr::crossing(auc_subject, MIC = mic_grid) |>
  group_by(alb_group, regimen, MIC) |>
  summarise(PTA = mean(auc24 / MIC >= 50), .groups = "drop")

ggplot(pta_efficacy, aes(MIC, PTA, colour = regimen)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0.9, linetype = "dashed") +
  scale_x_log10(breaks = mic_grid, labels = mic_grid) +
  scale_y_continuous(limits = c(0, 1)) +
  facet_wrap(~alb_group) +
  labs(
    x = "MIC (mg/L)", y = "PTA",
    colour = "Regimen",
    title = "Probability of target attainment for AUC/MIC >= 50",
    caption = paste(
      "Replicates Figure 3 of Liang 2023 (rows A-C / D-F / G-I; the",
      "age panels within each row are identical under the published model).",
      "Dashed line: PTA = 0.9."
    )
  ) +
  theme(legend.position = "bottom")
```

The paper states its Figure 3 conclusions in text; those statements are the
quantitative reference, and the simulated PTA values at MIC = 1 mg/L are
compared against them below.

```{r figure-3-table}
pta_mic1 <- pta_efficacy |>
  filter(MIC == 1) |>
  tidyr::pivot_wider(names_from = regimen, values_from = PTA) |>
  select(-MIC)

knitr::kable(
  pta_mic1 |> rename("Albumin stratum" = alb_group),
  digits  = 3,
  caption = paste(
    "Simulated PTA of AUC/MIC >= 50 at MIC = 1 mg/L.",
    "Liang 2023 recommends 100-75 mg q12h for the ultra-low stratum,",
    "100-60 mg q12h for the low stratum, and 100-50/100-60 mg q12h",
    "for the normal stratum."
  )
)
```

```{r figure-3-claims}
claims <- tibble::tribble(
  ~Claim,                                                          ~Source,
  "All regimens reach PTA >= 0.9 at MIC <= 0.5 mg/L",              "Results 3.4",
  "Ultra-low ALB, MIC = 1: 100 and 75 mg q12h reach PTA >= 0.9",   "Results 3.4",
  "Low ALB, MIC = 1: 100, 75 and 60 mg q12h reach PTA >= 0.9",     "Results 3.4",
  "Normal ALB, MIC = 1: down to 50-60 mg q12h reaches PTA >= 0.9", "Results 3.4",
  "No regimen reaches PTA >= 0.9 at MIC = 4 mg/L",                 "Discussion"
)

check <- tibble::tibble(
  Claim = claims$Claim,
  Source = claims$Source,
  Simulated = c(
    sprintf(
      "min PTA at MIC <= 0.5 = %.3f",
      min(pta_efficacy$PTA[pta_efficacy$MIC <= 0.5])
    ),
    sprintf(
      "100 mg %.3f, 75 mg %.3f, 60 mg %.3f",
      pta_mic1$`100 mg q12h`[1], pta_mic1$`75 mg q12h`[1],
      pta_mic1$`60 mg q12h`[1]
    ),
    sprintf(
      "100 mg %.3f, 75 mg %.3f, 60 mg %.3f, 50 mg %.3f",
      pta_mic1$`100 mg q12h`[2], pta_mic1$`75 mg q12h`[2],
      pta_mic1$`60 mg q12h`[2], pta_mic1$`50 mg q12h`[2]
    ),
    sprintf(
      "60 mg %.3f, 50 mg %.3f, 40 mg %.3f",
      pta_mic1$`60 mg q12h`[3], pta_mic1$`50 mg q12h`[3],
      pta_mic1$`40 mg q12h`[3]
    ),
    sprintf(
      "max PTA at MIC = 4 = %.3f",
      max(pta_efficacy$PTA[pta_efficacy$MIC == 4])
    )
  )
)

knitr::kable(
  check,
  caption = "Published Figure 3 / Results 3.4 claims vs. the simulated PTA."
)
```

Every published claim is reproduced in direction and in the regimens it names.
The simulated PTA is systematically a little *more* permissive than the paper's
at the low-dose end -- 60 mg q12h clears 0.9 in the ultra-low stratum here, and
40 mg q12h clears it in the normal stratum, where the paper stops one step
higher. This is the expected consequence of approximating the nonparametric CL
distribution by a log-normal, which has thinner tails; PTA is a
threshold-crossing probability and so is driven by exactly those tails. No
parameter was adjusted to close the gap.

# Replicate Figure 4 -- PTA of AUCss,24h > 100 (nephrotoxicity)

Methods 2.5 uses `AUCss,24h > 100 mg*h/L` as the nephrotoxicity limit.

```{r figure-4}
pta_tox <- auc_subject |>
  group_by(alb_group, regimen) |>
  summarise(PTA = mean(auc24 > 100), .groups = "drop")

ggplot(pta_tox, aes(regimen, PTA, fill = alb_group)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  geom_hline(yintercept = 0.9, linetype = "dashed") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    x = NULL, y = "PTA",
    fill = "Albumin stratum",
    title = "Probability of exceeding the nephrotoxicity limit AUCss,24h > 100",
    caption = paste(
      "Replicates Figure 4 of Liang 2023. Under the published model this",
      "quantity depends on albumin only, not on age. Dashed line: PTA = 0.9."
    )
  ) +
  theme(legend.position = "bottom")
```

```{r figure-4-table}
knitr::kable(
  pta_tox |>
    tidyr::pivot_wider(names_from = regimen, values_from = PTA) |>
    rename("Albumin stratum" = alb_group),
  digits  = 3,
  caption = paste(
    "Simulated PTA of AUCss,24h > 100 mg*h/L. Liang 2023 identifies",
    "100 mg q12h as the toxicity-limit dose for patients with normal",
    "albumin levels."
  )
)
```

The simulated toxicity PTA rises with albumin (higher albumin gives lower
clearance and therefore higher exposure) and is highest for 100 mg q12h,
reproducing the paper's identification of 100 mg q12h as the toxicity-limiting
regimen for normal-albumin patients. The simulated absolute probabilities are
higher than the paper's Figure 4 for the low-dose regimens; see Assumptions and
deviations.

# Assumptions and deviations

* **Nonparametric distribution approximated by a log-normal.** Pmetrics NPAG
  returns a discrete nonparametric joint parameter distribution, which nlmixr2
  cannot represent directly. The Table 3 **median** of each parameter is mapped
  onto `exp(l<param>)` (the log-normal median) and the Table 3 "%CV" is
  converted via `omega^2 = log(CV^2 + 1)`. This mapping reproduces the paper's
  reported population *means* well for CL (1.28 implied vs 1.24 reported), Q
  (2.92 vs 3.04), and Vc (17.03 vs 16.64). It does not for Vp, where the
  nonparametric distribution is left-skewed (median 76.53 > mean 66.20) and no
  log-normal can match both summaries; the implied mean Vp is 87.3 vs 66.20
  reported. Vp affects only the distribution phase, not steady-state AUC, so
  this does not propagate to any of the exposure targets validated above.
  The Table 3 mean +/- SD values quoted in the paper's Abstract and Discussion
  are recorded in the model file comments alongside each value.
* **Parameter correlations are not modelled.** Table 3 reports no correlations
  between CL, Q, Vc, and Vp, so the etas are diagonal. A nonparametric joint
  distribution generally does carry correlation structure that is lost here.
* **Covariate exponents encoded as `fixed()`.** The paper never states whether
  the +/-0.95 exponents were estimated or imposed, and reports no SE, RSE, or CI
  for either. They are encoded as `fixed()` because Table 2's
  information-criterion arithmetic settles the question: `AIC - (-2LL)` is
  11.0-11.1 for every two-compartment model (models 2 through 6) regardless of
  whether zero, one, or two covariate relationships are present, so adding a
  covariate added no estimated parameter. Table 3 correspondingly reports
  estimates for CL, Q, Vc, and Vp only.
* **Residual error is `fixed()`.** `C0` and `C1` were "scientifically set to 0.1,
  0.15" a priori (Methods 2.4) and the final `gamma` of 0.72 is reported without
  uncertainty (Results 3.2). Pmetrics scales the assay-error polynomial by the
  *observed* concentration; nlmixr2's `prop()` uses the *predicted*
  concentration. The two coincide at the fit and differ only in how residuals are
  weighted during estimation.
* **Age cannot affect AUC under the published model.** Age enters only through
  `Vc`, so `AUCss,24h = daily dose / CL` is independent of age. The paper's
  Figure 3 is consistent with this (the three age panels within each albumin row
  are visually indistinguishable), but the Conclusion's claim that "the
  possibility of obtaining an AUC/MIC >= 50 decreased with older age" is not
  supported by the model as published. The vignette therefore collapses each
  Figure 3 row to a single panel and verifies the age-invariance explicitly. The
  model file encodes the published structure without modification.
* **Albumin stratum bounds.** Methods 2.5 gives the ultra-low albumin group as
  "21.3-24.9 g/L", but Table 1 reports a cohort minimum albumin of 23.1 g/L and
  the Figure 3 caption gives the same group as "23.1-24.9 g/L". The Figure 3
  caption bound (23.1 g/L) is used here; the Methods value appears to be a
  typographical error.
* **Virtual covariate distributions.** Albumin is drawn uniformly within each
  stratum, because the paper does not report the within-stratum albumin
  distribution. Age is fixed at the cohort median of 68 years for the PTA
  cohorts (immaterial, per the age-invariance above) and set to 34 / 68 / 93
  years in the dedicated age-invariance check, matching the paper's 5th
  percentile, median, and 95th percentile.
* **No loading dose in the simulations.** The paper's clinical regimens included
  a 100-150 mg loading dose, but Methods 2.5 defines the Monte Carlo regimens as
  maintenance dosing only. The simulations here run 60 q12h doses (30 days,
  about 10 terminal half-lives) and evaluate the final 24 h window, so the
  loading dose is irrelevant to the steady-state comparison.
* **Infusion duration fixed at 1 h.** The source reports infusions "about 1-2 h"
  (Methods 2.2) without specifying what the simulations used; 1 h is used
  throughout. Infusion duration affects Cmax and Cmin but not AUC, so none of the
  AUC-based targets are sensitive to this choice.
* **Absolute PTA values are slightly more permissive than the published
  figures.** Both PTA endpoints are threshold-crossing probabilities on the AUC
  distribution, which under this model is driven entirely by the CL
  distribution. The log-normal approximation has thinner tails than the NPAG
  support points, so the simulated PTA clears 0.9 one dose step lower than the
  paper in the ultra-low and normal albumin strata (see the Figure 3 claims
  table). Every published claim reproduces in direction, and the regimen ordering,
  the albumin gradient, and the identification of 100 mg q12h as the
  toxicity-limiting regimen for normal-albumin patients all hold. Nothing was
  tuned to close the residual gap.
* **Total, not unbound, concentrations.** The model describes total polymyxin B
  (sum of B1, B2, and B1-Ile). The efficacy target was converted from an unbound
  target by the paper using a fixed unbound fraction of 0.42, and the Discussion
  notes that the effect of albumin on the *unbound* concentration could not be
  determined from these data.
* **No CRRT or ECMO patients.** The model was fitted after excluding all patients
  on continuous renal replacement therapy or extracorporeal membrane oxygenation
  and must not be extrapolated to them (Discussion, study limitations).
