---
title: "Ceftriaxone (Tsai 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ceftriaxone (Tsai 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("Tsai_2023_ceftriaxone")()$reference`
* Article: <https://doi.org/10.1093/jac/dkad190>
* Supplement (Tables S1-S3, Figure S1): <https://doi.org/10.1093/jac/dkad190> (JAC Online supplementary data)

Tsai and colleagues studied a novel 2 g three-times-weekly **post-dialysis**
ceftriaxone regimen in Indigenous Australian adults with end-stage renal
disease. The regimen is designed to align every ceftriaxone dose with an
existing dialysis session so that no additional venous cannulation is needed --
vein preservation matters greatly in a population that may require hemodialysis
for decades.

The model has three features that make it unusual for a beta-lactam popPK
model, and all three are reproduced here:

1. **Explicit albumin binding.** Rather than assuming a fixed unbound fraction,
   the model carries protein-bound ceftriaxone as its own state and exchanges it
   with unbound drug through second-order association (`k1`) and first-order
   dissociation (`k2`) rate constants. Both total and unbound plasma
   concentrations are model outputs. This matters because the measured unbound
   fraction in this cohort (median 0.29) was far higher than in healthy
   individuals (0.04-0.17).
2. **Dialysis replaces, rather than augments, clearance.** The Table S2 model
   file states `IF (HDx.EQ.1) CL = CL_HD`, so while a session runs the
   interdialytic clearance arm is switched out entirely for a > 10-fold larger
   dialytic clearance. This differs from the additive dialysis-arm convention
   used elsewhere in `nlmixr2lib` (`Veinstein_2013_gentamicin`,
   `Eyler_2014_ertapenem`, `Jacobs_2016_colistin`).
3. **Bilirubin, not creatinine, drives clearance.** In anuric patients the
   biliary route dominates, so clearance falls with serum bilirubin through an
   inverse-power relationship.

# Population

Sixteen Indigenous Australian adults (13 female, 81%) on three-times-weekly
intermittent high-flux hemodialysis contributed 122 plasma samples assayed for
both total and unbound ceftriaxone (Tsai 2023 Table 1). Median age was 57 years
(IQR 51-64) and median weight 71 kg (IQR 59-83). Median serum albumin was
36 g/L (IQR 33-39) and median total bilirubin 10 umol/L (IQR 6-14); seven
subjects had bilirubin above 10 umol/L, all of acute cause (two cholecystitis,
five severe infection). Infection sources were respiratory (11), urinary (1),
bacteremia (1), skin and soft tissue (1) and intra-abdominal (1). Dialyzers were
high-flux throughout: FX80 (3), FX100 (10), FX120 (3).

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

```{r population}
pop <- readModelDb("Tsai_2023_ceftriaxone")()$population
str(pop[c("species", "n_subjects", "n_samples", "age_median", "weight_median")])
```

# Source trace

Every `ini()` entry carries an in-file comment naming its source location. The
table below collects them for review. "Table S2" refers to the Pmetrics model
file reproduced in the supplementary data.

