---
title: "Raltegravir (Wang 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Raltegravir (Wang 2011)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

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

Wang 2011 was a pilot pharmacokinetic study of a single 400 mg oral
raltegravir dose in six healthy male Singaporean volunteers, designed
to determine whether sustained antiviral efficacy in plasma-PK-variable
HIV patients could be explained by intracellular drug accumulation in
peripheral blood mononuclear cells (PBMCs). The population PK model
links plasma and intracellular concentrations through a simple
empirical accumulation ratio (ACR) of 11.2%, supporting the authors'
conclusion that raltegravir does NOT accumulate intracellularly and
that PBMC concentrations can be predicted from plasma data by simple
multiplicative scaling.

# Population

Six healthy adult male volunteers were enrolled in Singapore (single
centre, National University Health System). Their median weight was
67.4 kg (range 56.5 to 87.1 kg) and median age 33.5 years (range 31 to
36). All subjects fasted overnight and received a single 400 mg oral
dose of raltegravir on day 1 at 8 a.m. Inclusion criteria required
smoking less than 10 cigarettes per day, hemoglobin > 10.9 g/dl,
platelet count >= 125,000/mm^3, creatinine clearance >= 60 mL/min, and
lipase or amylase < 1.1x ULN.

Paired blood samples for plasma and PBMCs were collected predose and at
4, 8, 12, 24, and 48 h post-dose; additional plasma-only samples were
collected at 0.5, 1, 1.5, 2, 3, 5, 6, and 10 h. PBMC concentrations
were computed assuming a cell volume of 0.4 pL.

Programmatic access:

```{r meta}
mod_meta <- readModelDb("Wang_2011_raltegravir")()$meta
mod_meta$population
```

# Source trace

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

| Equation / parameter | Value                     | Source location                                       |
|----------------------|---------------------------|-------------------------------------------------------|
| MAT                  | 0.887 h                   | Wang 2011 Table 1, Mean absorption time (RSE 27.6%)   |
| n transit            | 2 (structural; fixed)     | Wang 2011 Table 1, Number of transition compartments  |
| CL/F                 | 39.1 L/h                  | Wang 2011 Table 1, Clearance/bioavailability (RSE 20.8%) |
| V/F                  | 272 L                     | Wang 2011 Table 1, Volume distribution/bioavailability (RSE 21.8%) |
| ACR                  | 0.112                     | Wang 2011 Table 1, Accumulation ratio = 11.2% (RSE 35%) |
| IIV in MAT           | 66.2% CV                  | Wang 2011 Table 1, IIV in MAT (RSE 52.5%)             |
| IIV in CL/F          | 56.4% CV                  | Wang 2011 Table 1, IIV in CL/F (RSE 28.6%)            |
| IIV in V/F           | 56.7% CV                  | Wang 2011 Table 1, IIV in V/F (RSE 24.6%)             |
| IIV in ACR           | 104% CV                   | Wang 2011 Table 1, IIV in accumulation ratio (RSE 55.6%) |
| Correlation CL/F~V/F | 0.559                     | Wang 2011 Table 1, Correlation between CL/F and V/F (RSE 24.6%) |
| Exp. error plasma    | 86.7% (= sigma on log)    | Wang 2011 Table 1, Exponential error of plasma concns (RSE 15.4%) |
| Exp. error IC        | 64.9% (= sigma on log)    | Wang 2011 Table 1, Exponential error of intracellular concns (RSE 51.1%) |
| Transit-chain ktr    | 3 / MAT (derived)         | Wang 2011 Methods page 4091, MAT = (n+1)/ktr formula  |
| Plasma ODE chain     | depot -> T1 -> T2 -> central | Wang 2011 Fig. 1, Basic population PK model         |
| Intracellular model  | C_IC,obs = ACR * C_plasma * exp(eps) | Wang 2011 Methods page 4091, intracellular partition equation |

CV values are converted to internal log-scale variances using
`omega^2 = log(1 + CV^2)`. The CL/F~V/F covariance is
`0.559 * sqrt(omega^2_lcl * omega^2_lvc) = 0.155`.

# Molecular-weight unit conversion

Wang 2011 reports concentrations and AUC values in nmol/L (the units of
HIV-1 integrase IC50). The packaged model uses mass-balanced mg / mg/L
(matching `Lee_2016_raltegravir`). Conversions for cross-checking
against the paper's NCA table use raltegravir molecular weight
444.42 g/mol:

