---
title: "Levofloxacin in rat plasma and prostate ISF (Hurtado 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Levofloxacin in rat plasma and prostate ISF (Hurtado 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

* Citation: `r rxode2::rxode(readModelDb("Hurtado_2014_levofloxacin_rat"))$reference`
* Description: `r rxode2::rxode(readModelDb("Hurtado_2014_levofloxacin_rat"))$description`
* Article: <https://doi.org/10.1128/AAC.01884-13>

# Population

Hurtado et al. (2014) studied 7 male Wistar rats (FEPPS, Porto Alegre; 0.25-0.35
kg body weight) anaesthetised with urethane (1.25 g/kg IP) and instrumented with
a carotid artery cannula for arterial blood sampling and CMA 20 concentric
microdialysis probes (4 mm membrane, 20 kDa cutoff) implanted in the prostate.
Each rat received a single 7 mg/kg IV bolus of levofloxacin into the femoral
vein (calculated to scale from the human 500 mg daily dose for a 70 kg adult).
Plasma was sampled at 0.083, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 6, 8, and 12 h
post-dose; prostate microdialysate was collected at 20-min intervals for up to
12 h. Levofloxacin was quantified in both matrices by HPLC-fluorescence (LOQ 10
ng/mL plasma, 5 ng/mL dialysate). Plasma protein binding measured by in vitro
microdialysis gave fb_plasma = 45.5 +/- 9.4% (so fu_plasma = 0.545); in vivo
prostate microdialysis probe recovery by the no-net-flux method was RR_NNF =
17.5%, used to correct dialysate concentrations to free prostate ISF
concentrations (Hurtado 2014 Materials and Methods, "Animal experiments"
through "Levofloxacin plasma protein binding").

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

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Hurtado_2014_levofloxacin_rat.R`. The table
below collects the references for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| CL (total plasma clearance) | 0.22 L/h | Table 2 |
| V1 (central volume)         | 0.38 L  | Table 2 |
| k12 (central -> peripheral) | 2.27 /h | Table 2 (derived Q = k12 * V1 = 0.8626 L/h) |
| k21 (peripheral -> central) | 1.44 /h | Table 2 (derived Vp = Q / k21 = 0.5990 L) |
| k13 (central -> prostate)   | 0.69 /h | Table 2 |
| k31 (prostate -> central, linear)  | 3.67 /h | Table 2 |
| Vmax (saturable efflux)     | 7.19 ug/h (= 7.19e-3 mg/h) | Table 2 |
| kM (Michaelis constant)     | 0.35 ug/mL (= 0.35 mg/L)   | Table 2 |
| V3* = V_prostate / fu_prostate | 0.05 L | Table 2 |
| Inter-individual variance omega^2(V1)   | 21.0 %CV -> log(1 + 0.210^2) = 0.04316 | Table 2 |
| omega^2(CL)                  | 36.7 %CV -> 0.12640 | Table 2 |
| omega^2(Vmax)                | 41.6 %CV -> 0.15967 | Table 2 |
| omega^2(kM)                  | 76.0 %CV -> 0.45593 | Table 2 |
| Plasma residual: proportional | 10.2 %CV -> propSd  = 0.102  | Table 2 |
| Plasma residual: additive     | 0.085 ug/mL -> addSd = 0.085 | Table 2 |
| Tissue residual: proportional | 15.2 %CV -> propSd_Cprostate = 0.152 | Table 2 |
| Tissue residual: additive     | 0.015 ug/mL -> addSd_Cprostate = 0.015 | Table 2 |
| ODE system (equations 3a-3c) | -- | Results, "Population PK model"; final paragraph references equations 3a-3c with the saturable Vmax/kM term on prostate -> central efflux |
| Residual error structure (combined prop + add, separate per matrix) | -- | Methods, "PK data analysis (ii) Population PK model": "combined error... Separate residual error models were used for total plasma concentrations and free prostate concentrations." |

# Virtual cohort

Individual animal-level data are not publicly available. The figures below use a
virtual cohort that mirrors the Hurtado 2014 design (single 7 mg/kg IV bolus of
levofloxacin in 0.25-0.35 kg male Wistar rats). The VPC overlays use 200
simulated rats; the typical-value replication uses one typical animal.

```{r cohort}
set.seed(20140184L) # AAC.01884-13

