---
title: "Voriconazole (Han 2010)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Voriconazole (Han 2010)}
  %\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("Han_2010_voriconazole"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Han_2010_voriconazole"))$description`
* Article: <https://doi.org/10.1128/AAC.00504-10>

# Population

The model was developed from a prospective single-center observational study
at the University of Pittsburgh Medical Center of adult lung transplant
recipients who started voriconazole prophylactically immediately after
transplantation (Han 2010 Materials and Methods, "Patients").

Thirteen lung transplant recipients (7 male, 6 female; 12 of 13 Caucasian)
aged 19-70 years (mean 50.9 +/- 16.1) with body weight 46-91 kg
(mean 68.0 +/- 15.2) were enrolled. Primary diagnoses leading to transplant
were cystic fibrosis (n = 3), emphysema (n = 5), idiopathic pulmonary fibrosis
(n = 4), and scleroderma (n = 1). The day of the oral pharmacokinetic study
was on average 8.5 days post-transplant (range 3-19 days). All patients
received tacrolimus as primary immunosuppression. One patient did not
complete the oral study (Han 2010 Table 1).

Dosing: two 2-hour intravenous infusions of 6 mg/kg every 12 hours
immediately post-transplant, followed by oral 200 mg every 12 hours for
3 months. Blood samples were collected at pre-dose and at 0.5, 1, 1.5, 2, 4,
6, 8, and 12 h following the second intravenous dose and following an oral
dose (the 5th to 37th dose, mean 15th). Bioanalysis: HPLC with assay precision
1.3-9.0%, bias 0.7-3.1%, and linearity range 0.2-9 ug/mL (R^2 = 0.9998).
NONMEM 6.2.0 (GloboMax) with FOCE-I.

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

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| Two-compartment with first-order absorption + elimination | n/a | Results, "Population pharmacokinetic analysis" paragraph 1 |
| Exponential IIV; combined additive + proportional residual error | n/a | Methods, "Population pharmacokinetic analysis" |
| `lka` (absorption rate, 1/h)                    | 0.591  | Results, "Population pharmacokinetic analysis" |
| `lcl` (clearance, L/h)                          | 3.45   | Results, "Population pharmacokinetic analysis" |
| `lvc` (central volume, L)                       | 54.7   | Results, "Population pharmacokinetic analysis" |
| `lvp` (peripheral volume, L)                    | 143    | Results, "Population pharmacokinetic analysis" |
| `lq`  (inter-compartmental clearance, L/h)      | 22.6   | Results, "Population pharmacokinetic analysis" |
| `lfdepot` (oral bioavailability, fraction)      | 0.459  | Results, "Population pharmacokinetic analysis" |
| IIV ka  (CV%)                                   | 115.2% | Results, "Population pharmacokinetic analysis" |
| IIV CL  (CV%)                                   | 107%   | Results, "Population pharmacokinetic analysis" |
| IIV Vc  (CV%)                                   | 78.4%  | Results, "Population pharmacokinetic analysis" |
| IIV Vp  (CV%)                                   | 88.3%  | Results, "Population pharmacokinetic analysis" |
| IIV Q   (CV%)                                   | 50.1%  | Results, "Population pharmacokinetic analysis" |
| IIV F   (CV%)                                   | 82.9%  | Results, "Population pharmacokinetic analysis" |
| Proportional residual SD (fraction)             | 0.31   | Results, "Population pharmacokinetic analysis" |
| Additive residual SD (ug/mL)                    | 0.49   | Results, "Population pharmacokinetic analysis" |
| Reported Tmax (oral, mean +/- SD)               | 1.9 +/- 1.3 h | Results, "Noncompartmental analysis" |
| Reported Cmax (oral, mean +/- SD)               | 3.6 +/- 2.6 ug/mL | Results, "Noncompartmental analysis" |
| Reported Cmax (IV infusion, mean +/- SD)        | 5.9 +/- 2.2 ug/mL | Results, "Noncompartmental analysis" |
| Covariate sub-model 1 (CF on F)                 | `F = 0.107 + 0.72 * NCF` | Results, "Model 1: cystic fibrosis (CF)" |
| Covariate sub-model 2 (POT on F)                | `F = 0.619 * POT / (POT + 1.97)` | Results, "Model 2: postoperative time (POT)" |
| Covariate sub-model 3 (WT on Vp)                | `Vp = 148 * (WT/68)^3.56` | Results, "Model 3: body weight (WT)" |

# Virtual cohort

The original observed concentrations are not openly available. The virtual
cohort below mirrors the demographics in Han 2010 Table 1, with body weight
sampled from a truncated normal distribution matching the reported mean and
SD.

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

n_subjects <- 200L

wt <- pmin(pmax(rnorm(n_subjects, mean = 68.0, sd = 15.2), 46), 91)

age <- pmin(pmax(rnorm(n_subjects, mean = 50.9, sd = 16.1), 19), 70)

# Postoperative time on the day of the oral study (Han 2010 Table 1
# reports a mean 8.5 days +/- SD 4.4 days, range 3-19). Used for the
# Model 2 (POT on F) sub-model below; not used by the base model.
pot_days <- pmin(pmax(rnorm(n_subjects, mean = 8.5, sd = 4.4), 3), 19)

# Cystic-fibrosis indicator. 3 of 13 patients (~23%) had CF.
cf <- rbinom(n_subjects, size = 1, prob = 3 / 13)

demo <- tibble(
  id   = seq_len(n_subjects),
  WT   = wt,
  AGE  = age,
  POT_days = pot_days,
  CF   = cf
)
stopifnot(!anyDuplicated(demo$id))
```

# Simulation

The clinical regimen per Han 2010 is two 2-hour intravenous infusions of
6 mg/kg every 12 h, immediately followed by oral 200 mg every 12 h. The
simulation below builds an event table covering the first dosing interval of
each route (IV infusion 0-12 h and oral 0-12 h, simulated separately) so the
day-1 IV and steady-state-like oral profiles can be compared against Han
2010 Figure 1 panels (a)-(d).

```{r build-events}
# Two separate simulation cohorts: one IV-infusion-only and one oral-only,
# each with disjoint id ranges so they can be bind_rows-ed without
# rxSolve silently merging IDs.

make_iv_cohort <- function(demo, t_obs, id_offset = 0L) {
  d <- demo |>
    mutate(id = id + id_offset,
           iv_dose = WT * 6)  # 6 mg/kg
  dose <- d |>
    transmute(id, time = 0, amt = iv_dose, evid = 1L,
              cmt = "central", rate = iv_dose / 2,  # 2 h infusion
              WT, AGE, POT_days, CF,
              route = "IV 6 mg/kg over 2 h")
  obs <- d |>
    select(id, WT, AGE, POT_days, CF) |>
    tidyr::crossing(time = t_obs) |>
    mutate(amt = NA_real_, evid = 0L, cmt = NA_character_,
           rate = NA_real_, route = "IV 6 mg/kg over 2 h")
  bind_rows(dose, obs) |> arrange(id, time, desc(evid))
}

make_oral_cohort <- function(demo, t_obs, id_offset = 0L) {
  d <- demo |> mutate(id = id + id_offset)
  dose <- d |>
    transmute(id, time = 0, amt = 200, evid = 1L,
              cmt = "depot", rate = NA_real_,
              WT, AGE, POT_days, CF,
              route = "Oral 200 mg")
  obs <- d |>
    select(id, WT, AGE, POT_days, CF) |>
    tidyr::crossing(time = t_obs) |>
    mutate(amt = NA_real_, evid = 0L, cmt = NA_character_,
           rate = NA_real_, route = "Oral 200 mg")
  bind_rows(dose, obs) |> arrange(id, time, desc(evid))
}

t_obs <- sort(unique(c(seq(0, 12, by = 0.25))))
events <- bind_rows(
  make_iv_cohort(demo,   t_obs, id_offset = 0L),
  make_oral_cohort(demo, t_obs, id_offset = n_subjects)
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

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

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

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

# Replicate published figures

## Figure 1 -- concentration-time profiles by route

Han 2010 Figure 1 panels (a)-(d) show observed plasma voriconazole
concentration-time profiles over the first 12 h of an IV infusion and an
oral dose. Median observed Cmax for the IV infusion was 5.9 +/- 2.2 ug/mL
and for the oral dose was 3.6 +/- 2.6 ug/mL with Tmax around 1.9 +/- 1.3 h
(Han 2010 Results, "Noncompartmental analysis"). The simulated typical-value
profile and the 10-90 percentile band of the virtual cohort below reproduce
those qualitative features.

```{r figure-1, fig.cap = "Replicates Figure 1 of Han 2010: simulated plasma voriconazole concentration-time profiles over the first 12 h of an intravenous infusion (6 mg/kg over 2 h) and an oral dose (200 mg). The dashed line is the typical-value (no random effects) profile; the ribbon shows the 10-90 percentile of a 200-subject virtual cohort."}
plot_data <- sim |>
  filter(time > 0) |>
  group_by(route, time) |>
  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"
  )

typical_data <- sim_typical |>
  filter(time > 0) |>
  group_by(route, time) |>
  summarise(Cc_typ = median(Cc, na.rm = TRUE), .groups = "drop")

ggplot(plot_data, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20) +
  geom_line() +
  geom_line(data = typical_data, aes(time, Cc_typ),
            linetype = "dashed", colour = "darkred") +
  facet_wrap(~route) +
  labs(x = "Time after dose (h)", y = "Cc (ug/mL)",
       title = "Voriconazole concentration-time profiles by route",
       caption = "Replicates Figure 1 of Han 2010.")
```

# PKNCA validation

A standard NCA over the 12 h interval after each dose gives Cmax, Tmax,
and AUClast. The simulated NCA can be compared against the noncompartmental
values reported in Han 2010 Results, "Noncompartmental analysis."

```{r pknca}
nca_window <- sim |>
  filter(time >= 0, time <= 12) |>
  mutate(treatment = route) |>
  select(id, time, Cc, treatment)

dose_df <- events |>
  filter(evid == 1) |>
  mutate(treatment = route) |>
  select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(nca_window, 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      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE,
  clast.obs = 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 = "Simulated NCA over the first 12 h by route (200-subject virtual cohort).")
```

## Comparison against published NCA

Han 2010 reports the following noncompartmental values (Results,
"Noncompartmental analysis"):

* Oral Tmax: 1.9 +/- 1.3 h
* Oral Cmax: 3.6 +/- 2.6 ug/mL
* IV-infusion Cmax: 5.9 +/- 2.2 ug/mL

```{r typical-comparison}
nca_tbl <- as.data.frame(nca_res$result) |>
  group_by(treatment, PPTESTCD) |>
  summarise(median_value = median(PPORRES, na.rm = TRUE),
            Q10          = quantile(PPORRES, 0.10, na.rm = TRUE),
            Q90          = quantile(PPORRES, 0.90, na.rm = TRUE),
            .groups      = "drop") |>
  filter(PPTESTCD %in% c("cmax", "tmax"))

published <- tibble::tribble(
  ~treatment,             ~PPTESTCD, ~published,
  "IV 6 mg/kg over 2 h",  "cmax",    "5.9 +/- 2.2 ug/mL",
  "Oral 200 mg",          "cmax",    "3.6 +/- 2.6 ug/mL",
  "Oral 200 mg",          "tmax",    "1.9 +/- 1.3 h"
)

comparison <- published |>
  left_join(nca_tbl, by = c("treatment", "PPTESTCD")) |>
  mutate(simulated = sprintf("%.2f (%.2f-%.2f)", median_value, Q10, Q90)) |>
  select(treatment, PPTESTCD, published, simulated)

knitr::kable(comparison,
             caption = "Simulated NCA vs Han 2010 reported noncompartmental values.")
```

The simulated medians track the published means reasonably well given the
small (n = 13) source cohort and the large reported IIV (CV 107% on CL,
115% on ka). The IIV-driven 10-90 percentile spans straddle the published
mean +/- SD in all three comparisons.

# Covariate sub-models (Han 2010 Models 1, 2, 3)

The base model encoded in `Han_2010_voriconazole` carries no covariates.
Han 2010 reports three separate single-covariate sub-models. Each was
significant (P < 0.01) on its own; the paper notes that a fully combined
final model was also built but did not improve the goodness-of-fit
materially, and the parameter estimates of the combined model are not
reported (Han 2010 Discussion paragraph beginning "A final model was also
built using a standard forward addition and reverse removal approach").
The three sub-models below reproduce the reported equations using the
published parameter estimates. Because the model file does not encode the
covariate equations, the sub-model results are computed in R from the
published equations and applied as a post-simulation transform of the base
model's predictions where the covariate enters multiplicatively (CF and POT
on F; WT on Vp).

## Model 1 -- cystic fibrosis (CF) on bioavailability

Han 2010 Results, "Model 1: cystic fibrosis (CF)":
`F = F_CF + F' * NCF`, with `F_CF = 0.107` and `F' = 0.72`, and NCF = 1 for
non-CF patients, 0 for CF. Equivalent: F = 0.827 for non-CF, F = 0.107 for
CF. Compared to the base-model `F = 0.459`, the typical CF patient's oral
AUC is reduced approximately 4.3-fold and the typical non-CF patient's is
increased approximately 1.8-fold.

```{r model-1-cf}
# The base model uses F = 0.459. Model 1 partitions this into F = 0.107 for
# CF and F = 0.827 for non-CF; the per-subject Model-1 F can be applied as
# a multiplicative scaling of the base-model oral concentration:
#   Cc_model1 = Cc_base * (F_model1 / F_base).
F_base <- 0.459
F_cf      <- 0.107
F_noncf   <- 0.107 + 0.72  # = 0.827

model1_scale <- function(cf) ifelse(cf == 1, F_cf / F_base, F_noncf / F_base)

oral_model1 <- sim |>
  filter(route == "Oral 200 mg", time > 0) |>
  mutate(Cc_model1 = Cc * model1_scale(CF),
         cf_label  = ifelse(CF == 1, "CF", "non-CF"))

model1_plot <- oral_model1 |>
  group_by(cf_label, time) |>
  summarise(
    Q10 = quantile(Cc_model1, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc_model1, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc_model1, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(model1_plot, aes(time, Q50, colour = cf_label, fill = cf_label)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20, colour = NA) +
  geom_line() +
  labs(x = "Time after oral dose (h)", y = "Cc (ug/mL)",
       title = "Han 2010 Model 1: CF vs non-CF bioavailability on Cc",
       caption = "Approximates Figure 5(a) of Han 2010 (CF vs non-CF voriconazole concentration profiles).")
```

## Model 2 -- postoperative time (POT) on bioavailability

Han 2010 Results, "Model 2: postoperative time (POT)":
`F(POT) = F_max * POT / (POT + Fc)`, with `F_max = 0.619` and `Fc = 1.97 h`.
This Hill-form increase in F with POT reaches half its maximum at POT = 1.97 h
and approaches `F_max = 0.619` asymptotically. The paper notes that
bioavailability reaches steady levels within about one week posttransplant
(Han 2010 Figure 4).

```{r model-2-pot}
# Sweep POT from 0 to 14 days (=336 h) at the typical-value F equation.
F_max <- 0.619
Fc    <- 1.97  # hours, per Han 2010 Results

pot_grid <- tibble(POT_h = seq(0, 14 * 24, length.out = 200)) |>
  mutate(POT_days = POT_h / 24,
         F_model2 = F_max * POT_h / (POT_h + Fc))

ggplot(pot_grid, aes(POT_days, F_model2)) +
  geom_line() +
  geom_hline(yintercept = F_max, linetype = "dotted") +
  geom_hline(yintercept = F_base, linetype = "dashed", colour = "darkred") +
  labs(x = "Postoperative time (days)", y = "Bioavailability F",
       title = "Han 2010 Model 2: F(POT) Hill-form rise to F_max",
       caption = "Reproduces Han 2010 Figure 4 typical-value curve. Dashed red line is the base-model F = 0.459; dotted line is the Model-2 F_max = 0.619.")
```

## Model 3 -- body weight (WT) on peripheral volume

Han 2010 Results, "Model 3: body weight (WT)":
`Vp(WT) = TV(Vp) * (WT/68)^a` with `TV(Vp) = 148 L` and `a = 3.56`. The
reference weight 68 kg is the cohort mean (Han 2010 Table 1). Note that
TV(Vp) = 148 L in Model 3 differs slightly from the base-model Vp = 143 L
because each sub-model re-estimates the parameters it modifies; the
difference is within the IIV (88.3% CV on Vp).

```{r model-3-wt}
TV_vp_m3 <- 148
a_wt_vp  <- 3.56
ref_wt   <- 68

wt_grid <- tibble(WT = seq(46, 91, by = 0.5)) |>
  mutate(Vp_model3 = TV_vp_m3 * (WT / ref_wt)^a_wt_vp)

ggplot(wt_grid, aes(WT, Vp_model3)) +
  geom_line() +
  geom_hline(yintercept = 143, linetype = "dashed", colour = "darkred") +
  geom_hline(yintercept = TV_vp_m3, linetype = "dotted") +
  labs(x = "Body weight (kg)", y = "Peripheral volume Vp (L)",
       title = "Han 2010 Model 3: power-law WT effect on Vp",
       caption = "Dashed red line is the base-model Vp = 143 L; dotted line is the Model-3 TV(Vp) = 148 L at the 68 kg reference weight.")
```

# Assumptions and deviations

- **Base model only -- covariate sub-models reproduced in this vignette.**
  Han 2010 reports three separate single-covariate sub-models (CF, POT, WT)
  but does not publish a combined final-model parameter set. The model file
  encodes the base structural model (no covariates); the three sub-model
  equations are reproduced above using the published parameter estimates as
  post-simulation transforms (CF and POT on F as a multiplicative F-scaling
  of the base-model oral profile; WT on Vp shown as a typical-value sweep).
  This decision is consistent with the `extract-literature-model` skill's
  "replicate the author's structure" policy applied to a paper without a
  reported combined final model.

- **Diagonal IIV -- the paper estimated correlated IIV but did not publish
  the off-diagonal correlations.** Han 2010 Methods state "Correlations
  between pharmacokinetic parameters were always incorporated and
  estimated," but the off-diagonal correlations are not reported in any
  table. A diagonal IIV with the reported CV% values is used here as the
  most faithful extraction achievable without fabricating correlation
  magnitudes; downstream users who wish to add an OMEGA-BLOCK structure
  should regard the off-diagonals as unknown and either set them to zero
  (as encoded), elicit them from the corresponding author, or refit the
  block on a new dataset.

- **Residual error encoded as SDs.** Han 2010 reports
  "proportional and additive residual variability was 0.31 and 0.49 ug/mL,
  respectively." The Methods text describes the epsilon variances
  symbolically (sigma^2, sigma'^2), but the explicit ug/mL unit on the
  additive component (concentration scale rather than concentration-squared)
  and the conventional NONMEM popPK reporting practice indicate the printed
  values are standard deviations. The SD interpretation is encoded as
  `propSd = 0.31` and `addSd = 0.49`; downstream users who prefer the
  variance interpretation should set `propSd = sqrt(0.31) = 0.557` and
  `addSd = sqrt(0.49) = 0.700`.

- **CV% interpreted as the NONMEM approximate convention.** Han 2010 reports
  the IIVs as percent-CV values (e.g., 107% on CL). Following the standard
  NONMEM popPK reporting convention, the internal log-normal omega^2 is
  taken as `(CV/100)^2` (e.g., 1.07^2 = 1.1449 for CL). The exact log-normal
  identity `omega^2 = log(1 + CV^2)` (e.g., 0.7561 for CL) gives slightly
  lower values; the approximate convention is used here.

- **IV infusion modelled as a 2 h zero-order infusion to central.** Han 2010
  administered the IV doses as 2 h infusions; this is encoded explicitly in
  the simulation event table via `rate = amt / 2` on the IV dose record.
  Users who want bolus IV dosing should drop the `rate` column.

- **Cohort defined at 200 subjects.** The virtual cohort uses 200 subjects
  (versus the source paper's 13) so the per-time-point percentiles in the
  Figure 1 reproduction and the NCA summary tables are smooth. The vignette
  renders within the 5-minute pkgdown gate at this cohort size.

- **Covariate-distribution assumptions.** The virtual cohort samples body
  weight, age, and postoperative-time-on-day-of-oral-study from truncated
  normal distributions matching the Han 2010 Table 1 mean +/- SD and
  reported range. The CF indicator is sampled as a Bernoulli with
  P(CF) = 3/13 to match the source proportion. None of these covariates
  enter the base-model `model()`; they are carried only for the post-hoc
  covariate sub-model demonstrations above.

- **No combined-final-model parameter set.** Han 2010 mentions that a fully
  combined final model was built (Discussion) but reports neither its
  parameter estimates nor its OFV improvement vs the single-covariate
  sub-models. The base structural model and the three single-covariate
  sub-models above are therefore the only reproducible content of the
  paper.
