---
title: "Cefuroxime axetil (Bulitta 2009)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cefuroxime axetil (Bulitta 2009)}
  %\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 rxode2::rxode(readModelDb("Bulitta_2009_cefuroxime_axetil"))$reference`
* Description: `r rxode2::rxode(readModelDb("Bulitta_2009_cefuroxime_axetil"))$description`
* Article: <https://doi.org/10.1128/AAC.00054-09>

# Population

Bulitta 2009 enrolled 24 healthy male Caucasian volunteers (mean age 24.5
years, range 18-31; mean weight 73.8 kg, range 58.2-93.6; mean height 179
cm, range 166-193) in a single-dose, single-centre cefuroxime axetil
bioavailability study at IBMP (Nurnberg, Germany). Each volunteer received
a 300.72 mg cefuroxime axetil oral suspension (equivalent to 250 mg
cefuroxime) with 240 mL water immediately after a standardized high-fat
breakfast (four slices of crisp bread, margarine, cheese, jam, fruit tea,
and 100 mL of 3.5% milk; described in 'Study design and drug
administration'). Plasma samples for cefuroxime were drawn at predose and
30, 60, 90 min and 2, 2.33, 2.67, 3, 3.33, 3.67, 4, 4.5, 5, 6, 8, 10, and
12 h post-dose, then assayed by validated LC-MS/MS (LLOQ 0.00900 mg/L;
'Drug analysis'). The semiphysiological PK model was fitted in three
programs (NONMEM VI FOCE-I, S-ADAPT 1.55 with importance-sampling
parametric Monte Carlo EM, and NPAG via USC*PACK 12.00); the packaged
model reproduces the S-ADAPT parametric fit because that program is the
one the authors recommend as the best parametric solution (Discussion,
'parametric VPC using estimates from S-ADAPT had the best predictive
performance among the parametric VPCs').

The same information is available programmatically via the model's
`population` metadata
(`readModelDb("Bulitta_2009_cefuroxime_axetil")$population`).

# Source trace

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

| Equation / parameter                                              | Value                  | Source location                                               |
|-------------------------------------------------------------------|------------------------|---------------------------------------------------------------|
| `lcl`     (CL/F, apparent clearance)                              | log(21.7 L/h)          | Table 2 'S-ADAPT Population mean'                             |
| `lvc`     (V/F, apparent central volume)                          | log(38.7 L)            | Table 2 'S-ADAPT Population mean'                             |
| `lvmax`   (Vmax0/dose, max gastric-release rate as dose fraction) | log(0.505 1/h)         | Table 2 'S-ADAPT Population mean'                             |
| `lkm`     (Km/dose, Michaelis constant as dose fraction)          | log(0.426)             | Table 2 'S-ADAPT Population mean' (42.6% of dose)             |
| `ltabs`   (T_abs, absorption half-life from intestine)            | log(9.34 min)          | Table 2 'S-ADAPT Population mean' (footnote f for ka)         |
| `ltc50`   (TC50, time past meal at 50% Vmax change)               | log(1.61 h)            | Table 2 'S-ADAPT Population mean'                             |
| `lgtemax` (Lg_Emax, logit-transformed max fractional Vmax change) | -0.762                 | Table 2 'S-ADAPT Population mean' (footnote on logit transform) |
| `lhill`   (gamma, Hill coefficient for time-dependent Vmax)       | fixed(log(10))         | Table 2 footnote g                                            |
| Full 7x7 BSV variance-covariance matrix                           | see Table 3            | Table 3 'Variance-covariance matrix ... S-ADAPT'              |
| `propSd`  (proportional residual SD)                              | 0.0851                 | Table 2 'S-ADAPT' final-row                                   |
| `addSd`   (additive residual SD, mg/L)                            | 0.0029                 | Table 2 'S-ADAPT' final-row                                   |
| ODE dA1/dt: stomach -> intestine via Michaelis-Menten release     | n/a                    | 'Structural model', differential equations                    |
| ODE dA2/dt: intestine -> central via first-order absorption       | n/a                    | 'Structural model', differential equations                    |
| ODE dA3/dt: central elimination CL/V                              | n/a                    | 'Structural model', differential equations                    |
| Time-dependent Vmax(TPM) = Vmax0 * (1 + Emax * TPM^g / (TC50^g + TPM^g)) | n/a              | 'Structural model', Vmax equation                             |
| Logit transform Emax = 10 * exp(Lg_Emax) / (1 + exp(Lg_Emax)) - 1  | n/a                   | 'Parameter variability and observation model'                 |
| Dose-fraction scaling Vmax = (Vmax0/dose) * dose and similarly Km  | n/a                   | Results, 'Vmax0 and Km were estimated and are reported as fractions of the cefuroxime dose' |

# Virtual cohort

Bulitta 2009 administered a single 250 mg cefuroxime equivalent dose with a
standardized meal. The virtual cohort here mirrors that design: 24 subjects
(matching the published N) each receive a single oral dose into the
`stomach` compartment at t = 0. No covariates enter the model, so the only
inter-subject differences come from the 7-dimensional log-normal IIV block
in Table 3.

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

n_subj <- 24L
t_end  <- 12          # h, matches the longest published sample (12 h)
dt_obs <- 0.1         # 6-minute grid

doses <- tibble(
  id   = seq_len(n_subj),
  time = 0,
  amt  = 250,        # mg cefuroxime, dosed into the stomach compartment
  cmt  = "stomach",
  evid = 1L,
  ii   = 0,
  ss   = 0,
  addl = 0
)

obs <- tibble(
  id   = rep(seq_len(n_subj), each = length(seq(0, t_end, by = dt_obs))),
  time = rep(seq(0, t_end, by = dt_obs), times = n_subj),
  amt  = 0,
  cmt  = NA_character_,
  evid = 0L,
  ii   = 0,
  ss   = 0,
  addl = 0
)

events <- bind_rows(doses, obs) |>
  arrange(id, time, desc(evid))

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

# Simulation

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

# Stochastic simulation with the full 7x7 IIV block (reproduces the
# published VPCs). rxode2 draws etas from the correlated log-normal
# specified by the model's ini() block.
sim <- rxode2::rxSolve(mod, events = events) |>
  as.data.frame()

# Typical-value (zero IIV) -- useful for replicating the population-mean
# concentration profile and the Vmax(TPM) modulator in Figure 1.
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical, events = events) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2 -- individual curve fits

Bulitta 2009 Figure 2 shows 24 individual curve-fit panels (NONMEM, S-ADAPT
and NPAG overlays) for the cefuroxime concentration-time profile after a
single oral dose. The panel below superimposes the 24 stochastic
realisations from the packaged S-ADAPT model and the typical-value (zero-
IIV) line, on the same 0-12 h window. The diverse 'sharp-peak',
'plateau-like-peak' and 'dual-peak' shapes that the paper identifies as a
key motivation for the semiphysiological absorption structure (5 of 24
biphasic, 8 of 24 plateau-like; see Results, second paragraph) appear
naturally from the M-M release + time-dependent Vmax in the simulated
profiles.

```{r figure-2-individuals}
sim |>
  ggplot(aes(time, Cc, group = id)) +
  geom_line(alpha = 0.5, colour = "steelblue") +
  geom_line(
    data = sim_typical, mapping = aes(time, Cc), inherit.aes = FALSE,
    colour = "black", size = 0.9
  ) +
  labs(
    x = "Time (h)",
    y = "Cefuroxime plasma concentration (mg/L)",
    caption = paste(
      "Light blue: 24 simulated subjects from the packaged S-ADAPT model.",
      "Black: typical-value (zero-IIV) profile.",
      "Compare with Bulitta 2009 Figure 2 -- diverse 'sharp', 'plateau-like'",
      "and 'dual peak' shapes."
    )
  ) +
  theme_bw()
