---
title: "Voriconazole (Chen 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Voriconazole (Chen 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Chen_2015_voriconazole"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Chen_2015_voriconazole"))$description`
* Article: <https://doi.org/10.1248/bpb.b14-00768>

# Population

The model was developed from a prospective observational study at the First
Affiliated Hospital of Guangzhou Medical University (Guangzhou, China) between
March 2012 and May 2013 (Chen 2015 Materials and Methods, "Study Design").
Adult intensive-care-unit patients with pulmonary disease who received
intravenous voriconazole (VFEND, Pfizer Ireland) as therapeutic drug monitoring
were enrolled.

The final analysis included 62 adult Chinese ICU patients (42 male, 20 female)
aged 19-90 years (mean 59.7 +/- 16.7) with body weight 41-84 kg (mean 60.1 +/-
10.0); 240 plasma voriconazole concentration measurements were collected (Chen
2015 Table 1). Severity at enrolment: APACHE II 21.6 +/- 13.8, SOFA median 4
(range 0-20); observed in-hospital mortality 14.5%. The cohort had clinical
indication for invasive fungal infection, with culture-positive pathogens in 52
patients (aspergillosis n = 37, candidiasis n = 12, other n = 3). Pulmonary
infection was the dominant infected location (n = 45). Hepatotoxicity
attributable to voriconazole occurred in 24.2% of patients during therapy.

Dosing: 300 mg intravenous loading dose, then 200 mg intravenous infusion every
12 hours. Therapeutic drug monitoring sampling commenced at trough before the
sixth maintenance dose, treated as steady-state (72 hours after the loading
dose). Plasma samples were drawn at 0.5, 1, 1.5, 2, 4, 6, 9, and 12 h after the
start of an infusion, or at any two of those points per occasion.
Bioanalytical: HPLC-UV with LLOQ 70 ng/mL and calibration range 208-20800
ng/mL.

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

# Source trace

Every parameter in the model file carries an inline source-location comment.
The table below collects the entries in one place.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (CL at DBIL reference 2.6 umol/L) | 4.28 L/h | Table 2 (final model) |
| `lvc` (Vd)                              | 93.4 L   | Table 2 (final model) |
| `e_dbil_cl` (DBIL power exponent on CL) | -0.40    | Table 2 (final model); typical-value equation `CL = theta_CL * (DBIL/2.6)^theta_DBIL` |
| IIV CL (CV%)                            | 72.94%   | Table 2 row omega_CL |
| IIV Vd (CV%)                            | 26.50%   | Table 2 row omega_Vd |
| Proportional residual error             | 13.0%    | Table 2 row delta |
| 1-cmt IV first-order elimination (NONMEM ADVAN1 TRANS2) | n/a | Results, "Population Pharmacokinetic Analysis" paragraph |
| Exponential IIV; constant-CV (proportional) residual model | n/a | Results, "Population Pharmacokinetic Analysis" paragraph |
| Reported steady-state median Cmin       | 3.26 ug/mL at t = 0; 3.76 ug/mL at t = 12 h | Results, "Analytical Quantifications of VRC" |
| Reported steady-state AUC0-12 (mean +/- SD) | 47.08 +/- 26.28 ug.h/mL | Results, "Plasma Concentration-Effect Relationships" |

# Virtual cohort

Original observed concentrations are not openly available. The virtual cohort
below mirrors the demographics in Chen 2015 Table 1, with direct bilirubin
sampled from a log-normal distribution centred near the cohort mean of
3.16 umol/L (range observed 0.2-10.3 umol/L).

```{r cohort}
set.seed(20150427)

n_subjects <- 200L

# Body weight: cohort mean 60.13 kg +/- SD 10.03, range 41-84.
wt <- pmin(pmax(rnorm(n_subjects, mean = 60.13, sd = 10.03), 41), 84)

# Age (kept as a descriptive column; not used by the model).
age <- pmin(pmax(rnorm(n_subjects, mean = 59.71, sd = 16.67), 19), 90)

# Direct bilirubin distribution. Chen 2015 Table 1 reports DBIL mean 3.16 +/-
# SD 2.21 umol/L with range 0.2-10.3. The distribution is right-skewed in
# clinical data; a log-normal sampler with the cohort mean as the median is a
# reasonable approximation that respects the strictly positive support.
log_sd  <- sqrt(log(1 + (2.21 / 3.16)^2))
log_mu  <- log(3.16) - 0.5 * log_sd^2
dbil    <- pmin(pmax(rlnorm(n_subjects, log_mu, log_sd), 0.2), 10.3)

demo <- tibble(
  id   = seq_len(n_subjects),
  WT   = wt,
  AGE  = age,
  DBIL = dbil
)
stopifnot(!anyDuplicated(demo$id))
```

# Simulation

The clinical regimen is 300 mg loading then 200 mg q12h IV maintenance. Chen
2015 used NONMEM ADVAN1 TRANS2 (1-compartment first-order elimination) without
an explicit infusion duration, so the model treats each dose as a bolus. The
simulation runs five days (10 maintenance doses after the loading dose), with
dense sampling across the day-5 dosing interval that the paper reports for
its steady-state TDM observations.

