---
title: "Maraviroc IV (Weatherley 2009)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Maraviroc IV (Weatherley 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(tibble)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Weatherley_2009_maraviroc_iv"))$reference`
* Description: `r rxode2::rxode(readModelDb("Weatherley_2009_maraviroc_iv"))$description`
* Article: <https://doi.org/10.1111/j.1365-2125.2009.03455.x>

Weatherley & McFadyen 2009 (Br J Clin Pharmacol 68(5):648-657) report a
three-part modelling strategy for the CCR5 entry inhibitor maraviroc: (1) a
four-compartment population PK model of intravenous maraviroc fit to data from
20 healthy males in study A4001009; (2) a sigmoid Emax meta-regression of
dose-normalized NCA AUC pooled across five oral phase 1 studies (n = 134) to
describe the dose/food/formulation absorption non-linearity; and (3) a
closed-form mass balance model in S-PLUS that combines the IV disposition
with the oral absorption prediction and the published mass-balance study to
apportion clearance into renal / hepatic-metabolic / hepatic-secretion
components and derive a hepatic blood flow (101 L/h) for use downstream.

Only the first analysis is an ODE-form population PK model and is the focus
of this package entry. The second analysis is an algebraic dose-AUC
relationship (one observation per subject per dose, no time dimension) and
the third is a deterministic mass-balance calculation, neither of which fit
the `nlmixr2lib` model file shape. See *Assumptions and deviations* for
details and the cited supporting equations.

# Population

The intravenous analysis used 473 plasma concentrations from 20 healthy adult
male volunteers (ages 22 - 42 years, weights 65 - 93 kg) in Pfizer study
A4001009 (Weatherley 2009 Table 1, row "A4001009 / Population IV only"). Each
subject received 3, 10, or 30 mg of maraviroc as a one-hour intravenous
infusion in solution after an overnight fast. Maraviroc plasma concentrations
were measured by solid-phase extraction followed by LC-MS/MS detection
(LLOQ 0.5 ng/mL). No covariate effects were tested in the IV-only analysis
beyond a test of dose effect on clearance, which was found to be absent
across the 3 - 30 mg range.

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

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl`  (CL)          | log(48.0) L/h  | Table 3 row "CL (l h-1)"  (2.2% SE) |
| `lvc`  (V1)          | log(13.5) L    | Table 3 row "V1 (l)"     (14.4% SE) |
| `lvp`  (V2)          | log(16.1) L    | Table 3 row "V2 (l)"     (10.5% SE) |
| `lvp2` (V3)          | log(140) L     | Table 3 row "V3 (l)"     (5.8% SE) |
| `lvp3` (V4)          | log(36.4) L    | Table 3 row "V4 (l)"     (4.0% SE) |
| `lq`  (Q2)           | log(35.4) L/h  | Table 3 row "Q2 (l h-1)" (15.3% SE) |
| `lq2` (Q3)           | log(9.13) L/h  | Table 3 row "Q3 (l h-1)" (4.8% SE) |
| `lq3` (Q4)           | log(16.8) L/h  | Table 3 row "Q4 (l h-1)" (8.9% SE) |
| IIV variances        | log(CV^2 + 1)  | Standard log-normal mapping from Table 3 CV% values |
| IIV CL (hCL = 10.8%) | omega^2 = 0.01160 | Table 3 row "hCL (%)" |
| IIV V1 (hV1 = 25.6%) | omega^2 = 0.06348 | Table 3 row "hV1 (%)" |
| IIV V2 (hV2 = 18.0%) | omega^2 = 0.03189 | Table 3 row "hV2 (%)" |
| IIV Q3 (hQ3 = 18.8%) | omega^2 = 0.03473 | Table 3 row "hQ3 (%)" |
| IIV V3 (hV3 = 16.9%) | omega^2 = 0.02816 | Table 3 row "hV3 (%)" |
| `propSd`             | 0.109 (fraction)  | Table 3 row "Residual (%)" 10.9% (NONMEM EPS on log-transformed Y; equivalent to proportional in linear space) |
| 4-compartment ODE    | n/a               | Methods p.650 ("Two-, three- and four-compartment models were fitted... using the general linear compartmental NONMEM subroutine ADVAN7"); four-compartment model selected (Results p.651, "106 points lower in objective function value than the three-compartment model") |
| FOCE-I               | n/a               | Methods p.650 ("The FOCE method with INTERACTION was used") |
| Log-transformed data | n/a               | Methods p.650 ("Iv plasma concentration vs time data were logarithmically transformed prior to analysis") |
| Phase half-lives     | 0.08, 0.49, 2.03, 12.9 h | Results p.651 (empirical Bayesian estimates of individual parameter values, 4-phase decline) |
| Phase AUC fractions  | 36%, 23%, 23%, 18% | Results p.651 (percentages of the area under each phase) |