```

## Figure 3 -- visual predictive check

Bulitta 2009 Figure 3 reports VPCs (10th, 50th, 90th percentiles) for the
final population PK model in each of the three programs. The reproduction
below uses 1000 virtual subjects from the packaged S-ADAPT model so the
percentiles are well-estimated, on the linear and log10 y-axis as in the
paper's figure and inset.

```{r figure-3-vpc}
set.seed(20260604)
n_vpc <- 1000L

events_vpc <- bind_rows(
  tibble(
    id   = seq_len(n_vpc),
    time = 0, amt = 250, cmt = "stomach", evid = 1L,
    ii = 0, ss = 0, addl = 0
  ),
  tibble(
    id   = rep(seq_len(n_vpc), each = length(seq(0, t_end, by = dt_obs))),
    time = rep(seq(0, t_end, by = dt_obs), times = n_vpc),
    amt  = 0, cmt = NA_character_, evid = 0L,
    ii = 0, ss = 0, addl = 0
  )
) |>
  arrange(id, time, desc(evid))

sim_vpc <- rxode2::rxSolve(mod, events = events_vpc) |>
  as.data.frame()

vpc_summary <- sim_vpc |>
  filter(!is.na(Cc), time > 0) |>
  group_by(time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.25, fill = "steelblue") +
  geom_line(size = 0.9) +
  labs(
    x = "Time (h)",
    y = "Cefuroxime plasma concentration (mg/L)",
    caption = paste(
      "Simulated 10th-90th percentile prediction interval (shaded) and",
      "median (solid) from 1000 virtual subjects. Replicates the layout of",
      "Bulitta 2009 Figure 3 (S-ADAPT panel)."
    )
  ) +
  theme_bw()