| Equation / parameter | Value | Source location |
|---|---|---|
| `lcl` (CLnHD) | 0.83 L/h | Table 2, Mean column (SD 0.28, CV 33.33%, median 0.77) |
| `lcl_hemodialysis` (CLHD) | 8.76 L/h | Table 2, Mean column (SD 1.85, CV 21.10%, median 9.13) |
| `lvc` (Vc) | 2.25 L | Table 2, Mean column (SD 1.71, CV 75.81%, median 1.91) |
| `lk1` (Kon) | 1.18 L/mg/h | Table 2, Mean column (SD 0.33, CV 27.76%, median 1.18) |
| `lk2` (Koff) | 206.37 1/h | Table 2, Mean column (SD 48.38, CV 23.44%, median 202.36) |
| `lk12` (Kcp) | 15.37 1/h | Table 2, Mean column (SD 7.24, CV 47.10%, median 17.01) |
| `lk21` (Kpc) | 0.78 1/h | Table 2, Mean column (SD 0.50, CV 64.05%, median 0.55) |
| `e_tbili_cl` | 0.5 (fixed) | Table S2 secondary variables: `CL=CL_nHD*(14.1/Bili)**0.5`; same form in Results |
| bilirubin reference 14.1 umol/L | constant | Table S2 secondary variables (hard-coded; derivation not stated) |
| `bmax_per_g_alb` = 16.7 mg/g | constant | Table S2 secondary variables: `Bmax1=Alb*Vc*16.7` |
| All seven `eta` variances | omega^2 = log(CV^2+1) | Table 2, CV% column (see Assumptions) |
| `addSd`, `addSd_Cunbound` | 0.3 mg/L | Table S2 `#Error` block, C0 for both outputs |
| `propSd`, `propSd_Cunbound` | 0.1 | Table S2 `#Error` block, C1 for both outputs |
| `d/dt(central)`, `d/dt(complex)`, `d/dt(peripheral1)` | n/a | Table S2 `#Differential equations` (XP(1), XP(2), XP(3)) |
| `Cunbound`, `Cc` | n/a | Table S2 `#Output equations` (Y(1), Y(2)) |
| Dialysis switch | n/a | Table S2 secondary variables: `&IF (HDx.EQ.1) CL=CL_HD` |
| Compartment diagram | n/a | Figure S1 |

The `16.7 mg` ceftriaxone bound per gram of albumin encodes **2 binding sites
per albumin molecule**, which is a useful independent check on the transcription:

```{r stoichiometry}
# 2 sites * MW(ceftriaxone) / MW(albumin) * 1000 mg/g
2 * 554.58 / 66437 * 1000
```

# Model structure

```{r model}
mod <- readModelDb("Tsai_2023_ceftriaxone")
mod
```

# Virtual cohort

Original observed data are not publicly available. The cohort below reproduces
the Table 1 covariate distributions: serum albumin centred on the observed
median of 36 g/L and total bilirubin log-normally distributed with median
10 umol/L and an IQR close to the observed 6-14 umol/L.

```{r cohort}
set.seed(20230801)
n_sub <- 200L  # 200 per arm is the vignette cap

cohort <- tibble(
  id    = seq_len(n_sub),
  TBILI = pmax(1, rlnorm(n_sub, meanlog = log(10), sdlog = 0.55)),
  ALB   = pmax(20, rnorm(n_sub, mean = 36, sd = 5))
)

quantile(cohort$TBILI, c(0.25, 0.5, 0.75))  # paper: 10 (IQR 6-14) umol/L
quantile(cohort$ALB,   c(0.25, 0.5, 0.75))  # paper: 36 (IQR 33-39) g/L
```

A helper builds event tables. Two points matter:

* This model declares **two** endpoints (`Cc` and `Cunbound`), each with its own
  residual-error term, and **neither is an ODE state** -- both are algebraic
  observables built from `central`, `complex` and `peripheral1`. rxode2 injects a
  compartment slot for each endpoint after the ODE states and then requires the
  `dvid` -> `cmt` map to be satisfied, so observation rows carry
  `cmt = NA_character_` plus an explicit `dvid = 1L`. Pointing `cmt` at an ODE
  state (`cmt = "central"`) fails with
  `'dvid'->'cmt' ... undefined compartment`. `dvid` is set on the dose rows too so
  the column is not `NA`-typed. Both observables come back as columns regardless
  of which endpoint `dvid` names.
* `RRT_HEMODIAL_ACTIVE` is genuinely time-varying, so the event table is built
  as a plain data frame and the covariate column is set per row. Assigning
  covariates onto an `rxEt` object instead would silently drop them.