# Virtual cohort

We simulate the three IV cohorts from study A4001009 (3 mg, 10 mg, 30 mg
single-dose, 1 h infusion in solution) using 60 virtual subjects per dose
level (180 subjects total). Subject IDs are made disjoint across the three
cohorts so that bind_rows-ed data sets do not silently merge into
"Frankenstein" subjects in rxSolve.

```{r cohort, message = FALSE, warning = FALSE}
set.seed(20260613L)

mod_full <- rxode2::rxode(readModelDb("Weatherley_2009_maraviroc_iv"))
mod_typ  <- rxode2::zeroRe(mod_full)

infusion_h <- 1
obs_times <- sort(unique(c(
  seq(0, 1, by = 0.05),
  seq(1, 6, by = 0.1),
  seq(6, 48, by = 0.5)
)))

make_cohort <- function(n, dose_mg, id_offset = 0L) {
  obs <- tibble::tibble(
    id    = id_offset + rep(seq_len(n), each = length(obs_times)),
    time  = rep(obs_times, times = n),
    amt   = 0,
    rate  = 0,
    evid  = 0L,
    cmt   = NA_character_,
    dose_mg = dose_mg,
    cohort  = sprintf("%g mg IV (1 h)", dose_mg)
  )
  dose <- tibble::tibble(
    id    = id_offset + seq_len(n),
    time  = 0,
    amt   = dose_mg,
    rate  = dose_mg / infusion_h,
    evid  = 1L,
    cmt   = "central",
    dose_mg = dose_mg,
    cohort  = sprintf("%g mg IV (1 h)", dose_mg)
  )
  dplyr::bind_rows(dose, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(
  make_cohort(60L,  3, id_offset =   0L),
  make_cohort(60L, 10, id_offset = 100L),
  make_cohort(60L, 30, id_offset = 200L)
)

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

# Simulation

```{r simulate}
sim <- rxode2::rxSolve(mod_full, events, keep = c("cohort", "dose_mg")) |>
  as.data.frame()

# Typical (zero-eta) profile for line overlay
events_typ <- dplyr::bind_rows(
  make_cohort(1L,  3, id_offset = 1000L),
  make_cohort(1L, 10, id_offset = 1100L),
  make_cohort(1L, 30, id_offset = 1200L)
)
sim_typ <- rxode2::rxSolve(mod_typ, events_typ, keep = c("cohort", "dose_mg")) |>
  as.data.frame()