n_sim <- 200L

# Per-animal body weight uniform across the source 0.25-0.35 kg window.
bw_kg <- runif(n_sim, min = 0.25, max = 0.35)

cohort <- tibble(
  id      = seq_len(n_sim),
  WT_kg   = bw_kg,
  dose_mg = 7.0 * bw_kg # 7 mg/kg IV bolus per rat
)

# Build the event table. Each rat gets one IV-bolus dose into central at t = 0,
# then observations at the paper's plasma sampling times and at the midpoints
# of the 20-min prostate microdialysis collection intervals.
plasma_times <- c(0.083, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 6, 8, 12)
prostate_mid <- seq(from = 10/60, by = 20/60, length.out = 36) # 0.167, 0.500, ..., 11.83 h

obs_times <- sort(unique(c(plasma_times, prostate_mid)))

events <- cohort |>
  group_by(id) |>
  reframe(
    WT_kg   = WT_kg,
    dose_mg = dose_mg,
    rec = list(
      bind_rows(
        # IV bolus into central
        tibble(time = 0, evid = 1L, amt = unique(dose_mg), cmt = "central",
               DV = NA_real_),
        # observation rows -- rxode2 returns Cc and Cprostate on every row
        tibble(time = obs_times, evid = 0L, amt = NA_real_, cmt = "Cc",
               DV = NA_real_)
      )
    )
  ) |>
  tidyr::unnest(rec)

# Sanity: ids are unique within this single-cohort design.
stopifnot(!anyDuplicated(events[, c("id", "time", "evid")]))
```

# Simulation

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

# Stochastic simulation (full IIV + residual error) for the VPC overlay.
sim <- rxode2::rxSolve(mod, events, keep = c("WT_kg")) |>
  as.data.frame()
```

For the typical-value replication of Figure 1 (one representative animal,
fixed effects only), suppress between-subject variability and residual error:

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

cohort_typical <- tibble(
  id      = 1L,
  WT_kg   = 0.30,                 # midpoint of the 0.25-0.35 kg range
  dose_mg = 7.0 * 0.30            # 2.10 mg total dose -- matches the paper's mean dose 2.198 mg
)

events_typical <- cohort_typical |>
  group_by(id) |>
  reframe(
    WT_kg   = WT_kg,
    dose_mg = dose_mg,
    rec = list(
      bind_rows(
        tibble(time = 0,                              evid = 1L, amt = unique(dose_mg), cmt = "central", DV = NA_real_),
        tibble(time = seq(0, 12, length.out = 241L), evid = 0L, amt = NA_real_,        cmt = "Cc",       DV = NA_real_)
      )
    )
  ) |>
  tidyr::unnest(rec)

sim_typical <- rxode2::rxSolve(mod_typical, events_typical, keep = c("WT_kg")) |>
  as.data.frame()
```

# Replicate published figures

```{r figure-1-typical}
# Replicates the structure of Hurtado 2014 Figure 1 (typical-value population
# fits, dashed line PRED in the paper). The left panel shows total plasma
# levofloxacin; the right panel shows the free prostate ISF concentration.
typical_long <- sim_typical |>
  pivot_longer(c(Cc, Cprostate), names_to = "matrix", values_to = "conc") |>
  filter(!is.na(conc), time > 0)

ggplot(typical_long, aes(time, conc)) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~matrix, scales = "free_y",
             labeller = labeller(matrix = c(
               Cc        = "Total plasma (ug/mL)",
               Cprostate = "Free prostate ISF (ug/mL)"))) +
  scale_y_log10() +
  labs(
    x = "Time after 7 mg/kg IV bolus (h)",
    y = "Levofloxacin concentration (ug/mL, log scale)",
    title = "Figure 1 -- typical-value population fits",
    caption = "Replicates the PRED (dashed) line in Hurtado 2014 Figure 1."
  )
