---
title: "Piperacillin + tazobactam in infants (Cohen-Wolkowiez 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Piperacillin + tazobactam in infants (Cohen-Wolkowiez 2014)}
  %\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

Cohen-Wolkowiez et al. 2014 developed independent one-compartment population PK
models for piperacillin and tazobactam (coadministered as the 8:1 fixed-ratio
combination piperacillin-tazobactam) in 32 critically ill premature and term
infants under 61 days postnatal age (PNA). Both drugs share the same study
population, the same dosing regimen, and the same plasma sampling scheme; the
two structural models differ in their covariate composition on clearance. This
vignette walks the paper's narrative for both drugs and uses each `modellib()`
entry at the point in the narrative where the corresponding model is
introduced.

* Citation (piperacillin): `r readModelDb("CohenWolkowiez_2014_piperacillin")()$reference`
* Citation (tazobactam):   `r readModelDb("CohenWolkowiez_2014_tazobactam")()$reference`
* Article: <https://doi.org/10.1128/AAC.02139-13>

# Population

Thirty-two infants were enrolled 2010-2011 at four US neonatal intensive care
units in an open-label, prospective, early-phase PK and safety study of
intravenous piperacillin-tazobactam in young infants under 61 days postnatal
age with suspected systemic infection. Enrollment was stratified by gestational
age (GA) at birth and postnatal age (PNA) into four cohorts: cohort 1 (n = 12,
GA < 32 wk, PNA < 14 d), cohort 2 (n = 9, GA < 32 wk, PNA >= 14 d), cohort 3
(n = 8, GA >= 32 wk, PNA < 14 d), and cohort 4 (n = 3, GA >= 32 wk, PNA >= 14 d).
Cohorts 1-3 received piperacillin 80 mg/kg every 8 h IV (30-min infusion);
cohort 4 received 100 mg/kg every 8 h. Demographics across the pooled cohort
(Cohen-Wolkowiez 2014 Table 1): median GA 30 weeks (range 23-40), median PNA
8 days (1-60), median postmenstrual age (PMA) 32 weeks (25-48), median weight
1.439 kg (0.473-3.990), median serum creatinine 0.8 mg/dL (0.3-2.0). 63% of
infants received concurrent gentamicin, 13% dopamine, 6% epinephrine; 63% were
male. 128 timed plasma piperacillin and tazobactam samples (median 4 samples
per infant) were retained for the population PK analysis. Concurrent
dried-blood-spot samples were also collected; the structural PK models
extracted here were fit to plasma alone.

The same metadata is available programmatically via
`readModelDb("CohenWolkowiez_2014_piperacillin")()$population` (and the
equivalent for tazobactam).

# Source trace

The per-parameter origin is recorded inline next to each `ini()` entry in
`inst/modeldb/specificDrugs/CohenWolkowiez_2014_piperacillin.R` and
`inst/modeldb/specificDrugs/CohenWolkowiez_2014_tazobactam.R`. The table below
collects them in one place for review.

## Piperacillin (Cohen-Wolkowiez 2014 Table 3 final PMA-based model)

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (typical CL per kg body weight, L/kg/h, at PMA = 33 wk) | `log(0.080)` | Table 3: typical CL value 0.080 L/kg/h (RSE 7.9%) |
| `lvc` (typical V per kg body weight, L/kg) | `log(0.42)` | Table 3: typical V value 0.42 L/kg (RSE 9.6%) |
| `e_wt_cl` (linear WT exponent on CL, unitless) | `fixed(1.0)` | Population PK model building paragraph: WT exponent fixed at 1; free body-size and 3/4-power scaling explored and rejected |
| `e_wt_vc` (linear WT exponent on V, unitless) | `fixed(1.0)` | Population PK model building paragraph: WT exponent fixed at 1 |
| `e_page_cl` (power exponent on (PMA / 33 wk) for CL) | `1.76` | Table 3: PMA covariate 1.76 (RSE 33.6%) |
| `etalcl` (omega_CL on internal scale) | `0.12839` (37% CV) | Table 3: CL interindividual variability 37% CV (RSE 27.5%); omega^2 = log(1 + CV^2) |
| (no `etalvc`) | n/a | Population PK model building paragraph: between-subject variability on V fixed to 0 (boundary value) |
| `propSd` (proportional residual error) | `0.33` | Table 3: proportional residual error 33% CV (RSE 9.9%) |
| `addSd` (additive residual error, mg/L) | `6.90` | Table 3: plasma additive residual error 6.90 mg/L (RSE 42.6%) |
| `d/dt(central)` (one-compartment IV PK) | n/a | Population PK model building paragraph: one-compartment structural model selected |
| Piperacillin CL covariate equation `CL = theta_CL * WT * (PMA / 33)^COVPMA` | n/a | Table 2 final piperacillin model row and Population PK model building paragraph |

