---
title: "Tobramycin conventional intratracheal (Marier 2002) -- rat"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tobramycin conventional intratracheal (Marier 2002) -- rat}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Marier JF, Lavigne J, Ducharme MP. Pharmacokinetics and
  efficacies of liposomal and conventional formulations of tobramycin
  after intratracheal administration in rats with pulmonary
  *Burkholderia cepacia* infection. *Antimicrobial Agents and
  Chemotherapy*. 2002;46(12):3776-3781.
* Article: [doi:10.1128/aac.46.12.3776-3781.2002](https://doi.org/10.1128/aac.46.12.3776-3781.2002)

This vignette validates the **conventional** (Tobi inhalation
solution) formulation arm. The **liposomal** formulation arm is in the
sibling vignette
[`Marier_2002_tobramycin_rat_liposomal`](Marier_2002_tobramycin_rat_liposomal.html);
they share the same ODE structure and differ only in their fitted
parameter values. Marier 2002's primary scientific finding is the
contrast between the two formulations: the liposomal formulation
shows slower absorption (k_a 2.30 vs 4.68 /h), slower elimination
(k_20 0.155 vs 1.10 /h), and ~8-fold greater AUC (6,811 vs 821
ug.h / lungs) than the conventional formulation.

# Population

The model was developed from 36 evaluable male Sprague-Dawley rats
(weight 175-225 g; Charles River, Saint-Constant, Quebec; 39 dosed --
2 died during the surgical procedure 3 h post-dose and the infection
level could not be confirmed in 1 animal). All rats had chronic
pulmonary *B. cepacia* infection (strain BC 1368, genomovar III,
tobramycin MIC 128 ug/mL) established 6 days pre-dose by
intratracheal inoculation of agar beads containing 10^6 CFU.

A single 1,200 ug intratracheal dose of conventional Tobi tobramycin
inhalation solution (PathoGenesis Canada) was administered as 100 uL
via calibrated pipette followed by an air bolus. Sampling was as in
the liposomal arm (13 time points x 3 rats; cardiac-puncture
exsanguination under ketamine + xylazine; HPLC-UV quantitation of
homogenised lung tissue).

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

# Source trace

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| ka (`lka`)           | 4.68 /h (CV 38.7%) | Table 1 (Conventional column) |
| FL (`lfdepot`)       | 0.77 (CV 17.2%) | Table 1 |
| k_23 (`lk12`)        | 5.58 /h (CV 43.9%) | Table 1 |
| k_32 (`lk21`)        | 0.330 /h (CV 37.8%) | Table 1 |
| k_20 (`lkel`)        | 1.10 /h (CV 19.8%) | Table 1 |
| IIV (diagonal CV%)   | as in `ini()` comments | Table 1 (CV column) |
| Residual variability | 18.4% (combined) | Results paragraph after Table 1 |
| Dose                 | 1,200 ug intratracheal single bolus | Methods (Experimental infection and antibiotic treatment) |
| Sampling             | 0.5, 1, 2, 3, 5, 7, 9, 11, 13, 15, 16, 18, 24 h | Methods (Sample collection and bacterial growth) |
| Derived t_1/2 alpha  | 0.12 h (CV 17.1%) | Table 1 (analytic, from individual Bayesian estimates) |
| Derived t_1/2 beta   | 12.9 h (CV 31.4%) | Table 1 |
| Derived AUC          | 821 ug.h / lungs (CV 15.0%) | Table 1 |

# Loading the model

```{r load}
mod <- readModelDb("Marier_2002_tobramycin_rat_conventional")
mod_typical <- rxode2::zeroRe(mod)
```

# Replicating Figure 1 (conventional panel)

Figure 1 of Marier 2002 plots the observed and individually-predicted
amounts of tobramycin in the lungs of rats over 24 h after a single
1,200 ug intratracheal dose. The conventional panel shows a fast
absorption peak followed by a fast terminal phase. The chunk below
builds an event table for a single virtual rat, simulates the
typical-value trajectory, and shows the profile on log scale.

```{r figure-1-typical}
sim_grid <- data.frame(
  id   = 1L,
  time = c(0, seq(0.05, 24, by = 0.05)),
  evid = 0L,
  amt  = 0,
  cmt  = "Cc"
)
dose_row <- data.frame(
  id   = 1L,
  time = 0,
  evid = 1L,
  amt  = 1200,
  cmt  = "depot"
)
events_typ <- dplyr::bind_rows(dose_row, sim_grid)
sim_typ <- as.data.frame(rxode2::rxSolve(mod_typical, events_typ))

ggplot(sim_typ, aes(time, Cc)) +
  geom_line(linewidth = 0.8, colour = "#e7298a") +
  scale_y_log10(limits = c(0.1, 2000)) +
  labs(x = "Time after intratracheal dose (h)",
       y = "Amount of tobramycin in lungs (ug, log scale)",
       title = "Conventional tobramycin -- typical-value profile",
       caption = paste("Replicates the conventional panel of Figure 1 of",
                       "Marier 2002. The terminal slope corresponds to the",
                       "published t_1/2 beta = 12.9 h.")) +
  theme_minimal(base_size = 11)
```

# Side-by-side comparison with the liposomal formulation

The paper's primary finding is the marked pharmacokinetic contrast
between the two formulations. The chunk below simulates the
typical-value profile of each and overlays them on the same axes.

```{r side-by-side, fig.height = 4.5}
mod_lip <- readModelDb("Marier_2002_tobramycin_rat_liposomal") |>
  rxode2::zeroRe()
mod_con <- mod_typical

sim_lip <- as.data.frame(rxode2::rxSolve(mod_lip, events_typ)) |>
  dplyr::mutate(formulation = "Liposomal")
sim_con <- as.data.frame(rxode2::rxSolve(mod_con, events_typ)) |>
  dplyr::mutate(formulation = "Conventional")
sim_both <- dplyr::bind_rows(sim_lip, sim_con)

ggplot(sim_both, aes(time, Cc, colour = formulation)) +
  geom_line(linewidth = 0.8) +
  scale_y_log10(limits = c(0.1, 2000)) +
  scale_colour_manual("Formulation",
                      values = c("Liposomal"   = "#1f78b4",
                                 "Conventional" = "#e7298a")) +
  labs(x = "Time after intratracheal dose (h)",
       y = "Amount of tobramycin in lungs (ug, log scale)",
       title = "Typical-value profiles: liposomal vs conventional",
       caption = paste("Replicates the comparison in Figure 1 of Marier 2002.",
                       "Liposomal: slower absorption, longer terminal phase.")) +
  theme_minimal(base_size = 11)
```

# Virtual cohort -- between-subject variability

```{r virtual-cohort}
set.seed(2002)
n_sub <- 50L
sim_grid_n <- expand.grid(
  id   = seq_len(n_sub),
  time = c(0, seq(0.1, 24, by = 0.1)),
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)
sim_grid_n$evid <- 0L
sim_grid_n$amt  <- 0
sim_grid_n$cmt  <- "Cc"

dose_rows_n <- data.frame(
  id   = seq_len(n_sub),
  time = 0,
  evid = 1L,
  amt  = 1200,
  cmt  = "depot"
)
events_n <- dplyr::bind_rows(dose_rows_n, sim_grid_n)
stopifnot(!anyDuplicated(unique(events_n[, c("id", "time", "evid")])))

sim_n <- as.data.frame(rxode2::rxSolve(mod, events_n))

sim_summary <- sim_n |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    lo  = stats::quantile(Cc, 0.05, na.rm = TRUE),
    med = stats::median(Cc, na.rm = TRUE),
    hi  = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(sim_summary, aes(time, med)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.25, fill = "#e7298a") +
  geom_line(linewidth = 0.8, colour = "#e7298a") +
  scale_y_log10() +
  labs(x = "Time after intratracheal dose (h)",
       y = "Amount of tobramycin in lungs (ug, log scale)",
       title = "Conventional tobramycin -- virtual cohort (n = 50)",
       caption = paste("Median and 5-95% envelope across 50 simulated rats with",
                       "the published diagonal BSV (CV: ka 38.7%, FL 17.2%,",
                       "k_23 43.9%, k_32 37.8%, k_20 19.8%).")) +
  theme_minimal(base_size = 11)
```

# PKNCA validation against published Table 1

Marier 2002 Table 1 reports three derived NCA-style summaries for the
conventional arm:

* `t_1/2 alpha` = 0.12 h (CV 17.1%)
* `t_1/2 beta`  = 12.9 h (CV 31.4%)
* `AUC`         = 821 ug.h / lungs (CV 15.0%)

```{r pknca-setup}
set.seed(20021)
n_sub_pknca <- 200L
t_grid_pknca <- c(0.1, 0.25, 0.5, 1, 2, 3, 5, 7, 9, 11, 13, 15, 16,
                  18, 24, 30, 36, 48)
sim_grid_pknca <- expand.grid(
  id   = seq_len(n_sub_pknca),
  time = t_grid_pknca,
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)
sim_grid_pknca$evid <- 0L
sim_grid_pknca$amt  <- 0
sim_grid_pknca$cmt  <- "Cc"

dose_rows_pknca <- data.frame(
  id   = seq_len(n_sub_pknca),
  time = 0,
  evid = 1L,
  amt  = 1200,
  cmt  = "depot"
)
events_pknca <- dplyr::bind_rows(dose_rows_pknca, sim_grid_pknca)
stopifnot(!anyDuplicated(unique(events_pknca[, c("id", "time", "evid")])))

sim_pknca <- as.data.frame(rxode2::rxSolve(mod, events_pknca))
```

```{r pknca}
sim_nca <- sim_pknca |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::mutate(treatment = "conventional") |>
  dplyr::select(id, time, Cc, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)

dose_df <- events_pknca |>
  dplyr::filter(evid == 1L) |>
  dplyr::mutate(treatment = "conventional") |>
  dplyr::select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

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

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

nca_df <- as.data.frame(nca_res$result)

nca_medians <- nca_df |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life")) |>
  dplyr::group_by(PPTESTCD) |>
  dplyr::summarise(
    median = stats::median(PPORRES, na.rm = TRUE),
    cv_pct = 100 * stats::sd(PPORRES, na.rm = TRUE) /
                   stats::median(PPORRES, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(nca_medians, digits = 3,
             caption = "Simulated NCA medians (across 200 virtual rats).")
```

## Comparison against published Table 1

```{r pknca-comparison}
get_med <- function(metric) {
  nca_medians$median[nca_medians$PPTESTCD == metric]
}

compare <- data.frame(
  metric      = c("t_1/2 beta (h)", "AUCinf (ug.h)"),
  simulated   = c(get_med("half.life"), get_med("aucinf.obs")),
  published   = c(12.9, 821),
  published_cv_pct = c(31.4, 15.0)
)
compare$pct_diff <- 100 * (compare$simulated - compare$published) /
  compare$published

knitr::kable(compare, digits = c(NA, 2, 1, 1, 1),
             caption = paste("Simulated vs published Table 1 values.",
                             "PKNCA reports the terminal-phase slope only,",
                             "so the t_1/2 alpha comparison is not made here."))
```

The conventional formulation's terminal half-life and AUCinf are
markedly smaller than the liposomal arm, consistent with Table 1's
~8-fold AUC ratio and the paper's primary scientific finding that the
liposomal formulation prolongs pulmonary tobramycin exposure.

# Assumptions and deviations

* **Cc carries amount units, not concentration.** As in the liposomal
  vignette, Marier 2002 fits the amount of tobramycin in lung tissue
  (ug) rather than a concentration. Volumes of distribution were
  intentionally not fitted by the authors; the central compartment in
  this model carries amount in ug, and `Cc` carries amount units, not
  concentration units. PKNCA above treats `Cc` dimensionally as a
  concentration, so the reported `aucinf.obs` and `cmax` are in ug.h
  and ug respectively -- numerically equivalent to the paper's
  "ug.h / lungs" and Fig 1 vertical-axis amounts.
* **Residual error structure.** The combined population residual
  variability reported in Results (18.4% for the conventional arm) is
  encoded as proportional only (`propSd = 0.184`); the dropped
  additive component is documented above. See the liposomal vignette
  for the parallel discussion.
* **NONMEM compartment numbering.** Paper's `k_23` -> canonical `k12`,
  `k_32` -> canonical `k21`, `k_20` -> canonical `kel`. Values
  identical; only the numbering convention differs. See in-`ini()`
  comments for the per-parameter mapping.
* **Random-effect correlations.** Table 1 reports per-parameter CV%
  only and no between-parameter correlations. All etas are encoded as
  uncorrelated diagonal omegas.
* **N = 36 vs n = 39 dosed.** Two animals died during the surgical
  procedure 3 h post-dose and the infection level could not be
  confirmed in a third animal, leaving n = 36 evaluable. The published
  parameter estimates are from the n = 36 final analysis (per Results,
  paragraph 1).
* **Errata search.** No published erratum was located for
  doi:10.1128/aac.46.12.3776-3781.2002 at the time of extraction
  (2026-05-31).