```{r events-helper}
# Dialysis sessions are 4.5 h and each dose is given immediately after one.
SESSION_H <- 4.5

make_events <- function(subjects, dose_mg, dose_times, obs_times,
                        sessions = NULL, id_offset = 0L, label = NA_character_) {
  purrr_rows <- lapply(seq_len(nrow(subjects)), function(i) {
    sid <- id_offset + subjects$id[i]
    dose_rows <- data.frame(
      id = sid, time = dose_times, amt = dose_mg, evid = 1L,
      cmt = "central", dvid = 1L
    )
    obs_rows <- data.frame(
      id = sid, time = obs_times, amt = 0, evid = 0L,
      cmt = NA_character_, dvid = 1L
    )
    ev <- dplyr::arrange(dplyr::bind_rows(dose_rows, obs_rows), time, dplyr::desc(evid))
    ev$TBILI     <- subjects$TBILI[i]
    ev$ALB       <- subjects$ALB[i]
    ev$treatment <- label
    on <- rep(FALSE, nrow(ev))
    for (s in sessions) on <- on | (ev$time >= s[1] & ev$time < s[2])
    ev$RRT_HEMODIAL_ACTIVE <- as.numeric(on)
    ev
  })
  dplyr::bind_rows(purrr_rows)
}
```

# Simulation

The primary scenario is one 2 g post-dialysis dose followed by a 72 h
inter-dialysis interval -- the window the paper uses for its target-attainment
analysis. No dialysis runs during this window, because the dose is given at the
end of a session.

```{r simulate}
obs_grid <- sort(unique(c(seq(0, 72, by = 0.5), seq(0, 2, by = 0.1))))

ev_72 <- make_events(cohort, dose_mg = 2000, dose_times = 0,
                     obs_times = obs_grid, label = "2 g post-dialysis")
stopifnot(!anyDuplicated(unique(ev_72[, c("id", "time", "evid")])))

sim_72 <- rxode2::rxSolve(
  mod, events = ev_72, sigma = NA,
  keep = c("TBILI", "ALB", "treatment", "RRT_HEMODIAL_ACTIVE")
) |>
  as.data.frame()

nrow(sim_72)
```

Typical-value (no between-subject variability) profiles use `omega = NA` rather
than `zeroRe()`, because `zeroRe()` mutates shared model state. A single-subject
cohort at the cohort-median covariates serves for the figure replications below.

```{r simulate-typical}
typ <- tibble(id = 1L, TBILI = 10, ALB = 36)  # cohort medians
```

# Replicate published figures

## Figure 1 -- total and unbound concentration-time profiles

Figure 1 of Tsai 2023 shows observed total (grey) and unbound (black)
ceftriaxone with the final model's predicted lines for each patient. Here the
typical-value profile is shown across a full weekly cycle (doses on Day 1, Day 4
and Day 6, each immediately after a 4.5 h dialysis session), which exposes both
the slow interdialytic decline and the sharp intradialytic drops.

```{r figure-1, fig.cap = "Replicates the structure of Figure 1 of Tsai 2023: typical-value total and unbound ceftriaxone across one weekly cycle of the 2 g three-times-weekly post-dialysis regimen. Shaded bands are dialysis sessions."}
dose_t <- c(0, 72, 120)                                    # Day 1, Day 4, Day 6
sess   <- lapply(c(dose_t, 168), function(d) c(d - SESSION_H, d))
sess   <- Filter(function(s) s[2] > 0, sess)

week_obs <- sort(unique(c(seq(0, 168, by = 0.25), unlist(sess), dose_t)))
ev_week  <- make_events(typ, dose_mg = 2000, dose_times = dose_t,
                        obs_times = week_obs, sessions = sess, label = "typical")
sim_week <- rxode2::rxSolve(mod, events = ev_week, omega = NA, sigma = NA) |>
  as.data.frame()

sess_df <- do.call(rbind, lapply(sess, function(s) data.frame(xmin = max(0, s[1]), xmax = s[2])))

sim_week |>
  select(time, Total = Cc, Unbound = Cunbound) |>
  pivot_longer(-time, names_to = "Analyte", values_to = "conc") |>
  ggplot(aes(time, conc, colour = Analyte)) +
  geom_rect(data = sess_df, inherit.aes = FALSE,
            aes(xmin = xmin, xmax = xmax, ymin = 0.1, ymax = Inf),
            fill = "grey70", alpha = 0.4) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  scale_colour_manual(values = c(Total = "grey40", Unbound = "black")) +
  labs(x = "Time (h)", y = "Ceftriaxone concentration (mg/L)",
       title = "Weekly cycle, 2 g three-times-weekly post-dialysis",
       caption = "Replicates the structure of Figure 1 of Tsai 2023.")
```