## Tazobactam (Cohen-Wolkowiez 2014 Table 4 final model)

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (typical CL per kg body weight, L/kg/h, at PMA = 33 wk, SCR = 0.5 mg/dL, no concurrent gentamicin) | `log(0.088)` | Table 4: typical CL value 0.088 L/kg/h (RSE 8.9%) |
| `lvc` (typical V per kg body weight, L/kg) | `log(0.57)` | Table 4: typical V value 0.57 L/kg (RSE 5.4%) |
| `e_wt_cl` (linear WT exponent on CL, unitless) | `fixed(1.0)` | Population PK model building paragraph: WT exponent fixed at 1 |
| `e_wt_vc` (linear WT exponent on V, unitless) | `fixed(1.0)` | Population PK model building paragraph: WT exponent fixed at 1 |
| `e_page_cl` (power exponent on (PMA / 33 wk) for CL) | `1.35` | Table 4: PMA covariate 1.35 (RSE 27.3%) |
| `e_creat_cl` (power exponent on (CREAT / 0.5 mg/dL) for CL) | `-0.35` | Table 4: SCR covariate -0.35 (RSE 37.0%) |
| `e_conmed_gen_cl` (multiplicative factor on CL for GENT = 1) | `1.52` | Table 4: GENT covariate 1.52 (RSE 10.1%) |
| `etalcl` (omega_CL on internal scale) | `0.05154` (23% CV) | Table 4: CL interindividual variability 23% CV (RSE 34.0%); omega^2 = log(1 + CV^2) |
| (no `etalvc`) | n/a | Population PK model building paragraph: between-subject variability on V fixed to 0 (boundary value) |
| `propSd` (proportional residual error) | `0.24` | Table 4: proportional residual error 24% CV (RSE 19.8%) |
| `addSd` (additive residual error, mg/L) | `1.43` | Table 4: plasma additive residual error 1.43 mg/L (RSE 51.2%) |
| `d/dt(central)` (one-compartment IV PK) | n/a | Population PK model building paragraph: one-compartment structural model selected |
| Tazobactam CL covariate equation `CL = theta_CL * WT * (PMA / 33)^COVPMA * (SCR / 0.5)^COVSCR * COVGENT^GENT` | n/a | Table 2 final tazobactam model row and Population PK model building paragraph |

# Virtual cohort

Original observed concentrations are not publicly available. The simulations
below use a virtual cohort whose PMA / WT / SCR distributions mirror the
Cohen-Wolkowiez 2014 Table 1 demographics and whose dosing follows the per-
cohort protocol assignments (80 mg/kg q8h piperacillin for cohorts 1-3; 100
mg/kg q8h for cohort 4). The tazobactam component is dosed at one-eighth the
piperacillin dose (fixed 8:1 piperacillin:tazobactam ratio).

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

T_INF         <- 0.5   # hour; 30-min IV infusion
DOSE_INTERVAL <- 8     # hours
N_DOSES       <- 5     # five doses to approach steady state

# Per-cohort typical demographics from Cohen-Wolkowiez 2014 Table 1
# (median values; PMA in weeks, WT in kg, SCR in mg/dL, dose in mg/kg).
cohorts <- tibble::tribble(
  ~cohort,            ~n,  ~PMA_wk, ~WT_kg, ~CREAT, ~CONMED_GEN, ~dose_pip_mg_per_kg,
  "cohort 1",         50,  29,      1.109,  0.8,    0,           80,
  "cohort 2",         50,  33,      1.397,  0.8,    1,           80,
  "cohort 3",         50,  39,      3.248,  0.8,    0,           80,
  "cohort 4",         50,  43,      3.050,  0.6,    0,           100
)