```

# Replicate published figure

Weatherley 2009 Figure 5 is a predictive check for the 30 mg IV dose
(median + 95% CI of 100 simulated profiles plus the original observations).
We reproduce the same VPC shape from our virtual cohort. Observed
concentrations are not redistributed with the article and cannot be overlaid.

```{r figure-5}
vpc_df <- sim |>
  dplyr::filter(time > 0, !is.na(Cc)) |>
  dplyr::group_by(cohort, time) |>
  dplyr::summarise(
    Q025 = quantile(Cc, 0.025, na.rm = TRUE),
    Q500 = quantile(Cc, 0.500, na.rm = TRUE),
    Q975 = quantile(Cc, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

vpc_df |>
  dplyr::filter(cohort == "30 mg IV (1 h)") |>
  ggplot(aes(time, Q500)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue", linewidth = 0.7) +
  scale_x_continuous(limits = c(0, 48)) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Maraviroc plasma concentration (ng/mL)",
       title = "Replicates Figure 5 of Weatherley 2009 (30 mg IV, 1 h infusion)",
       caption = "Median and 95% prediction interval from 60 virtual subjects.")
```

All three cohorts on a single panel, log-linear y-axis:

```{r figure-all-doses}
vpc_df |>
  ggplot(aes(time, Q500, colour = cohort, fill = cohort)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.7) +
  scale_x_continuous(limits = c(0, 48)) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Maraviroc plasma concentration (ng/mL)",
       title = "Predicted IV maraviroc profiles across the three A4001009 dose levels",
       caption = "Median and 95% prediction interval per dose, 60 virtual subjects each.")
```

# PKNCA validation

PKNCA is used to compute Cmax, Tmax, AUC0-inf, and terminal half-life for
each cohort on the typical (zero-eta) profile. The paper does not print an
NCA table per dose for the IV cohort -- the structural model parameters and
the four-phase half-life decomposition (0.08, 0.49, 2.03, 12.9 h with AUC
fractions 36%, 23%, 23%, 18%; Weatherley 2009 Results p.651) are the published
checkpoints. Dose-linear AUC (AUC0-inf = dose / CL = dose / 48 L/h) gives
62.5, 208, and 625 ng.h/mL for 3, 10, and 30 mg respectively, and the
PKNCA result reproduces these within rounding.

```{r pknca, message = FALSE, warning = FALSE}
sim_nca <- sim_typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort)

sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, cohort) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, cohort, time, .keep_all = TRUE) |>
  dplyr::arrange(id, cohort, time)

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

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | cohort + id,
                             concu = "ng/mL", timeu = "h")
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
)

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

## Comparison against published values

The paper reports total clearance CL = 48 L/h (Table 3) and a terminal
elimination half-life of 12.9 h (Results p.651). AUC0-inf for a single
IV dose is therefore dose / CL by linear PK. Tmax is at the end of infusion
(1 h).

```{r pknca-compare}
nca_df <- as.data.frame(nca_res$result) |>
  dplyr::select(cohort, PPTESTCD, PPORRES) |>
  dplyr::group_by(cohort, PPTESTCD) |>
  dplyr::summarise(PPORRES = dplyr::first(PPORRES), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

published <- tibble::tribble(
  ~cohort,           ~cmax,  ~tmax, ~aucinf.obs, ~half.life,
  "3 mg IV (1 h)",   NA_real_, 1.0,  62.5,  12.9,
  "10 mg IV (1 h)",  NA_real_, 1.0, 208.3,  12.9,
  "30 mg IV (1 h)",  NA_real_, 1.0, 625.0,  12.9
)

cmp <- nca_df |>
  dplyr::left_join(
    published |> dplyr::rename_with(~ paste0(.x, "_ref"), -cohort),
    by = "cohort"
  ) |>
  dplyr::mutate(
    `AUC0-inf sim`   = sprintf("%.1f", aucinf.obs),
    `AUC0-inf ref`   = sprintf("%.1f", aucinf.obs_ref),
    `Tmax sim`       = sprintf("%.2f", tmax),
    `Tmax ref`       = sprintf("%.2f", tmax_ref),
    `t1/2 sim`       = sprintf("%.1f", half.life),
    `t1/2 ref`       = sprintf("%.1f", half.life_ref)
  ) |>
  dplyr::select(cohort,
                `AUC0-inf sim`, `AUC0-inf ref`,
                `Tmax sim`, `Tmax ref`,
                `t1/2 sim`, `t1/2 ref`)

knitr::kable(
  cmp,
  caption = paste(
    "Simulated NCA vs. Weatherley 2009 reference values.",
    "Reference AUC0-inf computed as dose / CL with CL = 48 L/h (Table 3);",
    "reference Tmax = end of infusion;",
    "reference terminal half-life = 12.9 h (Results p.651)."
  ),
  align = c("l", rep("r", 6))
)
```