```{r mw}
mw_ral <- 444.42  # g/mol
nm_to_mgL <- function(nm)   nm * mw_ral * 1e-6           # nmol/L -> mg/L
nmh_to_mghL <- function(nmh) nmh * mw_ral * 1e-6         # nmol*h/L -> mg*h/L

# Published geometric-mean NCA (Wang 2011 Table; Pharmacokinetic parameters)
paper_nca <- tibble::tribble(
  ~param,                    ~unit_paper, ~value_paper, ~value_mgL,
  "Plasma Cmax",             "nmol/L",    2246,         nm_to_mgL(2246),
  "Plasma AUC0-12",          "nmol*h/L",  10776,        nmh_to_mghL(10776),
  "Plasma AUC0-inf",         "nmol*h/L",  13119,        nmh_to_mghL(13119),
  "Plasma C12",              "nmol/L",    338.7,        nm_to_mgL(338.7),
  "Plasma half-life (NCA)",  "h",         7.8,          NA,
  "Intracellular Cmax",      "nmol/L",    383,          nm_to_mgL(383),
  "Intracellular AUC0-12",   "nmol*h/L",  2073,         nmh_to_mghL(2073),
  "Intracellular AUC0-inf",  "nmol*h/L",  2435,         nmh_to_mghL(2435),
  "IC half-life (NCA)",      "h",         4.5,          NA
)
knitr::kable(paper_nca, digits = 4,
             caption = "Geometric-mean NCA values from Wang 2011 (Pharmacokinetic parameters paragraph), shown in the paper's nmol/L units and as the equivalent mg/L used by the packaged model.")
```

# Virtual cohort

Original observed data are not publicly available. The figures below
use a virtual population with the published Wang 2011 covariate
distribution (six healthy adult males, all in their early thirties,
roughly two-thirds of subjects between 60 and 75 kg, no other
covariates retained in the final model).

```{r cohort}
set.seed(202611)
n_sub <- 200L
events <- data.frame(
  id = seq_len(n_sub),
  time = 0,
  amt = 400,
  cmt = "depot",
  evid = 1L
)
times <- sort(unique(c(0,
                       seq(0.25, 12, by = 0.25),
                       seq(13, 24, by = 1),
                       seq(28, 48, by = 4),
                       c(0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10, 12, 24, 48))))
obs_plasma <- expand.grid(id = seq_len(n_sub), time = times)
obs_plasma$amt <- NA_real_
obs_plasma$cmt <- "Cc"
obs_plasma$evid <- 0L
obs_ic <- obs_plasma
obs_ic$cmt <- "Cic"
events <- dplyr::bind_rows(events, obs_plasma, obs_ic)
events <- events[order(events$id, events$time, events$evid > 0), ]
```

# Simulation

```{r simulate}
mod <- readModelDb("Wang_2011_raltegravir")
# rxSolve produces one row per (id, time, dvid). Since Cc and Cic are
# both computed as derived model variables on every row, each row carries
# both columns -- we deduplicate to one row per (id, time) for plotting
# and NCA.
sim <- rxode2::rxSolve(mod, events = events) |>
  as.data.frame() |>
  dplyr::distinct(id, time, .keep_all = TRUE)
nrow(sim)
```

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical, events = events) |>
  as.data.frame() |>
  dplyr::distinct(id, time, .keep_all = TRUE)
```

# Replicate Figure 2 (mean intracellular and plasma profiles)

```{r figure-2}
# Replicates Figure 2 of Wang 2011: mean plasma and intracellular
# raltegravir concentration vs time after a single 400 mg dose.
sim_long <- sim |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    plasma_median = median(Cc,  na.rm = TRUE),
    plasma_q05    = quantile(Cc,  0.05, na.rm = TRUE),
    plasma_q95    = quantile(Cc,  0.95, na.rm = TRUE),
    ic_median     = median(Cic, na.rm = TRUE),
    ic_q05        = quantile(Cic, 0.05, na.rm = TRUE),
    ic_q95        = quantile(Cic, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

# Reshape to long format for ggplot
sim_fig2 <- dplyr::bind_rows(
  dplyr::transmute(sim_long, time, compartment = "Plasma",
                   median = plasma_median, q05 = plasma_q05, q95 = plasma_q95),
  dplyr::transmute(sim_long, time, compartment = "Intracellular (PBMC)",
                   median = ic_median, q05 = ic_q05, q95 = ic_q95)
)

# Convert to nmol/L to match the paper's y-axis
sim_fig2 <- sim_fig2 |>
  dplyr::mutate(
    median_nM = median * 1e6 / mw_ral,
    q05_nM    = q05    * 1e6 / mw_ral,
    q95_nM    = q95    * 1e6 / mw_ral
  )

ic50 <- 7  # nmol/L (top end of HIV-1 integrase IC50 range, Wang 2011 Fig 2 caption)

ggplot(sim_fig2, aes(time, median_nM, colour = compartment, fill = compartment)) +
  geom_ribbon(aes(ymin = q05_nM, ymax = q95_nM), alpha = 0.15, colour = NA) +
  geom_line() +
  geom_hline(yintercept = ic50, linetype = "dashed", colour = "grey50") +
  annotate("text", x = 40, y = ic50, label = "IC50 = 7 nM", vjust = -0.5, colour = "grey50") +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Raltegravir concentration (nmol/L)",
       colour = NULL, fill = NULL,
       title = "Figure 2 -- Plasma and intracellular raltegravir after a single 400 mg dose",
       caption = "Replicates Figure 2 of Wang 2011. Lines are simulated medians; ribbons are 5-95 percentile bands.")