# Helper to build one cohort's event table. id_offset shifts subject IDs so
# bind_rows() across cohorts does not collide on `id` (rxode2 silently merges
# duplicate IDs across cohorts).
make_cohort <- function(spec, id_offset, dose_mg_per_kg, drug,
                        ratio_piperacillin_to_tazobactam = 8L) {
  WT_kg <- spec$WT_kg
  per_dose_mg <- if (drug == "piperacillin") {
    dose_mg_per_kg * WT_kg
  } else {
    dose_mg_per_kg * WT_kg / ratio_piperacillin_to_tazobactam
  }
  rows <- vector("list", spec$n)
  for (k in seq_len(spec$n)) {
    sid <- id_offset + k
    ev <- rxode2::et(
      amt   = per_dose_mg,
      rate  = per_dose_mg / T_INF,
      ii    = DOSE_INTERVAL,
      addl  = N_DOSES - 1L,
      cmt   = "central"
    ) |>
      rxode2::et(seq(0, N_DOSES * DOSE_INTERVAL, by = 0.25))
    df <- as.data.frame(ev)
    df$id         <- sid
    df$cohort     <- spec$cohort
    df$WT         <- WT_kg
    df$PAGE       <- spec$PMA_wk / 4.35   # PMA (wk) -> PAGE canonical (months)
    df$CREAT      <- spec$CREAT
    df$CONMED_GEN <- spec$CONMED_GEN
    rows[[k]] <- df
  }
  dplyr::bind_rows(rows)
}

build_events <- function(drug) {
  offsets <- c(0L, cumsum(head(cohorts$n, -1)))
  parts <- lapply(seq_len(nrow(cohorts)), function(i) {
    make_cohort(
      cohorts[i, ],
      id_offset      = offsets[[i]],
      dose_mg_per_kg = cohorts$dose_pip_mg_per_kg[[i]],
      drug           = drug
    )
  })
  ev <- dplyr::bind_rows(parts)
  stopifnot(!anyDuplicated(unique(ev[, c("id", "time", "evid")])))
  ev
}

events_pip <- build_events("piperacillin")
events_taz <- build_events("tazobactam")

# rxode2::et(..., addl = ) encodes the entire repeated-dose schedule as a
# single event row. For PKNCA we need one row per administered dose, so
# expand the schedule analytically using the known cohort + dose layout.
expand_dose_schedule <- function(events, drug,
                                 ratio_piperacillin_to_tazobactam = 8L) {
  one_per_subj <- events |>
    dplyr::filter(!is.na(amt), amt > 0) |>
    dplyr::distinct(id, cohort, .keep_all = TRUE) |>
    dplyr::select(id, cohort)
  dose_times <- seq(0, by = DOSE_INTERVAL, length.out = N_DOSES)
  per_cohort_dose_pip <- setNames(cohorts$dose_pip_mg_per_kg, cohorts$cohort)
  per_cohort_wt       <- setNames(cohorts$WT_kg,              cohorts$cohort)
  tidyr::expand_grid(one_per_subj, time = dose_times) |>
    dplyr::mutate(
      amt = {
        per_kg <- per_cohort_dose_pip[as.character(cohort)]
        if (drug == "tazobactam") per_kg <- per_kg / ratio_piperacillin_to_tazobactam
        per_kg * per_cohort_wt[as.character(cohort)]
      }
    )
}

doses_pip <- expand_dose_schedule(events_pip, "piperacillin")
doses_taz <- expand_dose_schedule(events_taz, "tazobactam")
```

# Simulation

## Piperacillin (one-compartment, PMA-based final model)

```{r simulate-pip}
mod_pip <- readModelDb("CohenWolkowiez_2014_piperacillin")
sim_pip <- rxode2::rxSolve(mod_pip, events = events_pip,
                           keep = c("cohort", "WT", "PAGE"))