## Figure 2 -- inverse-power clearance versus bilirubin

```{r figure-2, fig.cap = "Replicates Figure 2 of Tsai 2023: the inverse-power relationship between interdialytic ceftriaxone clearance and serum bilirubin."}
tibble(TBILI = seq(2, 90, by = 0.5)) |>
  mutate(CL = 0.83 * (14.1 / TBILI)^0.5) |>
  ggplot(aes(TBILI, CL)) +
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = 14.1, linetype = "dashed", colour = "grey50") +
  annotate("text", x = 14.1, y = 1.9, label = "reference 14.1 umol/L",
           hjust = -0.05, size = 3, colour = "grey30") +
  labs(x = "Serum bilirubin (umol/L)", y = "Interdialytic clearance CL (L/h)",
       title = "CL = CLnHD * (14.1 / bilirubin)^0.5",
       caption = "Replicates Figure 2 of Tsai 2023 (reported r^2 = 0.74).")
```

# Protein binding

The paper reports a median unbound fraction of 0.29 (IQR 0.20-0.40) -- much
higher than the 0.04-0.17 typical of healthy individuals. The explicit binding
model reproduces this without any fitted unbound-fraction parameter.

```{r fu, fig.cap = "Simulated unbound fraction over the interdialytic interval versus the values reported by Tsai 2023."}
fu_df <- sim_72 |>
  filter(time >= 12, Cc > 0) |>
  mutate(fu = Cunbound / Cc)

fu_summary <- tibble(
  Source = c("Simulated", "Tsai 2023 (observed)"),
  Median = c(median(fu_df$fu), 0.29),
  Q1     = c(quantile(fu_df$fu, 0.25), 0.20),
  Q3     = c(quantile(fu_df$fu, 0.75), 0.40)
)

fu_summary |>
  rename("Unbound fraction" = Source, "Q1" = Q1, "Median" = Median, "Q3" = Q3) |>
  knitr::kable(digits = 3, caption = "Unbound fraction: simulated vs Tsai 2023 Results.")

ggplot(fu_df, aes(fu)) +
  geom_histogram(bins = 50, fill = "grey60", colour = "white") +
  geom_vline(xintercept = c(0.20, 0.29, 0.40), linetype = c(3, 1, 3)) +
  labs(x = "Unbound fraction", y = "Count",
       title = "Simulated unbound fraction",
       caption = "Solid line = Tsai 2023 median 0.29; dotted = reported IQR 0.20-0.40.")
```

# Dialysis effect

Two reported quantities probe the dialysis switch: clearance is > 10-fold higher
during a session, and unbound concentrations fall by a median 70% (IQR 64-74%)
across each session.

```{r dialysis}
# Clearance ratio at the reference bilirubin.
c(CLnHD = 0.83, CLHD = 8.76, ratio = 8.76 / 0.83)

# 4.5 h session beginning 48 h after a 2 g dose.
ev_hd <- make_events(cohort, dose_mg = 2000, dose_times = 0,
                     obs_times = sort(unique(c(seq(0, 72, by = 0.5), 48, 48 + SESSION_H))),
                     sessions = list(c(48, 48 + SESSION_H)), label = "with dialysis")
sim_hd <- rxode2::rxSolve(mod, events = ev_hd, sigma = NA, keep = "treatment") |>
  as.data.frame()

reduction <- sim_hd |>
  filter(time %in% c(48, 48 + SESSION_H)) |>
  group_by(id) |>
  summarise(pct = 100 * (1 - Cunbound[time == 48 + SESSION_H] / Cunbound[time == 48]),
            .groups = "drop")

tibble(
  Quantity = c("Simulated", "Tsai 2023 (observed)"),
  Median   = c(median(reduction$pct), 70),
  Q1       = c(quantile(reduction$pct, 0.25), 64),
  Q3       = c(quantile(reduction$pct, 0.75), 74)
) |>
  rename("Unbound reduction per session (%)" = Quantity) |>
  knitr::kable(digits = 1,
               caption = "Unbound concentration reduction across one 4.5 h dialysis session.")
```