# Assumptions and deviations

- **Companion analyses 2 and 3 not packaged as ODE models.** Weatherley 2009
  describes three distinct analyses (the IV popPK, the oral sigmoid-Emax
  NAUC meta-regression, and the closed-form S-PLUS mass balance). Only the
  IV popPK is structurally an ODE model and is shipped here. The sigmoid
  Emax model (Methods p.650, Equation 1) operates on dose-normalized NCA AUC
  (one observation per subject per dose, no time dimension) and would need a
  static `~ predict(dose, food, formulation, dosing)` shape that the
  `nlmixr2lib` `ini()` / `model()` template is not designed for. The mass
  balance calculations (Appendix B1 and B2) are closed-form arithmetic over
  the IV CL, the published mass-balance recovery (Table 2), and the absorbed
  fraction derived from the NAUC model; they yield the published hepatic
  extraction E = 0.602 and hepatic blood flow QH,B = 101 L/h, but they are
  not an ODE either. Furthermore, the paper reports only partial parameter
  values for the NAUC model in narrative form (ED50 fasted = 59.7 mg,
  ED50 fed = 365 mg; full Hill coefficient, NAUC0 and NAUCmax estimates,
  and the per-covariate effect estimates are not tabulated in the published
  article), so faithful reproduction of analysis 2 from the on-disk source
  is not possible without additional information.

- **No covariates in this model.** Weatherley 2009 explicitly tested only a
  dose effect on CL across the 3 - 30 mg IV range (Methods p.650, Results
  p.651) and found none. No body weight, age, sex, or renal-function effect
  was modelled because the cohort (20 healthy males, ages 22 - 42, weights
  65 - 93 kg) lacks the spread to identify such effects from the IV-only
  data set. `covariateData = list()` reflects the paper's reporting.

- **Residual error encoded as `propSd = 0.109`.** The paper fits residual
  variability with NONMEM EPS on the log-transformed dependent variable
  (Methods p.650, "proportional residual variability was modelled using
  additive form on log transformed data"). The Table 3 value of 10.9% is
  the on-log-scale CV%, which equals the linear-scale proportional residual
  to first order. Encoding as `Cc ~ prop(propSd)` with `propSd = 0.109`
  matches the paper to that approximation. For the 10.9% level the
  difference between exp-normal (`~ lnorm(0.109)`) and proportional
  (`~ prop(0.109)`) on the predicted concentration is negligible (<1%).

- **IIV declared only where reported.** Table 3 reports IIV on CL, V1, V2,
  Q3, and V3 but not on V4, Q2, or Q4. The model file declares etas only on
  the five reported parameters; the other three are taken as fixed-effect
  population values without between-subject variability, matching the
  paper's NONMEM `$OMEGA` structure.

- **Infusion duration is supplied by the data, not the model.** The source
  study used a 1 h IV infusion. The model code does not impose a structural
  infusion duration; the simulation in this vignette supplies `rate` on the
  dosing record to reproduce the 1 h infusion explicitly.

- **No NCA table per cohort published in the paper.** The published
  validation in Weatherley 2009 Figure 5 is a visual predictive check, not
  an NCA-table comparison. The reference column in the PKNCA comparison
  table above uses the published total clearance (48 L/h) and terminal
  half-life (12.9 h) to derive expected AUC0-inf = dose / CL by linear
  PK; the reference Tmax is the end of infusion. Cmax was not pinned by a
  published value and is omitted from the reference column.