```

```{r figure-3-vpc-log}
ggplot(vpc_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.25, fill = "steelblue") +
  geom_line(size = 0.9) +
  scale_y_log10() +
  labs(
    x = "Time (h)",
    y = "Cefuroxime plasma concentration (mg/L, log scale)",
    caption = "Same VPC on a semilogarithmic axis (replicates the Figure 3 inset)."
  ) +
  theme_bw()
```

## Figure 1 inset -- time-dependent Vmax modulator

Figure 1 of Bulitta 2009 shows the V_max/V_max0 ratio over time-past-meal
(TPM) for representative Emax values (3, 1, -0.5) with TC50 = 2.5 h and
gamma = 10. The reproduction below plots V_max(TPM)/V_max0 = 1 + Emax *
TPM^g / (TC50^g + TPM^g) for those three Emax scenarios using the
paper's figure-1-inset parameters (TC50 = 2.5 h, gamma = 10), and for the
final S-ADAPT typical Emax to show what the actual modulator looks like.

```{r figure-1-modulator}
typical_emax <- with(list(lg = -0.762), 10 * exp(lg) / (1 + exp(lg)) - 1)
tc50_fig1    <- 2.5
gamma_fig1   <- 10
tpm          <- seq(0, 8, by = 0.05)

modulator <- function(tpm, emax, tc50, gamma) {
  1 + emax * tpm^gamma / (tc50^gamma + tpm^gamma)
}

mod_df <- bind_rows(
  tibble(label = "Emax = 3 (paper Fig. 1)",      tpm = tpm,
         value = modulator(tpm,  3.0, tc50_fig1, gamma_fig1)),
  tibble(label = "Emax = 1 (paper Fig. 1)",      tpm = tpm,
         value = modulator(tpm,  1.0, tc50_fig1, gamma_fig1)),
  tibble(label = "Emax = -0.5 (paper Fig. 1)",   tpm = tpm,
         value = modulator(tpm, -0.5, tc50_fig1, gamma_fig1)),
  tibble(label = sprintf("Typical S-ADAPT (Emax = %.2f, TC50 = 1.61 h)", typical_emax),
         tpm = tpm,
         value = modulator(tpm, typical_emax, 1.61, gamma_fig1))
)