```{r build-events}
loading_dose <- 300
maint_dose   <- 200
n_maint      <- 10L                        # 5 days x 2 doses/day
ii           <- 12

build_events <- function(demo, sim_hours = 120) {
  loading <- demo |>
    mutate(amt  = loading_dose,
           evid = 1L,
           cmt  = "central",
           ii   = NA_real_,
           addl = NA_integer_,
           time = 0) |>
    select(id, time, amt, evid, cmt, ii, addl, WT, AGE, DBIL)

  maintenance <- demo |>
    mutate(amt  = maint_dose,
           evid = 1L,
           cmt  = "central",
           ii   = ii,
           addl = n_maint - 1L,
           time = 12) |>
    select(id, time, amt, evid, cmt, ii, addl, WT, AGE, DBIL)

  obs_times <- sort(unique(c(seq(0,  24, by = 1),
                             seq(96, sim_hours, by = 0.5))))
  obs <- demo |>
    select(id, WT, AGE, DBIL) |>
    tidyr::crossing(time = obs_times) |>
    mutate(amt  = NA_real_,
           evid = 0L,
           cmt  = NA_character_,
           ii   = NA_real_,
           addl = NA_integer_)

  bind_rows(loading, maintenance, obs) |>
    arrange(id, time, desc(evid))
}

events <- build_events(demo)
```

```{r simulate}
mod <- rxode2::rxode2(readModelDb("Chen_2015_voriconazole"))

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("WT", "AGE", "DBIL")
) |> as.data.frame()

mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(
  mod_typical, events = events,
  keep = c("WT", "AGE", "DBIL")
) |> as.data.frame()
```

# Replicate published figures

## Figure 1 -- steady-state mean concentration-time profile

Chen 2015 Figure 1B shows the mean plasma concentration-time profile over the
12 h interval at steady state. The simulated typical-value profile (no random
effects) tracks an approximate Cmin near 3 ug/mL and Cmax near 5 ug/mL,
consistent with the reported observed median Cmin of 3.26-3.76 ug/mL.

```{r figure-1b, fig.cap = "Replicates Figure 1B of Chen 2015: mean steady-state plasma concentration vs. time over the 12 h dosing interval (data here shown for hours 96-108). The solid line is the typical-value (no random effects) prediction; the ribbon shows the simulated 10th-90th percentile of the full virtual cohort."}
ss_profile <- sim |>
  filter(time >= 96, time <= 108) |>
  mutate(time_after_dose = time - 96) |>
  group_by(time_after_dose) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ss_typ <- sim_typical |>
  filter(time >= 96, time <= 108) |>
  mutate(time_after_dose = time - 96) |>
  group_by(time_after_dose) |>
  summarise(Cc_typ = median(Cc, na.rm = TRUE), .groups = "drop")

ggplot(ss_profile, aes(time_after_dose, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20) +
  geom_line() +
  geom_line(data = ss_typ, aes(time_after_dose, Cc_typ),
            linetype = "dashed", colour = "darkred") +
  labs(x = "Time after dose (h)", y = "Cc (ug/mL)",
       title = "Steady-state day-5 concentration profile",
       caption = "Replicates Figure 1B of Chen 2015 (mean steady-state profile).")
```

# PKNCA validation

A standard NCA over the day-5 dosing interval (steady state, 12 h interval)
gives Cmax, Cmin / Ctrough, and AUClast. The simulated NCA can be compared
against the steady-state AUC0-12 and Cmin reported in the paper.

```{r pknca}
last_dose_time  <- 96  # 9th dose (5th day morning); window 96-108 h
nca_window <- sim |>
  filter(time >= last_dose_time, time <= last_dose_time + 12) |>
  mutate(time_after_dose = time - last_dose_time,
         treatment       = "200 mg q12h IV") |>
  select(id, time = time_after_dose, Cc, treatment)

dose_df <- demo |>
  mutate(time      = 0,
         amt       = maint_dose,
         treatment = "200 mg q12h IV") |>
  select(id, time, amt, treatment)

conc_obj  <- PKNCA::PKNCAconc(nca_window, Cc ~ time | treatment + id)
dose_obj  <- PKNCA::PKNCAdose(dose_df,    amt ~ time | treatment + id)
intervals <- data.frame(start = 0, end = 12,
                        cmax = TRUE, tmax = TRUE,
                        auclast = TRUE,
                        cmin = TRUE, ctrough = TRUE)
nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- suppressMessages(suppressWarnings(PKNCA::pk.nca(nca_data)))
nca_summary <- summary(nca_res)
knitr::kable(nca_summary,
             caption = "Day-5 NCA on the simulated cohort (steady-state 12 h interval, 200 mg q12h IV after a 300 mg loading dose).")
```

## Comparison against published NCA

