---
title: "Ceftazidime (Bulitta 2010)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ceftazidime (Bulitta 2010)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Bulitta_2010_ceftazidime"))$reference`
* Description: `r rxode2::rxode(readModelDb("Bulitta_2010_ceftazidime"))$description`
* Article: [Antimicrob Agents Chemother 2010;54(3):1275-1282](https://doi.org/10.1128/AAC.00936-09)

# Population

Bulitta 2010 was a single-dose, single-centre, open, parallel-group study
in 8 cystic fibrosis (CF) patients (4 male, 4 female; age 10-45 years;
total body weight 14.2-73.5 kg; fat-free mass 13.9-57.7 kg) and 7 healthy
volunteers (HV; 4 male, 3 female; age 19-33 years; total body weight
56-71 kg; fat-free mass 46.4-60.5 kg). All 15 subjects were Caucasian.
CF patients were notably smaller and leaner than the HV cohort (Bulitta
2010 Table 1). All subjects received a single 5-minute intravenous
infusion of 2 g ceftazidime, with three exceptions (one CF patient
received 1 g, one 1.5 g, and one 3 g per clinical judgment). 21 plasma
concentrations per subject were obtained between 0 and 12 h after the
end of infusion via HPLC with an internal standard (linear range
0.6-200 mg/L). The study was performed in 1983.

The population PK model is a three-compartment model with linear IV
elimination, allometric scaling on clearance (fixed exponent 0.75) and
linear scaling on volumes (exponent 1.0) by fat-free mass, with a
disease-group scale factor for the CF-vs-healthy contrast on total
clearance and on volume of distribution at steady state. Fat-free mass
was calculated per subject from the Janmahasatian et al. 2005 formula.
The same model was fit in three estimation engines (NONMEM, S-ADAPT,
NPAG); this packaged model reproduces the NONMEM column of Bulitta 2010
Table 3.

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

# Source trace

Every numeric value in `ini()` carries an in-file comment pointing to
the Bulitta 2010 source location. The table below collects them in one
place for review.

| Equation / parameter             | Value           | Source location                                                |
|----------------------------------|-----------------|----------------------------------------------------------------|
| `lcl`  (CL, CF at 53 kg FFM)     | 7.82 L/h        | Table 3, "CL" row, NONMEM CF                                   |
| `lvc`  (V1, CF at 53 kg FFM)     | 5.73 L          | Table 3, "V1" row, NONMEM CF                                   |
| `lvp`  (V2, CF at 53 kg FFM)     | 3.92 L          | Table 3, "V2" row, NONMEM CF                                   |
| `lvp2` (V3, CF at 53 kg FFM)     | 3.16 L          | Table 3, "V3" row, NONMEM CF                                   |
| `lq`   (CLic_shallow, both grp)  | 27.9 L/h        | Table 3, "CLic shallow" row, NONMEM                            |
| `lq2`  (CLic_deep, both grp)     | 2.57 L/h        | Table 3, "CLic deep" row, NONMEM                               |
| `e_ffm_cl`                       | 0.75 (fixed)    | Methods, "Body size model": clearance exponent fixed to 0.75   |
| `e_ffm_vc`                       | 1.00 (fixed)    | Methods, "Body size model": volume exponent fixed to 1.0       |
| `e_healthy_cl`                   | log(1/1.17)     | Table 4 row 5: FCYFCL = 1.17 for FFM allometric                |
| `e_healthy_vc / vp / vp2`        | log(1/1.01)     | Table 4 row 5: FCYFVSS = 1.01 for FFM allometric               |
| `etalcl`  (28% CV on CL)         | 0.07555         | Table 3, NONMEM CV(CL) -> log(0.28^2 + 1)                      |
| `etalvc`  (45% CV on V1)         | 0.18438         | Table 3, NONMEM CV(V1) -> log(0.45^2 + 1)                      |
| `etalvp`  (25% CV on V2)         | 0.06062         | Table 3, NONMEM CV(V2) -> log(0.25^2 + 1)                      |
| `etalvp2` (29% CV on V3)         | 0.08075         | Table 3, NONMEM CV(V3) -> log(0.29^2 + 1)                      |
| Correlation r(V1, V2) = -0.67    | -0.07084 cov    | Table 3 footnote e                                             |
| Correlation r(V1, V3) =  0.56    |  0.06834 cov    | Table 3 footnote e                                             |
| Correlation r(V2, V3) = -0.59    | -0.04128 cov    | Table 3 footnote e                                             |
| `propSd` (CV C, proportional)    | 0.122           | Table 3 footer, "CV C" row, NONMEM column                      |
| `addSd`  (SD C, additive)        | 0.059 mg/L      | Table 3 footer, "SD C (mg/liter)" row, NONMEM column           |
| FFM reference                    | 53 kg           | Table 3 footnote b: "subjects of standard size, FFM = 53 kg"   |
| Combined add + prop residual     | n/a             | Methods, "Between-subject variability and observation model"   |
| Three-compartment IV structural  | n/a             | Methods, "Structural model"; Results paragraph 2               |

IIV variance derivation. The Bulitta 2010 Methods section states an
exponential (log-normal) IIV model; for log-normal etas, `omega^2 =
log(CV^2 + 1)`. The covariance entries in the `etalvc + etalvp + etalvp2`
block are `cov = r * sqrt(var_i * var_j)`.

# Virtual cohort

Original observed data are not publicly available. The cohort below
mirrors the Bulitta 2010 demographics: an 8-patient CF arm and a
7-volunteer HV arm scaled up to 200 simulated subjects each by sampling
fat-free mass from a log-normal distribution matched to the per-group
medians and ranges in Table 1. All simulated subjects receive a single
2 g ceftazidime IV infusion over 5 minutes.

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

n_sub <- 200L
dose_mg <- 2000
infusion_h <- 5 / 60

# FFM distributions matched to Table 1: CF median 35.9 kg (range
# 13.9-57.7); HV median 54.0 kg (range 46.4-60.5).
sample_ffm <- function(n, median_kg, sd_log) {
  exp(rnorm(n, mean = log(median_kg), sd = sd_log))
}

build_arm <- function(label, dis_healthy, ffm_kg, id_offset) {
  ids <- id_offset + seq_len(length(ffm_kg))

  dose_rows <- tibble(
    id          = ids,
    time        = 0,
    evid        = 1L,
    amt         = dose_mg,
    cmt         = "central",
    rate        = dose_mg / infusion_h,
    cohort      = label,
    DIS_HEALTHY = dis_healthy,
    FFM         = ffm_kg
  )

  obs_times <- c(0, 5, 10, 15, 20, 30, 45, 60, 90) / 60
  obs_times <- c(obs_times, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12)
  # Add a denser early grid for Cmax extraction at end of infusion.
  obs_times <- sort(unique(c(obs_times, seq(0, 0.5, by = 0.025))))

  obs_rows <- tidyr::expand_grid(id = ids, time = obs_times) |>
    mutate(
      evid        = 0L,
      amt         = 0,
      cmt         = NA_character_,
      rate        = 0,
      cohort      = label,
      DIS_HEALTHY = dis_healthy,
      FFM         = ffm_kg[match(id, ids)]
    )

  bind_rows(dose_rows, obs_rows) |> arrange(id, time, desc(evid))
}

ffm_cf <- sample_ffm(n_sub, median_kg = 35.9, sd_log = 0.35)
ffm_hv <- sample_ffm(n_sub, median_kg = 54.0, sd_log = 0.05)

events <- bind_rows(
  build_arm("CF", dis_healthy = 0L, ffm_kg = ffm_cf, id_offset =      0L),
  build_arm("HV", dis_healthy = 1L, ffm_kg = ffm_hv, id_offset = n_sub)
)

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

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

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("cohort", "DIS_HEALTHY", "FFM")
) |> as.data.frame()
```