ggplot(mod_df, aes(tpm, value, colour = label, linetype = label)) +
  geom_line(size = 0.9) +
  labs(
    x = "Time past meal (h)",
    y = "Vmax(TPM) / Vmax0",
    colour = NULL, linetype = NULL,
    caption = paste(
      "Replicates the Bulitta 2009 Figure 1 inset modulator curves",
      "(top three lines) and adds the S-ADAPT population-mean curve",
      "(bottom line) so the deployed modulator can be inspected next",
      "to the paper's illustrative values."
    )
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

# PKNCA validation

Bulitta 2009 Table 1 reports noncompartmental PK parameters for the same
250 mg single oral cefuroxime dose. The block below runs PKNCA on the
24-subject stochastic simulation (matching the published study size) and
compares the simulated mean +/- SD to the published values. A treatment
grouping is included for consistency with the skill convention even though
only one regimen is studied in the paper.

The integration window is extended to 48 h (about 35 cefuroxime half-lives
beyond the published 1.34 h terminal half-life) so the AUC0-inf estimator
in PKNCA has a clean terminal phase to fit. The full 7-dimensional log-
normal IIV in this model is wide enough that a handful of stochastic
draws still produce sluggish absorption profiles whose lambda.z fits
are unreliable; the comparison table reports the published values
alongside the mean +/- SD across the subjects whose terminal phase
PKNCA could resolve.

```{r pknca-events}
t_nca       <- 48
dt_nca      <- 0.25

events_nca <- bind_rows(
  tibble(
    id   = seq_len(n_subj),
    time = 0, amt = 250, cmt = "stomach", evid = 1L,
    ii = 0, ss = 0, addl = 0,
    treatment = "250 mg single oral dose"
  ),
  tibble(
    id   = rep(seq_len(n_subj), each = length(seq(0, t_nca, by = dt_nca))),
    time = rep(seq(0, t_nca, by = dt_nca), times = n_subj),
    amt  = 0, cmt = NA_character_, evid = 0L,
    ii = 0, ss = 0, addl = 0,
    treatment = "250 mg single oral dose"
  )
) |>
  arrange(id, time, desc(evid))

sim_nca_raw <- rxode2::rxSolve(
  mod, events = events_nca,
  keep = c("treatment")
) |>
  as.data.frame()
```

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

dose_df <- events_nca |>
  filter(evid == 1L, cmt == "stomach") |>
  select(id, time, amt, treatment)

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

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

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

# Mark subjects whose terminal-phase fit has a span ratio < 2 or
# r2 < 0.9 (standard PKNCA-quality thresholds). The Bulitta 2009 BSV
# on absorption is wide enough that a few stochastic draws produce
# subjects whose simulated tail is still absorption-dominated;
# their lambda.z, AUCinf, half-life, and CL/F are unreliable and
# are dropped from the published-comparison table downstream.
quality_tbl <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("span.ratio", "adj.r.squared")) |>
  pivot_wider(names_from = PPTESTCD, values_from = PPORRES, id_cols = id)
unreliable_ids <- quality_tbl |>
  filter(span.ratio < 2 | adj.r.squared < 0.9 | is.na(span.ratio)) |>
  pull(id)
knitr::kable(
  summary(nca_res),
  caption = "Simulated NCA parameters (PKNCA) from 24 virtual subjects."
)
```

## Comparison against the paper's NCA table

```{r pknca-comparison}
nca_tbl <- as.data.frame(nca_res$result)

# Cmax / Tmax are independent of the terminal-phase fit; AUCinf,
# half-life, and CL/F rely on lambda.z and so are gated on the
# unreliable_ids list computed above.
agg <- function(code, drop_unreliable = FALSE) {
  rows <- nca_tbl[nca_tbl$PPTESTCD == code, c("id", "PPORRES")]
  if (drop_unreliable) {
    rows <- rows[!(rows$id %in% unreliable_ids), , drop = FALSE]
  }
  vals <- rows$PPORRES
  vals <- vals[is.finite(vals) & !is.na(vals)]
  if (length(vals) == 0) return(c(mean = NA_real_, sd = NA_real_, n = 0))
  c(mean = mean(vals), sd = sd(vals), n = length(vals))
}

cmax_sim <- agg("cmax")
tmax_sim <- agg("tmax")
auc_sim  <- agg("aucinf.obs", drop_unreliable = TRUE)
hl_sim   <- agg("half.life",  drop_unreliable = TRUE)
clf_sim  <- agg("cl.obs",     drop_unreliable = TRUE)

comparison <- tibble(
  Parameter = c(
    "Cmax (mg/L)",
    "Tmax (h)",
    "AUC0-inf (mg.h/L)",
    "Terminal half-life (h)",
    "Apparent total clearance (L/h)"
  ),
  `Bulitta 2009 Table 1 (mean +/- SD)` = c(
    "2.64 +/- 0.64",
    "2.98 +/- 0.73",
    "11.9 +/- 2.49",
    "1.34 +/- 0.13",
    "21.8 +/- 4.29"
  ),
  `Simulated (mean +/- SD)` = sprintf(
    "%.2f +/- %.2f (n = %d)",
    c(cmax_sim["mean"], tmax_sim["mean"], auc_sim["mean"], hl_sim["mean"], clf_sim["mean"]),
    c(cmax_sim["sd"],   tmax_sim["sd"],   auc_sim["sd"],   hl_sim["sd"],   clf_sim["sd"]),
    c(cmax_sim["n"],    tmax_sim["n"],    auc_sim["n"],    hl_sim["n"],    clf_sim["n"])
  )
)