Chen 2015 reports the steady-state AUC0-12 as 47.08 +/- 26.28 ug.h/mL across
17 patients with complete sampling profiles (Results, "Plasma
Concentration-Effect Relationships"), and median observed Cmin of 3.26 ug/mL
at the start of the dosing interval and 3.76 ug/mL at the end (Results,
"Analytical Quantifications of VRC"). The expected typical-value AUC0-12 at
steady state from this model is `Dose / CL = 200 / 4.28 = 46.7 ug.h/mL`,
matching the published 47.08 ug.h/mL within < 1%.

```{r typical-comparison}
trough_sim <- sim |>
  filter(time == 108) |>
  summarise(median_trough = median(Cc),
            Q10           = quantile(Cc, 0.10),
            Q90           = quantile(Cc, 0.90))

auc_sim <- nca_res$result |>
  filter(PPTESTCD == "auclast") |>
  summarise(median_auc = median(PPORRES, na.rm = TRUE),
            Q10        = quantile(PPORRES, 0.10, na.rm = TRUE),
            Q90        = quantile(PPORRES, 0.90, na.rm = TRUE))

tbl <- tibble::tibble(
  metric = c("Chen 2015 steady-state AUC0-12 (mean +/- SD, ug.h/mL)",
             "Chen 2015 steady-state median Cmin (ug/mL)",
             "Simulated steady-state AUC0-12 (median, 10-90 percentile, ug.h/mL)",
             "Simulated steady-state Cmin at t = 108 h (median, 10-90 percentile, ug/mL)"),
  value  = c("47.08 +/- 26.28",
             "3.26 (start of interval) and 3.76 (end of interval)",
             sprintf("%.1f (%.1f-%.1f)",
                     auc_sim$median_auc, auc_sim$Q10, auc_sim$Q90),
             sprintf("%.2f (%.2f-%.2f)",
                     trough_sim$median_trough,
                     trough_sim$Q10, trough_sim$Q90))
)
knitr::kable(tbl, caption = "Simulated steady-state day-5 PK vs. Chen 2015 reported values.")
```

The simulated AUC0-12 matches the published value to within 1% (it is the
direct algebraic identity `Dose / CL` at steady state). The simulated Cmin
sits in the same range as the reported observed median trough of 3.26-3.76
ug/mL.

# Assumptions and deviations

- **IV infusion is treated as bolus (no infusion duration).** Voriconazole IV
  is clinically administered as a 1-2 hour infusion. Chen 2015 used NONMEM
  ADVAN1 TRANS2 (1-compartment IV with first-order elimination) without an
  explicit infusion duration, so the model treats each dose as a bolus.
  Users who need to represent the infusion explicitly can add
  `dur(central) <- ...` and pass `rate = -2` in event records.
- **No covariates on Vd.** Chen 2015's stepwise covariate analysis identified
  direct bilirubin as the only retained covariate on CL; no covariate was
  retained for Vd. The model applies the covariate effect to CL only.
- **DBIL is treated as a baseline / time-invariant covariate.** Chen 2015
  reports DBIL as a baseline laboratory value (Table 1) and does not describe
  it as time-varying within the dosing interval. Simulated patients carry a
  single DBIL value across the 5-day simulation.
- **Bioavailability not parameterised.** Voriconazole has high oral
  bioavailability in adults (~96%) but the Chen 2015 dataset is IV-only and
  Ka and F were excluded from estimation. The model file does not
  parameterise `lka` or `lfdepot`; users modelling oral voriconazole must
  use a different model.
- **No CYP2C19 covariate.** The Discussion notes that CYP2C19 genotype is a
  known driver of voriconazole PK variability, but the dataset did not
  include CYP2C19 genotyping. The paper explicitly lists this as a study
  limitation. The model has no CYP2C19 covariate.
- **Cohort defined at 200 subjects.** The virtual cohort uses 200 subjects
  (versus the source paper's 62) so the per-time-point percentiles in the
  Figure 1B reproduction and the Cmin / AUC summary tables are smooth. The
  vignette runs within the 5-minute pkgdown gate at this cohort size.
- **DBIL distribution assumed log-normal.** The source paper reports DBIL
  mean 3.16 +/- 2.21 umol/L (range 0.2-10.3) but does not publish a
  histogram. The virtual cohort samples DBIL from a log-normal distribution
  matching the reported mean and SD, truncated to the reported range.
- **Shrinkage on Vd.** Chen 2015 reports shrinkage on Vd of 48.32% (Results,
  "Population Pharmacokinetic Analysis"), indicating that the dataset
  provides limited information about between-subject Vd variability. The
  model retains the reported 26.5% IIV-Vd; users should be aware that
  individual Vd predictions in TDM contexts will be heavily shrunk toward
  the typical value.
- **Sampling design.** The source dataset is opportunistic TDM sampling
  with samples drawn at 0.5, 1, 1.5, 2, 4, 6, 9, and 12 h after the start
  of an infusion (or any two of those points per occasion). The 240
  observations were divided across 62 patients, so most patients
  contributed sparse profiles.