The simulated median reduction is smaller than the reported 70%; the cause is
worked through in *Assumptions and deviations* below (it stems from the paper's
own half-lives being one-compartment hand calculations, not from a transcription
difference).

# PKNCA validation

NCA is run on the 72 h interdialytic interval, separately for total and unbound
concentrations, stratified by treatment.

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

# Guarantee a time-zero record per subject so PKNCA can anchor AUC.
sim_nca <- bind_rows(
  sim_nca,
  sim_nca |> distinct(id, treatment) |> mutate(time = 0, Cc = 0, Cunbound = 0)
) |>
  distinct(id, treatment, time, .keep_all = TRUE) |>
  arrange(id, treatment, time)

dose_df <- ev_72 |>
  filter(evid == 1) |>
  select(id, time, amt, treatment)

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, auclast = TRUE, half.life = TRUE
)

nca_total <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id),
  PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id),
  intervals = intervals
))

nca_unbound <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(sim_nca, Cunbound ~ time | treatment + id),
  PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id),
  intervals = intervals
))
```

```{r pknca-summary}
summarise_nca <- function(res, analyte) {
  as.data.frame(res) |>
    filter(!is.na(PPORRES), start == 0, end == Inf) |>
    group_by(PPTESTCD) |>
    summarise(Median = median(PPORRES), Q1 = quantile(PPORRES, 0.25),
              Q3 = quantile(PPORRES, 0.75), .groups = "drop") |>
    mutate(Analyte = analyte, .before = 1)
}

bind_rows(summarise_nca(nca_total, "Total"),
          summarise_nca(nca_unbound, "Unbound")) |>
  rename("NCA parameter" = PPTESTCD) |>
  knitr::kable(digits = 2,
               caption = "Simulated NCA over the 72 h interdialytic interval (2 g dose).")
```

## Comparison against published values

The paper reports no conventional NCA table (no Cmax or AUC), so the only
directly comparable NCA parameter is the interdialytic half-life. It is compared
here with `ncaComparisonTable()`; the remaining reported quantities -- troughs,
unbound fraction and dialysis reduction -- are compared in the table that
follows.

```{r nca_compare}
published <- tibble::tibble(
  treatment = "2 g post-dialysis",
  half.life = 38.31   # Tsai 2023 Table 2, T(1/2)nHD mean (median 30.40)
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_unbound,
  reference     = published,
  by            = "treatment",
  params        = "half.life",
  units         = c(half.life = "h"),
  tolerance_pct = 20
)

knitr::kable(cmp, caption = paste(
  "Simulated vs published interdialytic half-life.",
  "* marks a difference of more than 20% from the reference."
))
```

```{r other-compare}
trough_48 <- sim_72 |> filter(time == 48) |> pull(Cunbound)
trough_72 <- sim_72 |> filter(time == 72) |> pull(Cunbound)

tibble::tribble(
  ~Quantity,                                        ~Simulated,                ~`Tsai 2023 reported`,
  "Unbound trough, 48 h interval (mg/L)",           median(trough_48),         18.2,
  "Unbound trough, 72 h interval (mg/L)",           median(trough_72),         8.8,
  "Unbound fraction",                               median(fu_df$fu),          0.29,
  "Unbound reduction per dialysis session (%)",     median(reduction$pct),     70,
  "Interdialytic clearance at bilirubin 14.1 (L/h)", 0.83,                     0.83,
  "Dialytic / interdialytic clearance ratio",        8.76 / 0.83,              10
) |>
  knitr::kable(digits = 2, caption = paste(
    "Simulated medians versus values reported in Tsai 2023 Results.",
    "Reported IQRs: 48 h trough 9.7-25.9; 72 h trough 7.1-17.7;",
    "unbound fraction 0.20-0.40; session reduction 64-74%."
  ))