```

```{r figure-3-vpc-plasma}
# Replicates Hurtado 2014 Figure 3 (plasma panel): VPC of total plasma
# levofloxacin with median and 90% prediction interval over 1000 simulated
# rats (this vignette uses 200 to keep the render fast).
plasma_vpc <- sim |>
  filter(!is.na(Cc), time > 0) |>
  group_by(time) |>
  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(plasma_vpc, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(linewidth = 0.8, color = "steelblue") +
  scale_y_log10() +
  labs(
    x = "Time after 7 mg/kg IV bolus (h)",
    y = "Total plasma levofloxacin (ug/mL, log scale)",
    title = "Figure 3 -- plasma VPC",
    caption = "Replicates the plasma panel of Hurtado 2014 Figure 3 (median + 90% PI)."
  )
```

```{r figure-3-vpc-prostate}
# Replicates Hurtado 2014 Figure 3 (tissue panel): VPC of free prostate ISF
# concentration. The saturable efflux dominates at high prostate concentrations
# during the distribution phase (0-2 h) and the linear k31 efflux dominates
# during the elimination phase, producing the characteristic faster apparent
# decline in tissue than in plasma (paper terminal t1/2: plasma 5.0 h vs
# prostate 2.3 h).
prostate_vpc <- sim |>
  filter(!is.na(Cprostate), time > 0) |>
  group_by(time) |>
  summarise(
    Q05 = quantile(Cprostate, 0.05, na.rm = TRUE),
    Q50 = quantile(Cprostate, 0.50, na.rm = TRUE),
    Q95 = quantile(Cprostate, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(prostate_vpc, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "tomato") +
  geom_line(linewidth = 0.8, color = "tomato") +
  scale_y_log10() +
  labs(
    x = "Time after 7 mg/kg IV bolus (h)",
    y = "Free prostate ISF levofloxacin (ug/mL, log scale)",
    title = "Figure 3 -- prostate ISF VPC",
    caption = "Replicates the tissue panel of Hurtado 2014 Figure 3 (median + 90% PI)."
  )
```

# PKNCA validation

PKNCA is run on the stochastic simulation, one block per output (total plasma
`Cc` and free prostate ISF `Cprostate`). All animals are in a single dose group
(7 mg/kg IV bolus), so the grouping variable `treatment = "IV_7mgkg"` is set
to a constant column.

```{r pknca-plasma}
# Deduplicate so PKNCA sees one (id, time) row per subject (rxode2 returns
# Cc and Cprostate on every observation row regardless of which cmt the
# row requested).
plasma_nca <- sim |>
  filter(!is.na(Cc), time > 0) |>
  distinct(id, time, .keep_all = TRUE) |>
  mutate(treatment = "IV_7mgkg") |>
  select(id, time, Cc, treatment)

dose_df <- events |>
  filter(evid == 1, cmt == "central") |>
  mutate(treatment = "IV_7mgkg") |>
  select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(plasma_nca, Cc ~ time | treatment + id,
                             concu = "ug/mL", 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
)

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

plasma_summary <- as.data.frame(nca_plasma$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life")) |>
  group_by(PPTESTCD) |>
  summarise(
    median_value = median(PPORRES, na.rm = TRUE),
    sd_value     = sd(PPORRES,     na.rm = TRUE),
    .groups      = "drop"
  )

knitr::kable(
  plasma_summary,
  caption = "Simulated plasma NCA (n = 200 rats, 7 mg/kg IV). Units: Cmax ug/mL, Tmax h, AUCinf.obs ug/mL*h, t1/2 h."
)
```

```{r pknca-prostate}
prostate_nca <- sim |>
  filter(!is.na(Cprostate), time > 0) |>
  distinct(id, time, .keep_all = TRUE) |>
  mutate(treatment = "IV_7mgkg") |>
  select(id, time, Cprostate, treatment) |>
  rename(Cc = Cprostate) # PKNCAconc needs a single conc column name

conc_obj_p <- PKNCA::PKNCAconc(prostate_nca, Cc ~ time | treatment + id,
                               concu = "ug/mL", timeu = "h")

nca_prostate <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_p, dose_obj, intervals = intervals)
)

prostate_summary <- as.data.frame(nca_prostate$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life")) |>
  group_by(PPTESTCD) |>
  summarise(
    median_value = median(PPORRES, na.rm = TRUE),
    sd_value     = sd(PPORRES,     na.rm = TRUE),
    .groups      = "drop"
  )

knitr::kable(
  prostate_summary,
  caption = "Simulated free-prostate-ISF NCA (n = 200 rats, 7 mg/kg IV). Units: Cmax ug/mL, Tmax h, AUCinf.obs ug/mL*h, t1/2 h."
)
```

## Comparison against published NCA (Table 1)

Hurtado 2014 Table 1 reports per-rat (n = 7) NCA means +/- SD. Free plasma
values in the paper are derived from total plasma multiplied by
fu_plasma = 0.545 (so e.g. free-plasma Cmax = 2.52 = 5.74 * 0.455 -- note the
paper's free fraction is reported as 1 - fb_plasma = 0.545 but the C0 entry
implies a slightly different rounding). The packaged model simulates **total
plasma** Cc and **free prostate ISF** Cprostate, so the natural comparison is:

| Parameter | Source | Published (Table 1, mean +/- SD, n = 7) | Simulated median (n = 200) |
|-----------|--------|-----------------------------------------|----------------------------|
| Plasma Cmax (ug/mL)     | total      | 4.58 +/- 0.57  | (see Cmax in plasma table above) |
| Plasma C0 (ug/mL)       | total      | 5.74 +/- 0.82  | dose / V1 = 2.10 / 0.38 = 5.53 ug/mL (typical) |
| Plasma AUC0-inf (ug/mL*h) | total    | 11.7 +/- 4.3   | dose / CL = 2.10 / 0.22 = 9.55 ug/mL*h (typical); paper text computes 9.9 from V1 = 0.38, CL = 0.22 |
| Plasma t1/2 (h)         | total      | 5.0 +/- 1.7    | (see half.life in plasma table above) |
| Plasma MRT (h)          | total      | 6.1 +/- 2.7    | not computed by PKNCA in default intervals |
| Plasma CL (L/h)         | total      | 0.21 +/- 0.08  | model lcl = log(0.22) -> 0.22 L/h |
| Plasma Vss (L)          | total      | 1.2 +/- 0.4    | Vc + Vp = 0.38 + 0.599 = 0.979 L (typical) |
| Prostate Cmax (ug/mL)   | free       | 2.31 +/- 0.65  | (see Cmax in prostate table above) |
| Prostate AUC0-inf (ug/mL*h) | free   | 4.8  +/- 1.9   | (see aucinf.obs in prostate table above) |
| Prostate t1/2 (h)       | free       | 2.3  +/- 0.7   | (see half.life in prostate table above) |

The simulated typical-value plasma AUC of 9.55 ug/mL*h is ~18% below the NCA
mean (11.7 ug/mL*h), within the typical bias of comparing a population-model
typical value against an NCA arithmetic mean computed over a small (n = 7)
sample with appreciable inter-rat variability (CL %CV = 37%, Vmax %CV = 42%).
Hurtado 2014 Discussion makes the same back-calculation: "From the relationship
AUC0-inf = dose/CL, the plasma AUC0-inf was calculated to be 9.9 ug/mL*h, being
in close agreement with the results presented in Table 1." The simulated
prostate Cmax / AUC / t1/2 should likewise land within ~20% of the Table 1
free-prostate values, with the simulated prostate t1/2 noticeably shorter than
plasma (the model's k31 + saturable efflux drives a faster apparent decline in
prostate, matching the paper's narrative).

# Assumptions and deviations

* **Asymmetric prostate transport requires lk13 and lk31 as primary
  parameters.** The Hurtado 2014 NONMEM 6 ADVAN6 TRANS1 model parameterises
  the central -> peripheral arm using rate constants k12 / k21 (algebraically
  equivalent to a single Q with derived Vp, which the model file does compute
  in canonical lcl / lvc / lq / lvp form), but the central -> prostate arm is
  genuinely asymmetric: k13 (uptake) and k31 (linear efflux) plus a saturable
  Michaelis-Menten efflux do not collapse to a single inter-compartmental
  clearance. The model file therefore keeps lk13 and lk31 as primary
  log-transformed parameters (precedent: `Clinckers_2008_MHD_rat.R`
  lk23/lk32). `checkModelConventions()` may flag these as outside the
  canonical lcl / lq / lvp set; the deviation is intentional and is the
  cheapest faithful encoding of the paper's structure.

* **Effect compartment used for the prostate ISF.** The third compartment
  carries the prostate interstitial-fluid free concentration. The canonical
  `effect` compartment role (PK/PD biophase) is used here for a microdialysis-
  measured tissue site, following the precedent set by
  `Clinckers_2008_MHD_rat.R` (hippocampal ECF, also microdialysis). The
  output variable name `Cprostate` makes the anatomical site explicit; the
  underlying compartment is `effect`.

* **Prostate apparent volume V3*.** Because the microdialysis probe measures
  unbound prostate concentrations, the paper can only estimate the apparent
  volume V3* = V_prostate / fu_prostate, not V_prostate directly. The model
  file's `lveff` is this apparent volume (0.05 L); a downstream user wanting
  to recover the physical prostate volume would need an independent estimate
  of fu_prostate (the paper hypothesises fu_prostate ~ 1 because Hurtado
  Discussion notes "the volume of the prostate tissue compartment (V3) was
  found to be 8-fold smaller than V1", consistent with the physical prostate
  volume in a 0.3 kg rat).

* **Vmax unit conversion.** Hurtado 2014 Table 2 reports Vmax = 7.19 ug/h
  while the model's internal amount unit is mg (consistent with the input
  dose in mg and volumes in L producing concentrations in mg/L = ug/mL).
  `lvmax` is therefore set to `log(7.19e-3)` so that the saturable-efflux
  ODE term `vmax * Cprostate / (km + Cprostate)` has the correct mg/h units.
  kM is unit-invariant (0.35 ug/mL = 0.35 mg/L numerically).

* **CV% -> internal variance for IIV.** Table 2 reports IIV in %CV; the
  model file converts each value to the internal lognormal-variance scale
  via omega^2 = log(1 + CV^2). The annotations in `ini()` show the per-
  parameter arithmetic.

* **Residual error magnitude in %CV / ug/mL on the linear scale.** Table 2
  labels the four residual-error rows as sigma^2 with units %CV
  (proportional) or ug/mL (additive). The values are interpreted as
  residual standard deviations on the linear concentration scale (matching
  Hurtado's prose: "different residual error models (e.g., constant
  coefficient of variation [CV], combined error) were tested"). propSd =
  0.102 corresponds to the 10.2 %CV plasma proportional term, etc.

* **No covariates.** Hurtado 2014 did not investigate covariate effects --
  the n = 7 design is too small. The model file's `covariateData` list is
  empty; weight is carried as a per-animal scalar only to compute the per-
  animal mg dose from the 7 mg/kg dose level.

* **Single-dose IV bolus only.** The original study used a single 7 mg/kg
  IV bolus. The packaged model can be re-purposed to simulate multiple-
  dose regimens (the Hurtado 2014 Discussion explicitly invites this) by
  supplying additional `evid = 1` dose records into the `central`
  compartment in the event table.