```

# Replicate Figure 4 (visual predictive check)

```{r figure-4}
# Replicates Figure 4 of Wang 2011: VPC of plasma and intracellular RAL
# with observed concentrations overlaid against simulated percentiles.
# Wang 2011 reported 97% of observations within the 95% prediction
# interval. We reproduce the plasma and intracellular panels separately.
ggplot(sim_long, aes(time, plasma_median)) +
  geom_ribbon(aes(ymin = plasma_q05, ymax = plasma_q95), alpha = 0.2) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Plasma raltegravir (mg/L)",
       title = "Figure 4A -- Plasma raltegravir VPC",
       caption = "Replicates Figure 4A of Wang 2011 (median and 5-95 percentile band).")

ggplot(sim_long, aes(time, ic_median)) +
  geom_ribbon(aes(ymin = ic_q05, ymax = ic_q95), alpha = 0.2) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Intracellular raltegravir (mg/L)",
       title = "Figure 4B -- Intracellular raltegravir VPC",
       caption = "Replicates Figure 4B of Wang 2011 (median and 5-95 percentile band).")
```

# PKNCA validation

We run PKNCA separately for plasma and intracellular concentrations,
each with the cohort label as a single-level treatment grouping (the
study had one treatment arm: 400 mg single dose).

```{r pknca-plasma}
# Plasma NCA
# rxSolve returns both Cc and Cic on every observation row; deduplicate to
# one row per (id, time) before passing to PKNCA.
sim_plasma <- sim |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::transmute(id, time, Cc, treatment = "400 mg single oral dose")

dose_df <- tibble::tibble(
  id = seq_len(n_sub),
  time = 0,
  amt = 400,
  treatment = "400 mg single oral dose"
)

conc_obj <- PKNCA::PKNCAconc(sim_plasma, 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,
  auclast = TRUE, aucinf.obs = TRUE, half.life = TRUE
)

nca_plasma <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
plasma_summary <- summary(nca_plasma)
knitr::kable(plasma_summary, caption = "Simulated plasma NCA (200-subject virtual cohort, 400 mg single dose).")
```

```{r pknca-intracellular}
# Intracellular NCA (treating Cic as the observation)
sim_ic <- sim |>
  dplyr::filter(!is.na(Cic), time > 0) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::transmute(id, time, Cc = Cic, treatment = "400 mg single oral dose (intracellular)")

dose_df_ic <- dose_df |>
  dplyr::mutate(treatment = "400 mg single oral dose (intracellular)")

conc_obj_ic <- PKNCA::PKNCAconc(sim_ic, Cc ~ time | treatment + id,
                                concu = "mg/L", timeu = "h")
dose_obj_ic <- PKNCA::PKNCAdose(dose_df_ic, amt ~ time | treatment + id, doseu = "mg")

nca_ic <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj_ic, dose_obj_ic, intervals = intervals))
ic_summary <- summary(nca_ic)
knitr::kable(ic_summary, caption = "Simulated intracellular NCA (200-subject virtual cohort, 400 mg single dose).")
```

## Comparison against published NCA

```{r pknca-compare}
# Pull the simulated geometric mean (or median) values back as numeric
# from the typical-value (no-IIV, no-residual) simulation, and compare
# to the published geometric means converted from nmol/L to mg/L.
get_typ <- function(df, col) df[[col]][which.max(df[[col]])]
typ <- sim_typical |> dplyr::filter(!is.na(Cc))
typ_plasma_cmax <- max(typ$Cc,  na.rm = TRUE)
typ_ic_cmax     <- max(typ$Cic, na.rm = TRUE)
typ_plasma_c12  <- typ$Cc[abs(typ$time - 12) < 1e-6][1]

cmp <- tibble::tribble(
  ~metric,             ~paper_mgL,                ~typical_value_mgL,
  "Plasma Cmax",       nm_to_mgL(2246),            typ_plasma_cmax,
  "Plasma C12",        nm_to_mgL(338.7),           typ_plasma_c12,
  "Intracellular Cmax", nm_to_mgL(383),            typ_ic_cmax
)
cmp$ratio_typ_over_paper <- cmp$typical_value_mgL / cmp$paper_mgL
knitr::kable(cmp, digits = 4,
             caption = "Wang 2011 geometric-mean NCA (converted to mg/L) vs. typical-value simulation from the packaged model. Differences arise because the typical value (no IIV) does not reproduce a small-N geometric mean exactly, and because the one-compartment fit gives an apparent half-life t1/2 = ln(2) * V/CL = 4.82 h compared to the 7.8 h NCA-derived half-life reported in the paper (see Assumptions and deviations).")
