---
title: "Tobramycin liposomal intratracheal (Marier 2002) -- rat"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tobramycin liposomal 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 **liposomal** formulation arm. The
**conventional** formulation arm is in the sibling vignette
[`Marier_2002_tobramycin_rat_conventional`](Marier_2002_tobramycin_rat_conventional.html);
they share the same ODE structure and differ only in their fitted
parameter values.

# Population

The model was developed from 39 male Sprague-Dawley rats (weight
175-225 g; Charles River, Saint-Constant, Quebec) with established
chronic pulmonary *B. cepacia* infection (strain BC 1368, genomovar
III, tobramycin MIC 128 ug/mL). Infection was established 6 days
pre-dose by intratracheal inoculation of agar beads containing 10^6
CFU; pre-dose throat-swab cultures confirmed infection in every
animal.

A single 1,200 ug intratracheal dose of liposomal tobramycin
(DPPC:DMPG 10:1 phospholipids, extruded to 230-400 nm with Tobi
rehydration) was administered as 100 uL via calibrated pipette
followed by an air bolus. Three rats per time point were
cardiac-puncture exsanguinated under ketamine + xylazine anaesthesia
at 0.5, 1, 2, 3, 5, 7, 9, 11, 13, 15, 16, 18, and 24 h post-dose
(13 sampling times x 3 rats = 39 animals total). HPLC-UV quantitated
total tobramycin in homogenised lung tissue (LOQ ~10 ug per lungs);
the assay reports the sum of free and liposome-encapsulated
tobramycin.

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

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| ka (`lka`)           | 2.30 /h (CV 30.8%) | Table 1 (Liposomal column) |
| FL (`lfdepot`)       | 0.85 (CV 12.2%) | Table 1 |
| k_23 (`lk12`)        | 1.19 /h (CV 38.5%) | Table 1 |
| k_32 (`lk21`)        | 0.407 /h (CV 40.2%) | Table 1 |
| k_20 (`lkel`)        | 0.155 /h (CV 30.2%) | Table 1 |
| IIV (diagonal CV%)   | as in `ini()` comments | Table 1 (CV column) |
| Residual variability | 18.1% (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.40 h (CV 11.0%) | Table 1 (analytic, from individual Bayesian estimates) |
| Derived t_1/2 beta   | 19.7 h (CV 24.2%) | Table 1 |
| Derived AUC          | 6,811 ug.h / lungs (CV 19.7%) | Table 1 |

# Loading the model

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

# Replicating Figure 1 (liposomal 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 liposomal panel shows a slow
absorption and a long terminal phase. The chunk below builds an event
table for a single virtual rat, simulates the typical-value (no-IIV)
trajectory, and overlays the published mean values reported in the
text and Fig 1 caption.

```{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 = "#1f78b4") +
  scale_y_log10(limits = c(10, 2000)) +
  labs(x = "Time after intratracheal dose (h)",
       y = "Amount of tobramycin in lungs (ug, log scale)",
       title = "Liposomal tobramycin -- typical-value profile",
       caption = paste("Replicates the liposomal panel of Figure 1 of Marier 2002.",
                       "The terminal slope corresponds to the published t_1/2 beta =",
                       "19.7 h.")) +
  theme_minimal(base_size = 11)
```

# Virtual cohort -- between-subject variability

A virtual cohort of 50 rats with the same 1,200 ug intratracheal dose
illustrates the between-subject variability implied by the diagonal
IIV terms in Table 1. The shaded ribbon is the 5-95% prediction
envelope; the solid line is the median trajectory.

```{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 = "#1f78b4") +
  geom_line(linewidth = 0.8, colour = "#1f78b4") +
  scale_y_log10() +
  labs(x = "Time after intratracheal dose (h)",
       y = "Amount of tobramycin in lungs (ug, log scale)",
       title = "Liposomal tobramycin -- virtual cohort (n = 50)",
       caption = paste("Median and 5-95% envelope across 50 simulated rats with",
                       "the published diagonal BSV (CV: ka 30.8%, FL 12.2%,",
                       "k_23 38.5%, k_32 40.2%, k_20 30.2%).")) +
  theme_minimal(base_size = 11)
```

# PKNCA validation against published Table 1

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

* `t_1/2 alpha` = 0.40 h (CV 11.0%) -- first-phase half-life
* `t_1/2 beta`  = 19.7 h (CV 24.2%) -- terminal-phase half-life
* `AUC`         = 6,811 ug.h / lungs (CV 19.7%)

PKNCA computes Cmax, Tmax, AUCinf, and a single terminal half-life
(corresponding to t_1/2 beta). The cohort below uses a virtual
population (n = 200, dense sampling) so that the per-subject NCA
medians are stable. The amount-in-lung values from the model are
passed as the `Cc` column; PKNCA treats them dimensionally as
concentration, so the reported `aucinf.obs` is in ug.h units and is
directly comparable to the paper's "ug.h / lungs" value.

```{r pknca-setup}
set.seed(20020)
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, 60, 72)
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 = "liposomal") |>
  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 = "liposomal") |>
  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(19.7, 6811),
  published_cv_pct = c(24.2, 19.7)
)
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.",
                             "t_1/2 alpha is the rapid distribution half-life;",
                             "PKNCA reports the terminal slope only, so the",
                             "alpha-phase comparison is not made here."))
```

PKNCA estimates a single terminal half-life from log-linear regression
on the last few sampling points (default lambda-z selection); this
corresponds to the paper's t_1/2 beta. Both the terminal half-life and
AUCinf reproduce the published Table 1 medians to within ~10%, which
is the expected agreement given that the published values are derived
analytically from individual Bayesian POSTHOC estimates whereas the
PKNCA values are computed from sparse-grid simulated concentrations.

# Assumptions and deviations

* **Cc carries amount units, not concentration.** Marier 2002 fits the
  amount of tobramycin in lung tissue (ug, computed by the assay as
  homogenate concentration times lung volume per animal) rather than a
  concentration. Volumes of distribution were intentionally not fitted
  by the authors (Methods, "Pharmacokinetic analysis": "On the other
  hand, volumes of distribution were not fitted because of this."), so
  the central compartment in this model carries amount in ug, and the
  observation variable `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.** Marier 2002 Methods declare a combined
  proportional + additive residual error model, but the Results
  paragraph after Table 1 reports only a single combined population
  residual variability (18.1% for the liposomal arm). The individual
  proportional and additive SDs are not separately reported. The model
  file encodes the residual as proportional only (`propSd = 0.181`); a
  non-zero additive component was either dropped during model
  selection or estimated as negligible and not reported. Downstream
  users requiring a separate additive component should consult the
  original NONMEM control stream (not on disk for this extraction).
* **NONMEM compartment numbering.** The paper's `k_23`, `k_32`, `k_20`
  labels follow NONMEM ADVAN4 numbering (depot=1, central=2,
  peripheral=3). The canonical nlmixr2 micro-constant names map as:
  paper `k_23` -> canonical `k12` (central to peripheral1); paper
  `k_32` -> canonical `k21` (peripheral1 to central); paper `k_20` ->
  canonical `kel` (elimination from central). The numerical values are
  the same; only the numbering convention differs. See the
  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; this matches the most plausible
  diagonal-OMEGA reading of the paper's parameter table.
* **Encapsulation status of the observation.** The HPLC assay measures
  total tobramycin in the lung tissue homogenate -- the sum of
  liposome-encapsulated and free tobramycin -- per Methods
  ("Analytical assay": "This analytical method for the liposomal
  formulation gives pulmonary tobramycin concentrations that represent
  the summation of the concentration achieved with the encapsulated
  form and the concentration achieved with the free form."). The model
  output therefore does not distinguish encapsulated from free
  tobramycin; downstream PK-PD work that needs the free-drug fraction
  must add a separate sub-model.
* **Lung volume per animal is implicit.** The amounts the paper reports
  for individual rats were calculated by multiplying the measured
  homogenate concentration by the per-animal lung volume. The model
  output is per-animal amount, so any downstream "concentration" or
  "AUC per gram of lung" calculation requires the user to specify the
  lung mass / volume for their virtual cohort.
* **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). If a future correction surfaces, the affected `ini()`
  comments and the `reference` field should be updated in a follow-up
  PR.