sim_pip <- as.data.frame(sim_pip)
```

## Tazobactam (one-compartment with PMA + SCR + GENT covariates)

```{r simulate-taz}
mod_taz <- readModelDb("CohenWolkowiez_2014_tazobactam")
sim_taz <- rxode2::rxSolve(mod_taz, events = events_taz,
                           keep = c("cohort", "WT", "PAGE", "CREAT", "CONMED_GEN"))
sim_taz <- as.data.frame(sim_taz)
```

# Replicate published figures

Cohen-Wolkowiez 2014 Figures 1 and 2 are diagnostic plots from the
NONMEM fit (observed vs individual predictions, VPC of concentration
vs time, CWRES vs population predictions / time). Reproducing the
VPC panel against simulated trajectories shows the per-cohort
piperacillin and tazobactam profiles across one dosing interval at
steady state.

```{r figure-1-vpc-pip}
# Reproduces the simulation side of Cohen-Wolkowiez 2014 Figure 1B
# (piperacillin VPC) at steady state: simulated piperacillin
# concentration vs time within the final dosing interval, with 90% prediction
# interval and median by cohort.
ss_window_pip <- sim_pip |>
  dplyr::filter(time >= (N_DOSES - 1) * DOSE_INTERVAL,
                time <= N_DOSES * DOSE_INTERVAL,
                !is.na(Cc)) |>
  dplyr::mutate(time_in_interval = time - (N_DOSES - 1) * DOSE_INTERVAL)

vpc_pip <- ss_window_pip |>
  dplyr::group_by(cohort, time_in_interval) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_pip, aes(time_in_interval, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~ cohort) +
  scale_y_log10() +
  labs(
    x = "Time within steady-state dosing interval (h)",
    y = "Piperacillin Cc (mg/L)",
    title = "Steady-state piperacillin VPC by cohort",
    caption = "Replicates the simulation side of Cohen-Wolkowiez 2014 Figure 1B."
  )
```

```{r figure-2-vpc-taz}
# Reproduces the simulation side of Cohen-Wolkowiez 2014 Figure 2B
# (tazobactam VPC) at steady state.
ss_window_taz <- sim_taz |>
  dplyr::filter(time >= (N_DOSES - 1) * DOSE_INTERVAL,
                time <= N_DOSES * DOSE_INTERVAL,
                !is.na(Cc)) |>
  dplyr::mutate(time_in_interval = time - (N_DOSES - 1) * DOSE_INTERVAL)

vpc_taz <- ss_window_taz |>
  dplyr::group_by(cohort, time_in_interval) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_taz, aes(time_in_interval, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~ cohort) +
  scale_y_log10() +
  labs(
    x = "Time within steady-state dosing interval (h)",
    y = "Tazobactam Cc (mg/L)",
    title = "Steady-state tazobactam VPC by cohort",
    caption = "Replicates the simulation side of Cohen-Wolkowiez 2014 Figure 2B."
  )
```

```{r figure-3-bayesian-cl-vs-pma}
# Reproduces the typical-value side of Cohen-Wolkowiez 2014 Figure 3:
# population-typical piperacillin clearance (per kg) vs PMA. With
# random effects zeroed, the curve shows the PMA maturation only
# (since linear WT scaling has no effect on per-kg CL).
typical_pip <- rxode2::zeroRe(mod_pip)
pma_grid <- seq(25, 50, by = 1)
typical_cl <- vapply(pma_grid, function(pma_wk) {
  ev1 <- rxode2::et(amt = 80, rate = 80 / T_INF, cmt = "central") |>
    rxode2::et(c(0.5, 1, 8))
  ev1$WT <- 1
  ev1$PAGE <- pma_wk / 4.35
  out <- rxode2::rxSolve(typical_pip, ev1)
  # CL = -V * ln(C(t2)/C(t1)) / (t2 - t1); use post-infusion samples
  out <- as.data.frame(out)
  post <- out[out$time %in% c(1, 8), ]
  V <- exp(log(0.42))   # vc per kg for the typical model
  cl <- -V * (log(post$Cc[post$time == 8]) - log(post$Cc[post$time == 1])) / (8 - 1)
  cl
}, numeric(1))