knitr::kable(
  comparison,
  caption = paste(
    "PKNCA estimates from a single stochastic 24-subject draw alongside",
    "the published noncompartmental analysis. Cmax and Tmax agree with",
    "the published means within about 10-30%; the wide BSV on the",
    "absorption parameters in Table 3 produces a wider Tmax",
    "distribution than the 24 actually-enrolled subjects showed in",
    "Table 1. For AUC0-inf, half-life, and CL/F the per-subject",
    "lambda.z fit can be unreliable when a stochastic draw produces a",
    "slow-absorption subject whose 'terminal' phase is still",
    "absorption-limited; the (n = ...) column reports how many of the",
    "24 subjects PKNCA could resolve."
  )
)
```

# Assumptions and deviations

- **Source of typical-value estimates.** The packaged model reproduces the
  S-ADAPT importance-sampling parametric Monte Carlo EM fit (Table 2,
  'S-ADAPT Population mean' column), and the full 7x7 BSV variance-
  covariance matrix from Table 3. Bulitta 2009 also reports NONMEM VI
  FOCE-I and NPAG fits; the authors note that the S-ADAPT VPC had the
  best predictive performance among the parametric fits (Results,
  'parametric VPCs') and that the NPAG nonparametric MCS had the overall
  best predictive performance, but the nonparametric support-point
  distribution is not a parametric (mean, variance-covariance)
  specification expressible in nlmixr2 -- so the packaged model uses the
  closest parametric counterpart, which is S-ADAPT.

- **Dose-fraction scaling.** Vmax0 and Km were estimated by the authors
  as fractions of the cefuroxime dose because the rate of stomach
  release is primarily a meal effect, not a dose effect (Discussion,
  'Vmax and Km are expressed as fractions of dose'). The packaged model
  recovers the absolute Vmax0 (mg/h) and Km (mg) inside `model()` by
  multiplying the fraction parameters by `podo(stomach)` (the
  most-recent dose amount into the stomach compartment). This preserves
  the dose-independent Tmax that the paper's analysis depends on.

- **Logit transform on Emax.** The paper transforms Emax through
  Lg_Emax with Emax = 10 * exp(Lg_Emax) / (1 + exp(Lg_Emax)) - 1 so that
  the individual Emax values stay in (-1, 9) -- Emax = -1 represents
  complete inhibition of gastric release, Emax = 0 means no change, and
  Emax > 0 means stimulation. The packaged model carries `lgtemax` (and
  its IIV `etalgtemax`) on the transformed scale and applies the back-
  transform inside `model()` so that the BSV variance from Table 3 maps
  directly without re-parameterisation.

- **Hill coefficient fixed at 10.** Table 2 footnote g reports that the
  Hill coefficient was initially estimated in the 10-15 range; for
  estimation stability the authors fixed it at 10 in the final fit
  with no notable effect on objective function or curve quality. The
  packaged model encodes this as `lhill <- fixed(log(10))`.

- **No covariates in the published model.** The 24-subject cohort was
  selected to be uniform (healthy male Caucasian adults age 18-31,
  weight 58.2-93.6 kg). No covariate effects on CL, V, or absorption
  parameters were tested or reported. The packaged model therefore
  has an empty `covariateData` list.

- **Single dose only.** The analysis used a single 250 mg cefuroxime
  oral dose with a standardized high-fat meal. Multi-dose simulations
  (e.g., for the q8h / q12h MCS that the paper performs for PTA
  prediction) rely on `tad(stomach)` to reset the time-past-meal
  modulator at each dose, which implicitly assumes each dose coincides
  with a comparable meal -- the paper makes the same assumption.

- **Time-of-meal vs time-of-dose.** Because the paper's protocol gives
  the dose immediately after the breakfast, TPM (time past meal) equals
  time since dose. The packaged model uses `tad(stomach)` as a direct
  proxy for TPM. For users who want to model a different meal-vs-dose
  offset (e.g. the fasted-state extrapolation discussed in the
  Discussion), this assumption no longer holds and the user must
  introduce an external TPM variable.