```

Both simulated troughs fall inside the reported interquartile ranges, and the
simulated unbound fraction sits just below the reported median but inside its
IQR. The half-life agrees closely with the reported mean.

# Target attainment (Table 3)

The paper's target is 100% *f*T > MIC over the final 24 h of a 72 h
inter-dialysis interval. Computing the minimum unbound concentration in that
window once per subject yields the whole MIC grid from a single simulation per
scenario.

```{r pta}
pta_scenarios <- tidyr::expand_grid(dose_mg = c(1000, 2000), bili = c(5, 10, 20))

pta_min <- function(dose_mg, bili, idx) {
  subj <- cohort |> mutate(TBILI = bili)
  ev <- make_events(subj, dose_mg = dose_mg, dose_times = 0,
                    obs_times = sort(unique(c(seq(0, 48, by = 2), seq(48, 72, by = 0.5)))),
                    id_offset = as.integer(idx * 1000L),
                    label = paste0(dose_mg / 1000, " g, bili ", bili))
  rxode2::rxSolve(mod, events = ev, sigma = NA, keep = "treatment") |>
    as.data.frame() |>
    filter(time >= 48, time <= 72) |>
    group_by(id, treatment) |>
    summarise(min_fC = min(Cunbound), .groups = "drop")
}

pta_raw <- bind_rows(lapply(seq_len(nrow(pta_scenarios)), function(i) {
  pta_min(pta_scenarios$dose_mg[i], pta_scenarios$bili[i], i) |>
    mutate(dose_mg = pta_scenarios$dose_mg[i], bili = pta_scenarios$bili[i])
}))

mic_grid <- c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16)

pta_tbl <- tidyr::expand_grid(pta_raw, MIC = mic_grid) |>
  group_by(dose_mg, bili, MIC) |>
  summarise(PTA = 100 * mean(min_fC > MIC), .groups = "drop") |>
  mutate(Regimen = paste0(dose_mg / 1000, " g")) |>
  select(Regimen, `Bili (umol/L)` = bili, MIC, PTA) |>
  pivot_wider(names_from = MIC, values_from = PTA)

knitr::kable(pta_tbl, digits = 1, caption = paste(
  "Simulated PTA (%) for 100% fT > MIC over the final 24 h of a 72 h interval.",
  "Compare Tsai 2023 Table 3; see Assumptions for the systematic offset."
))
```

Published values for the same cells (Tsai 2023 Table 3, MIC = 1 mg/L): 2 g at
bilirubin 5 / 10 / 20 = 97.9 / 100 / 100%; 1 g at bilirubin 5 / 10 / 20 =
86 / 99.5 / 100%. The simulation reproduces the qualitative findings the paper's
recommendations rest on -- PTA rises with bilirubin (because clearance falls),
2 g clears the 95% bar where 1 g does not at low bilirubin -- but is
systematically several points low. The cause is diagnosed below.

# Assumptions and deviations

* **Typical values are Table 2 means.** Table 2 reports both a mean and a median
  for every parameter; the mean column is used here, with the median recorded in
  each `ini()` comment. Pmetrics NPAG estimates a discrete non-parametric
  distribution, so neither column is a "typical value" in the parametric sense.

* **IIV is a log-normal approximation of a non-parametric distribution.** NPAG
  produces a set of joint support points, not an omega matrix. Table 2 summarises
  that distribution only by its per-parameter mean, SD and CV%, so the CV% is
  carried into independent log-normal random effects via
  `omega^2 = log(CV^2 + 1)`. **The joint structure is unrecoverable from the
  published summary**, and this is the single largest deviation in this
  extraction. It has a measurable consequence: independent sampling generates
  parameter combinations that the true support-point distribution would not
  contain. About 5% of simulated subjects pair a high clearance with a low
  central volume (median CL 1.90 vs 1.31 L/h and Vc 0.86 vs 2.37 L in that
  subgroup), giving an implausibly fast elimination rate and a near-zero 72 h
  concentration. Those subjects fail the PTA target at *every* MIC, which places
  a ceiling of roughly 93% on simulated PTA and explains essentially all of the
  offset from Table 3. The published PTA values should be preferred over the
  simulated ones; the simulation is included to show the model reproduces the
  qualitative dose- and bilirubin-dependence, not to restate Table 3.

* **Residual error is the assay polynomial, unscaled.** The Table S2 `#Error`
  block gives an assay-error polynomial of `0.3, 0.1, 0, 0` for each output
  (SD = 0.3 + 0.1 * concentration), which is encoded here as a 0.3 mg/L additive
  plus 10% proportional term on both `Cc` and `Cunbound`. Pmetrics multiplies
  that polynomial by an estimated noise-inflation factor gamma; the model file
  sets only the gamma *starting* value (`G=2`) and the final estimated gamma is
  not reported anywhere in the paper or supplement. The polynomial is therefore
  carried unscaled (equivalent to gamma = 1), the minimum-assumption reading. If
  the final gamma were near its starting value, the true residual SD would be
  about twice what is encoded here.