df_typical <- tibble::tibble(PMA_wk = pma_grid, CL_per_kg = typical_cl)
ggplot(df_typical, aes(PMA_wk, CL_per_kg)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = 33, linetype = "dashed", colour = "grey50") +
  labs(
    x = "Postmenstrual age (weeks)",
    y = "Typical piperacillin CL (L/kg/h)",
    title = "Typical-value piperacillin CL maturation curve vs PMA",
    caption = "Replicates the typical-value side of Cohen-Wolkowiez 2014 Figure 3 (PMA reference 33 weeks, dashed)."
  )
```

# PKNCA validation

Steady-state NCA per cohort. The fixed-effect (typical-value) trajectory is
the appropriate target for comparison against the empirical-Bayesian summary
the source paper reports in Table 5; full-IIV simulations are also computed
to illustrate the spread.

## Piperacillin (steady state, AUC0-tau, Cmax,ss, Cmin,ss)

```{r pknca-pip}
tau <- DOSE_INTERVAL
start_ss <- (N_DOSES - 1) * tau
end_ss   <- N_DOSES * tau

sim_nca_pip <- sim_pip |>
  dplyr::filter(!is.na(Cc), time >= start_ss, time <= end_ss) |>
  dplyr::select(id, time, Cc, cohort)

dose_nca_pip <- doses_pip |>
  dplyr::filter(time >= start_ss, time <= end_ss)

conc_obj_pip <- PKNCA::PKNCAconc(sim_nca_pip, Cc ~ time | cohort + id,
                                 concu = "mg/L", timeu = "hr")
dose_obj_pip <- PKNCA::PKNCAdose(dose_nca_pip, amt ~ time | cohort + id,
                                 doseu = "mg")
intervals_pip <- data.frame(
  start    = start_ss,
  end      = end_ss,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  auclast  = TRUE,
  cav      = TRUE,
  ctrough  = TRUE
)
nca_pip <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj_pip, dose_obj_pip,
                                          intervals = intervals_pip))
knitr::kable(summary(nca_pip),
             caption = "Simulated steady-state piperacillin NCA by cohort.")
```

## Tazobactam (steady state, AUC0-tau, Cmax,ss, Cmin,ss)

```{r pknca-taz}
sim_nca_taz <- sim_taz |>
  dplyr::filter(!is.na(Cc), time >= start_ss, time <= end_ss) |>
  dplyr::select(id, time, Cc, cohort)

dose_nca_taz <- doses_taz |>
  dplyr::filter(time >= start_ss, time <= end_ss)

conc_obj_taz <- PKNCA::PKNCAconc(sim_nca_taz, Cc ~ time | cohort + id,
                                 concu = "mg/L", timeu = "hr")
dose_obj_taz <- PKNCA::PKNCAdose(dose_nca_taz, amt ~ time | cohort + id,
                                 doseu = "mg")
intervals_taz <- data.frame(
  start    = start_ss,
  end      = end_ss,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  auclast  = TRUE,
  cav      = TRUE,
  ctrough  = TRUE
)
nca_taz <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj_taz, dose_obj_taz,
                                          intervals = intervals_taz))
knitr::kable(summary(nca_taz),
             caption = "Simulated steady-state tazobactam NCA by cohort.")
```

## Comparison against Cohen-Wolkowiez 2014 Table 5

Cohen-Wolkowiez 2014 Table 5 reports per-cohort empirical-Bayesian piperacillin
post-hoc summaries (CL, V, half-life, and predicted concentrations at 50%, 75%,
and the end of the dosing interval at steady state). The table below compares
the published per-cohort median Cmin,ss (predicted concentration at the end of
the dosing interval at steady state) against the simulated typical-value
Cmin,ss from the virtual cohort. Differences are expected because the published
values are empirical-Bayesian medians across heterogeneous patients within each
cohort whereas the simulated values represent a single typical-value subject at
the cohort median demographics.

```{r pip-cmin-comparison}
sim_typical_pip <- rxode2::rxSolve(rxode2::zeroRe(mod_pip),
                                   events = events_pip, keep = "cohort")