```

# Assumptions and deviations

- **Number of transit compartments fixed at 2 (structural).** Wang 2011
  reports n = 2 as the optimal number of transition compartments selected
  during model development, not as an estimated theta with a precision
  estimate. The model encodes n = 2 implicitly in the ODE chain
  (`depot -> transit1 -> transit2 -> central`) and ties the transit-chain
  rate constant to MAT through `ktr = (n + 1) / MAT = 3 / MAT`. There is
  no `fixed()` ini entry because n is not a theta; it is a structural
  choice that determines the number of `d/dt` lines.
- **Mass-balance unit choice.** The paper reports concentrations in
  nmol/L; the packaged model carries concentrations in mg/L
  (`units$concentration = "mg/L"`) and the dose is supplied as `amt = 400`
  with `units$dosing = "mg"`. The cross-checking conversion to nmol/L uses
  raltegravir MW = 444.42 g/mol. This matches `Lee_2016_raltegravir`'s
  unit convention.
- **Exponential residual error interpretation.** The paper Methods state
  "Interindividual and residual unexplained variability were estimated
  with an exponential error model" (Wang 2011 page 4091). In NONMEM this
  is Y = F * exp(eps) with eps ~ N(0, sigma^2). nlmixr2's `~ lnorm(expSd)`
  is the exact analogue with `expSd = sigma` on the log scale. Wang 2011
  Table 1 reports residual error as a percentage in the same column it
  reports IIV %CV; following the standard NONMEM `$SIGMA`-on-log-scale
  reporting convention, the percentage is read as `sigma * 100%`
  (`expSd = 0.867` for plasma; `expSd_Cic = 0.649` for intracellular). For
  IIV the percentages are read as linear-space %CV and converted via
  `omega^2 = log(1 + CV^2)`, which is the standard nlmixr2lib convention
  matching `Lee_2016_raltegravir`.
- **Apparent vs. NCA half-life mismatch.** The one-compartment fit gives
  an apparent terminal half-life `t1/2 = ln(2) * V / CL = 0.693 * 272 / 39.1
  ~= 4.82 h`, whereas the paper's NCA reports `t1/2 = 7.8 h` (95% CI
  5.5-11.3 h) using terminal-phase regression over 24 and 48 h time points.
  The paper Discussion acknowledges this: "our attempts to fit the data
  to a two-compartment model were not successful, most likely due to
  little data collected during the beta-elimination phase". The simulated
  trajectories therefore decline somewhat faster than the late-time
  observed data; this is a feature of the published 1-compartment model,
  not a transcription error. Larger studies cited by the authors
  (Dickinson et al., Iwamoto et al.) reported either two-compartment fits
  or different one-compartment parameter estimates -- consult those papers
  for a model intended for late-time extrapolation.
- **Accumulation ratio is a misnomer for raltegravir.** Wang 2011 names
  the parameter ACR (accumulation ratio) but the point estimate of 11.2%
  means raltegravir does NOT accumulate intracellularly. The packaged
  model maps this parameter onto the canonical paper-named bare
  parameter `frac` (a unitless fraction in (0, 1)); the label and inline
  comments preserve the paper's "accumulation ratio" name and emphasise
  the (0, 1) interpretation. The paper itself states "this result is
  consistent with simple diffusion of unbound RAL into cells" (Discussion,
  page 4094).
- **Intracellular model has no separate compartment kinetics.** Wang 2011
  Eq.: `C_IC,obs = ACR * C_plasma,pred * exp(eps_IC)` -- the
  intracellular concentration is an empirical multiplicative transform of
  the predicted plasma concentration, not a separate compartment with
  its own ODE. The packaged model implements this exactly as
  `Cic <- frac * Cc` followed by `Cic ~ lnorm(expSd_Cic)`. There is no
  intracellular ODE state, and PBMC isolation / processing effects
  (which the paper Discussion notes contribute to the high intracellular
  residual error of 64.9%) enter only through `expSd_Cic`.
- **Sample size is small.** Six subjects produce wide confidence
  intervals on every parameter (RSE 15-56%). The visual predictive
  check in Wang 2011 Fig. 4 covered 97% of observed concentrations, but
  the authors explicitly cautioned (Discussion, page 4094) that "large
  pharmacokinetic studies are required to further characterize the
  pharmacokinetics of RAL before the population pharmacokinetic model
  could be applied in a clinical setting." Treat the simulated NCA
  values as illustrative rather than predictive for new dosing
  scenarios.