* **Dialysis session duration is assumed to be 4.5 h.** The paper does not state
  it. 4.5 h is a standard prescription and is used for every figure and table
  above; no reported quantity was used to select it.

* **The reported per-session unbound reduction is not reproducible from Table 2's
  parameters, and this is a property of the paper rather than of the
  transcription.** Table 2 footnote *a* marks T(1/2)HD and T(1/2)nHD as
  "manually calculated", and they are consistent with the one-compartment
  formula `0.693 * Vd / CL`: `0.693 * 36.22 / 8.76 = 2.87 h` against a reported
  T(1/2)HD of 3.04 h (median 2.65). A 4.5 h session at a 2.65 h half-life gives
  a 69% drop, matching the reported 70%. But the *model* in Table S2 is a
  three-state system in which elimination acts only on unbound drug in the
  central compartment, while `Kcp / Kpc = 15.37 / 0.78 = 19.7` puts about 95% of
  the drug in the peripheral compartment. During a short session the return flux
  from the periphery (`Kpc = 0.78 1/h`) rate-limits the decline, so the model's
  true intradialytic half-life is about 5 h and a 4.5 h session removes roughly
  half the unbound drug rather than 70%. The one-compartment hand calculation
  therefore understates the model's intradialytic half-life. Nothing was tuned
  to close this gap. The same arithmetic explains why Table 2's Vd of 36.22 L
  differs from `Vc * (1 + Kcp/Kpc) = 46.6 L`: the derived rows were computed per
  support point with a reduced formula, not from the three-state model.

* **The "administered during dialysis" counterfactual is not reproduced.** The
  paper reports a PTA of 0.1% when the 2 g dose is given as a 30 min infusion
  during the last 30 min of a session. Simulating exactly that -- 30 min of
  dialysis running from the start of the dose, then the session ending -- gives a
  far higher PTA, because fast albumin binding (`k2 = 206 1/h`) and fast
  distribution (`k12 = 15.4 1/h`) move most of the dose out of the central
  compartment within minutes, where dialysis cannot reach it. Reproducing 0.1%
  would require the session to continue for several hours after the dose. Since
  the supplement does not show how `HDx` was set in that simulation, the scenario
  is omitted rather than guessed at. The paper's clinical recommendation against
  dosing during dialysis is not in question -- the model does predict substantial
  intradialytic loss -- only the specific 0.1% figure is unreproducible from the
  published information.

* **Covariates are time-fixed.** Serum bilirubin and albumin were treated as
  time-fixed per subject, as in the source analysis. `TBILI` enters as a
  denominator and must be strictly positive; the virtual cohort floors it at
  1 umol/L.

* **Bilirubin reference value.** The 14.1 umol/L reference in the clearance
  covariate is hard-coded in the Table S2 model file and the paper does not state
  its derivation. It is not the cohort median (10 umol/L). It is carried exactly
  as printed.

* **Dose administration.** Doses are simulated as bolus inputs into `central`.
  The paper gave each 2 g dose by slow injection within 5 min, which is
  negligible against a 48-72 h dosing interval. The Pmetrics model file uses
  `RATEIV(1)`, i.e. an infusion into compartment 1; infusion versus bolus is an
  event-table choice and requires no model change.