A typical-value simulation (random effects zeroed) is used for direct
comparison against the Bulitta 2010 Table 3 reference estimates at the
standard FFM of 53 kg.

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()

events_typ <- bind_rows(
  build_arm("CF_53kg", dis_healthy = 0L, ffm_kg = rep(53, 5L), id_offset =  0L),
  build_arm("HV_53kg", dis_healthy = 1L, ffm_kg = rep(53, 5L), id_offset = 10L)
)

sim_typ <- rxode2::rxSolve(
  mod_typical,
  events = events_typ,
  keep   = c("cohort", "DIS_HEALTHY", "FFM")
) |> as.data.frame()
```

# Replicate Bulitta 2010 Figure 1 (visual predictive check)

Bulitta 2010 Figure 1 shows the observed concentrations overlaid on the
10/25/50/75/90% prediction percentiles from the linear three-compartment
model based on FFM. The figure below reproduces the simulated percentile
envelope from the packaged model, split by disease group; observed data
are not available for overlay (no individual concentration data in the
public source).

```{r figure-1, fig.width = 7.5, fig.height = 4.5}
sim |>
  group_by(cohort, time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q25 = quantile(Cc, 0.25, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q75 = quantile(Cc, 0.75, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  filter(time > 0, time <= 12) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.25) +
  geom_line() +
  facet_wrap(~ cohort, labeller = labeller(
    cohort = c("CF" = "Cystic fibrosis patients",
               "HV" = "Healthy volunteers")
  )) +
  scale_y_log10() +
  labs(
    x = "Time after end of infusion (h)",
    y = "Ceftazidime Cc (mg/L, log scale)",
    title = "Simulated ceftazidime VPC by disease group",
    subtitle = "Single 2 g IV over 5 min; FFM sampled from Bulitta 2010 Table 1 distributions",
    caption = paste(
      "Light band: 10-90% percentile envelope;",
      "darker band: 25-75% percentile envelope.",
      "Compare against Bulitta 2010 Figure 1."
    )
  )
```

# Typical-value reference (Bulitta 2010 Table 3)

Reading the Table 3 reference estimates (NONMEM column) at FFM = 53 kg
back from the simulated typical-value curves confirms that the packaged
parameters round-trip to the published values.

```{r table3-check}
table3 <- tibble::tibble(
  param  = c("CL (L/h)", "V1 (L)", "V2 (L)", "V3 (L)",
             "Q (L/h)", "Q2 (L/h)"),
  CF_published = c(7.82, 5.73, 3.92, 3.16, 27.9, 2.57),
  HV_published = c(6.68, 5.67, 3.88, 3.13, 27.9, 2.57)
)

mod_params <- mod_typical$theta

ffm_ref <- 53
size_cl <- (ffm_ref / 53)^mod_params[["e_ffm_cl"]]
size_v  <- (ffm_ref / 53)^mod_params[["e_ffm_vc"]]

cf_calc <- c(
  CL = exp(mod_params[["lcl"]] )                            * size_cl,
  V1 = exp(mod_params[["lvc"]] )                            * size_v,
  V2 = exp(mod_params[["lvp"]] )                            * size_v,
  V3 = exp(mod_params[["lvp2"]])                            * size_v,
  Q  = exp(mod_params[["lq"]]  )                            * size_cl,
  Q2 = exp(mod_params[["lq2"]] )                            * size_cl
)
hv_calc <- c(
  CL = exp(mod_params[["lcl"]]  + mod_params[["e_healthy_cl"]] ) * size_cl,
  V1 = exp(mod_params[["lvc"]]  + mod_params[["e_healthy_vc"]] ) * size_v,
  V2 = exp(mod_params[["lvp"]]  + mod_params[["e_healthy_vp"]] ) * size_v,
  V3 = exp(mod_params[["lvp2"]] + mod_params[["e_healthy_vp2"]]) * size_v,
  Q  = exp(mod_params[["lq"]]   )                                * size_cl,
  Q2 = exp(mod_params[["lq2"]]  )                                * size_cl
)

table3_cmp <- table3 |>
  mutate(
    CF_packaged = round(cf_calc, 2),
    HV_packaged = round(hv_calc, 2)
  )
knitr::kable(
  table3_cmp,
  caption = "Bulitta 2010 Table 3 NONMEM reference estimates at FFM = 53 kg vs the packaged model's typical-value parameters."
)
```

# PKNCA validation

Run PKNCA over the 200-subject-per-arm stochastic cohort and compare the
per-group median NCA parameters against Bulitta 2010 Table 2
(noncompartmental analysis in the observed 15 subjects). Table 2 is not
size-normalised, so the simulated NCA must reflect the full FFM
distribution within each arm rather than a single typical FFM.

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

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

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

intervals <- data.frame(
  start       = 0,
  end         = Inf,
  cmax        = TRUE,
  tmax        = TRUE,
  aucinf.obs  = TRUE,
  half.life   = TRUE,
  cl.obs      = TRUE,
  mrt.iv.obs  = TRUE,
  vss.iv.obs  = TRUE
)

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

nca_tbl <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "aucinf.obs", "half.life",
                          "cl.obs", "vss.iv.obs", "mrt.iv.obs")) |>
  mutate(cohort = as.character(cohort)) |>
  group_by(cohort, PPTESTCD) |>
  summarise(
    median  = median(PPORRES, na.rm = TRUE),
    p_min   = min(PPORRES, na.rm = TRUE),
    p_max   = max(PPORRES, na.rm = TRUE),
    .groups = "drop"
  )

nca_wide <- nca_tbl |>
  tidyr::pivot_wider(
    names_from  = cohort,
    values_from = c(median, p_min, p_max),
    names_glue  = "{cohort}_{.value}"
  )

knitr::kable(
  nca_wide,
  digits  = 2,
  caption = "Simulated NCA per-group medians and min-max ranges across 200 subjects per arm. Compare against Bulitta 2010 Table 2 (CF median CL 5.37 L/h, Vss 9.14 L, peak 445 mg/L, t1/2 1.48 h, MRT 1.54 h; HV median CL 6.59 L/h, Vss 14.3 L, peak 210 mg/L, t1/2 1.94 h, MRT 2.10 h)."
)
```

## Comparison against Bulitta 2010 Table 2

```{r table2-comparison}
table2 <- tibble::tibble(
  PPTESTCD     = c("cl.obs", "vss.iv.obs", "cmax", "half.life", "mrt.iv.obs"),
  parameter    = c("Total clearance (L/h)",
                   "Vss (L)",
                   "Peak concentration (mg/L)",
                   "Terminal half-life (h)",
                   "Mean residence time (h)"),
  CF_published = c(5.37, 9.14, 445, 1.48, 1.54),
  HV_published = c(6.59, 14.3, 210, 1.94, 2.10)
)

sim_medians <- nca_tbl |>
  filter(PPTESTCD %in% table2$PPTESTCD) |>
  select(cohort, PPTESTCD, median) |>
  tidyr::pivot_wider(names_from = cohort, values_from = median,
                     names_glue = "{cohort}_simulated")

table2_cmp <- table2 |>
  left_join(sim_medians, by = "PPTESTCD") |>
  select(parameter,
         CF_published, CF_simulated,
         HV_published, HV_simulated) |>
  mutate(
    CF_pct_diff = round(100 * (CF_simulated - CF_published) / CF_published, 1),
    HV_pct_diff = round(100 * (HV_simulated - HV_published) / HV_published, 1)
  )
knitr::kable(
  table2_cmp,
  digits  = 2,
  caption = "Per-group median NCA from the packaged model vs Bulitta 2010 Table 2. CF n = 8, HV n = 7 in the source; here both arms are 200-subject virtual cohorts matching the per-arm FFM distributions."
)
```

The simulated medians track the published Table 2 medians: the
CF-vs-HV ordering (CF has higher peak, lower clearance per kg WT, shorter
MRT) is reproduced. Absolute differences > 20% relative to Table 2 are
expected for individual NCA values because the original Table 2 was
computed from 8 CF and 7 HV subjects with very wide ranges (e.g., peak
151-1200 mg/L in CF), and the virtual cohort uses a smooth log-normal
FFM distribution rather than the discrete actual demographics; per-group
medians from the virtual cohort therefore approximate the population
average without recovering each Table 2 cell exactly.

# Body-size scaling

Bulitta 2010 Table 4 / Table 5 emphasises that allometric scaling by FFM
reduced the unexplained between-subject variance by 32% on CL and by
18-26% on the peripheral volumes relative to linear scaling by total
weight. The figure below shows how typical-value CL changes with FFM in
each disease group under the packaged allometric-FFM model.

```{r ffm-scaling, fig.width = 7, fig.height = 4}
ffm_grid <- seq(15, 75, by = 1)

scaling_df <- tibble::tibble(
  FFM = rep(ffm_grid, 2),
  cohort = rep(c("CF", "HV"), each = length(ffm_grid))
) |>
  mutate(
    DIS_HEALTHY = as.integer(cohort == "HV"),
    CL = exp(mod_params[["lcl"]] + DIS_HEALTHY * mod_params[["e_healthy_cl"]]) *
         (FFM / 53)^mod_params[["e_ffm_cl"]]
  )

ggplot(scaling_df, aes(FFM, CL, colour = cohort)) +
  geom_line(linewidth = 1) +
  labs(
    x = "Fat-free mass (kg)",
    y = "Typical-value total clearance (L/h)",
    title = "Allometric FFM scaling of ceftazidime CL",
    subtitle = "Typical-value CL = CL_ref * (FFM/53)^0.75; FCYFCL = 1.17 separates CF from HV",
    caption = "At FFM = 53 kg: CL = 7.82 L/h (CF) and 6.68 L/h (HV); Bulitta 2010 Table 3 NONMEM."
  )
```

# Assumptions and deviations

- **Disease-group orientation re-encoded to `DIS_HEALTHY`.** Bulitta
  2010 parameterises the disease contrast through scale factors FCYFCL =
  1.17 and FCYFVSS = 1.01 with the healthy volunteer cohort serving as
  the structural reference (`CL_CF = CL_HV * 1.17`). The packaged model
  uses the canonical `DIS_HEALTHY` covariate (1 = healthy, 0 = patient)
  with CF as the structural reference (typical-value parameters equal
  the Table 3 CF column) and `e_healthy_<param> = log(1 / FCYF)`
  shifting to the HV estimates when `DIS_HEALTHY = 1`. This preserves
  the published parameter values exactly while matching the existing
  `DIS_HEALTHY` precedent (Goel 2016 sonidegib, Nikanjam 2019
  siltuximab, Yoneyama 2017 emicizumab).
- **Shared FCYFVSS encoded as three separate parameters.** Bulitta 2010
  estimates a single FCYFVSS = 1.01 applied to V1, V2, and V3. The
  packaged model assigns the same value to three canonical covariate
  effect parameters (`e_healthy_vc`, `e_healthy_vp`, `e_healthy_vp2`)
  to satisfy the `e_<cov>_<param>` naming convention. If users re-fit
  this model, the three coefficients should be tied to preserve the
  paper's parameterisation.
- **Intercompartmental clearances are not stratified by disease group.**
  Bulitta 2010 Table 3 reports identical Q and Q2 values for both
  cohorts (27.9 L/h and 2.57 L/h, NONMEM column); no FCYFCL is applied
  to Q or Q2. The packaged model preserves this structure.
- **FFM is the only structural body-size covariate.** Bulitta 2010
  Methods evaluated five body-size models (no scaling, linear WT,
  allometric WT, linear FFM, allometric FFM). The final model packaged
  here is the allometric FFM model (Table 4 row 5; Table 5 row 4); the
  other variants are not separately packaged. Users who want to
  benchmark a linear-WT scaling against the same compound need to
  rewrite the size-scaling lines in `model()`.
- **Allometric exponents fixed.** The clearance exponent (0.75) and
  volume exponent (1.0) were both fixed during estimation (Bulitta 2010
  Methods, "Body size model"); the packaged model wraps both `e_ffm_cl`
  and `e_ffm_vc` in `fixed()` to preserve this status.
- **Renal function was not retained as a covariate.** Bulitta 2010
  Discussion: "A potential limitation of our model is that it did not
  account for renal function." All 15 subjects had normal renal
  function. The packaged model has no `CRCL` covariate; users with
  renally impaired patients should consult a renal-impairment-specific
  ceftazidime model (e.g., dialysis cohorts) before applying.
- **Virtual cohort FFM distribution.** The vignette samples FFM from a
  log-normal centred on the Table 1 medians (CF 35.9 kg, HV 54.0 kg).
  The CF arm's log-SD of 0.35 broadens the right-tail to approximate
  the observed range 13.9-57.7 kg; the HV log-SD of 0.05 holds the
  cohort tight around 54 kg in line with the published 46.4-60.5 kg
  range. These distributional choices are convenience values for the
  validation simulation and do not affect the packaged model parameters.
- **Dose harmonised to 2 g.** Three of the eight CF patients received
  doses of 1, 1.5, or 3 g per clinical judgment (Methods). The virtual
  cohort uses 2 g for all subjects to simplify the NCA comparison
  against Table 2 (Table 2 footnote a normalises CF peak concentrations
  to 2 g for the same reason).
- **No published errata identified.** A search of the journal landing
  page (https://journals.asm.org/doi/10.1128/AAC.00936-09) for
  corrections / corrigenda did not return an erratum for Bulitta 2010
  doi:10.1128/AAC.00936-09. The packaged values are the original Table
  3 NONMEM estimates.