sim_typical_pip <- as.data.frame(sim_typical_pip)
cmin_sim <- sim_typical_pip |>
  dplyr::filter(time == end_ss, !is.na(Cc)) |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(Cmin_ss_sim_mgL = round(median(Cc), 1), .groups = "drop")

published <- tibble::tibble(
  cohort                  = c("cohort 1", "cohort 2", "cohort 3", "cohort 4"),
  Cmin_ss_pub_median_mgL  = c(38.3, 34.6, 28.0, 31.1),
  Cmin_ss_pub_range_mgL   = c("7.4-48.4", "5.5-48.4", "1.3-48.4", "1.5-60.2")
)
knitr::kable(
  dplyr::left_join(published, cmin_sim, by = "cohort"),
  caption = paste(
    "Published (Cohen-Wolkowiez 2014 Table 5) vs simulated typical-value",
    "piperacillin Cmin,ss by cohort."
  )
)
```

# Assumptions and deviations

- **PMA unit conversion.** The canonical covariate `PAGE` is in months; the
  source paper reports postmenstrual age in weeks. Both model files convert
  back to weeks inside `model()` via `pma_wk <- PAGE * 4.35` so the published
  exponent (1.76 for piperacillin, 1.35 for tazobactam) and the 33-week
  reference apply unchanged. The 4.35 weeks-per-month constant matches the
  canonical `PAGE` definition (`GA_weeks / 4.35 + postnatal_months`).
- **WT exponent fixed at 1.** The source paper explicitly held the linear
  body-weight exponent at 1 for both CL and V after rejecting free body-size
  estimation and 3/4-power allometric scaling for lack of model-fit
  improvement (Population PK model building paragraph). Both model files
  wrap the exponents in `fixed(1.0)` to preserve this provenance.
- **No between-subject variability on V.** The source paper fixed the between-
  subject variability on V to zero after observing a boundary-value estimate
  in the base model (Population PK model building paragraph). Neither model
  file declares an `etalvc` term; the only random effect is `etalcl`.
- **Tazobactam GENT effect is a confounded marker.** The Population PK model
  building paragraph notes that the concurrent gentamicin coadministration
  covariate was confounded by PMA / PNA differences between infants who
  received gentamicin (older, higher PMA / PNA) and those who did not
  (younger). The 1.52-fold multiplier on CL should therefore be interpreted as
  a marker of an older / sicker subset rather than a mechanistic drug-drug
  interaction. The source paper used the SCR + GENT + PMA tazobactam model on
  goodness-of-fit grounds but explicitly did **not** use it for dosing
  simulations or recommendations.
- **Virtual cohort approximates Table 1 medians.** The simulation cohort
  assigns per-cohort median PMA, WT, SCR, and gentamicin coadministration
  from Cohen-Wolkowiez 2014 Table 1; no covariate distribution is sampled
  (each cohort is represented by 50 copies of the same typical subject). This
  matches the figure scope (per-cohort VPC + typical-value Bayesian-CL curve)
  and the Cmin,ss comparison; a richer covariate-sampling cohort would be
  needed to reproduce the full per-cohort empirical-Bayesian Cmin,ss range.
- **Race / ethnicity not reported.** The source paper does not report race or
  ethnicity for the study cohort, so `population$race_ethnicity` is set to
  "Not reported in the source paper" for both model files and no race
  covariate is included.
- **Dried-blood-spot factor not extracted.** Cohen-Wolkowiez 2014 Tables 3
  and 4 also report a DBS matrix factor (0.38 for piperacillin, 0.48 for
  tazobactam) and a DBS additive residual SD from the joint plasma+DBS
  model. These are not included in either model file because the structural
  models extracted here are the **plasma-only** PMA-based piperacillin model
  (Table 3, "Point estimate" column) and the **plasma-only** SCR+GENT+PMA
  tazobactam model (Table 4, "Point estimate" column). The DBS-augmented
  variants would be a separate extraction (a future paired model file with
  an `assay_type` covariate to switch between plasma and DBS residual
  structures, or as a multi-output observation model).
