---
title: "Droperidol (Foo 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Droperidol (Foo 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Simulation size reduced for build time: this vignette took >300 s to render,
     against a 19 s median across the article set. Only the COHORT SIZE and/or
     OBSERVATION GRID DENSITY were reduced -- no model, parameter, dose or
     covariate value was touched, so the figures show the same behaviour. -->

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Foo_2016_droperidol"))$reference`
* Description: `r rxode2::rxode(readModelDb("Foo_2016_droperidol"))$description`
* Article: <https://doi.org/10.1111/bcp.13093>

# Population

Foo 2016 fit the model to a 41-patient subgroup of an Australian
emergency-department randomised controlled trial (ACTRN12607000527460)
comparing droperidol and midazolam for sedation of patients with acute
behavioural disturbance. Twenty-seven of the 41 (66%) were female, the
median age was 33 years (range 16-62), and the primary reasons for
presentation were threatened or deliberate self-harm (44%), alcohol
intoxication (41%), drug-induced delirium (10%), and psychosis (5%).
Seventeen received a single 5 mg intramuscular dose of droperidol and
24 received 10 mg; 5 of 41 received additional droperidol after the
index dose. Patient weight was not available, so allometric or linear
weight scaling could not be tested. A total of 128 plasma samples were
drawn (median 3 per patient, range 1-7) using either a
geometrically-spaced empirical design (5, 10, 30 min, 1, 2, 4, 8 h) or
a POPT-optimised design (5, 25, 40, 70 min, 2, 4, 10 h). Droperidol
was quantified by HPLC-UV with a lower limit of quantification of
5 ug/L; four observations were below this LLOQ (in four different
subjects) and were handled by the M6 method. See Foo 2016 Table 1 for
the full baseline summary; the same information is available
programmatically via `readModelDb("Foo_2016_droperidol")$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Foo_2016_droperidol.R`. The table below
collects them in one place. Table 2 of the paper duplicates the
parameter label "Vp" for the second and fifth rows; the abstract
disambiguates ("clearance of 41.9 l h-1 and volume of distribution of
the central compartment of, 73.6 l"), so the 73.6 L row is Vc and the
79.8 L row is Vp.

| Equation / parameter | Final value | Source location |
|----------------------|-------------|-----------------|
| `lka` (ka, FIXED) | 10 1/h | Table 2 (F = fixed); Results paragraph |
| `lcl` (CL) | 41.9 L/h | Table 2 (95% CI 34.8-49.0) |
| `lvc` (Vc) | 73.6 L | Table 2 row mislabelled "Vp" (95% CI 51.1-96.1); abstract confirms this is Vc |
| `lq` (Q) | 71.5 L/h | Table 2 (95% CI 42.3-100.7) |
| `lvp` (Vp) | 79.8 L | Table 2 (95% CI 58.8-100.8) |
| `etalka` (omega_ka^2, FIXED) | 1 (~100% CV) | Table 2 (F = fixed); Results paragraph |
| `etalcl` (shared IIV on CL and Vc) | 51% CV (95% CI 31.2-64.4%) | Table 2 + footnote a ("the same random effect was used for both Vc and CL") |
| `propSd` (proportional residual) | 22% CV | Table 2 (95% CI 8.5-30.3%) |
| `addSd` (additive residual, FIXED) | 0.0001 ug/L | Table 2 (F = fixed) |
| Equation: 2-compartment first-order input, first-order output | n/a | Results paragraph 2; Methods "Model building" |
| BLQ handling: M6 method | n/a | Methods "Handling data below the limit of quantitation" |

# Virtual cohort

Original observed concentrations are not publicly available. The
simulation below uses three virtual cohorts that mirror the three
dosing scenarios the published paper itself simulated (Figure 2 of
Foo 2016): a single 5 mg dose, a single 10 mg dose, and two 10 mg
doses separated by 15 min. Each cohort has 1000 simulated subjects
drawn with between-subject variability from the model's omega matrix.

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

n_per_arm <- 200
obs_grid  <- sort(unique(c(seq(0, 1, by = 0.05),
                           seq(1, 12, by = 0.1))))

make_cohort <- function(n, dose_schedule, label, id_offset) {
  doses <- dose_schedule |>
    tidyr::crossing(id = id_offset + seq_len(n)) |>
    dplyr::mutate(
      evid      = 1L,
      cmt       = "depot",
      treatment = label
    ) |>
    dplyr::select(id, time, amt, evid, cmt, treatment)
  obs <- tibble::tibble(id = id_offset + seq_len(n)) |>
    tidyr::crossing(time = obs_grid) |>
    dplyr::mutate(
      amt       = NA_real_,
      evid      = 0L,
      cmt       = "central",
      treatment = label
    )
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(
  make_cohort(
    n_per_arm,
    tibble::tibble(time = 0,             amt = 5),
    "5 mg single",          id_offset = 0L
  ),
  make_cohort(
    n_per_arm,
    tibble::tibble(time = 0,             amt = 10),
    "10 mg single",         id_offset =  n_per_arm
  ),
  make_cohort(
    n_per_arm,
    tibble::tibble(time = c(0, 0.25),    amt = c(10, 10)),
    "10 mg + 10 mg @ 15 min", id_offset = 2L * n_per_arm
  )
)

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

# Simulation

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

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

For a deterministic typical-value replication, zero out the random
effects:

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

events_typical <- events |>
  dplyr::filter(id %in% c(1L, n_per_arm + 1L, 2L * n_per_arm + 1L))

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

# Replicate published figures

## Figure 1: visual predictive check of the observed and simulated VPC

Foo 2016 Figure 1 plots the median and 5th / 95th percentiles of
observed and simulated droperidol concentrations after 5 mg and 10 mg
IM doses, with panel B zooming into the first hour. We mimic the
quantile band using the 5 mg and 10 mg single-dose cohorts.

```{r figure-1a}
band_single <- sim |>
  dplyr::filter(treatment %in% c("5 mg single", "10 mg single")) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(band_single, aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, colour = NA) +
  geom_line(linewidth = 0.8) +
  labs(
    x = "Time after dose (h)",
    y = "Droperidol concentration (ug/L)",
    title = "Replicates Figure 1A of Foo 2016",
    caption = paste(
      "Median and 5th-95th percentile band of simulated droperidol",
      "concentration over 12 h after 5 mg or 10 mg IM (n = 1000 per arm)."
    )
  )
```

```{r figure-1b}
ggplot(band_single |> dplyr::filter(time <= 1),
       aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, colour = NA) +
  geom_line(linewidth = 0.8) +
  labs(
    x = "Time after dose (h)",
    y = "Droperidol concentration (ug/L)",
    title = "Replicates Figure 1B of Foo 2016 (first hour)",
    caption = "Same simulation as Figure 1A, zoomed into the absorption phase."
  )
```

## Figure 2: probability of being above the lower limit of quantification

Foo 2016 Figure 2 reports the simulated probability that the
droperidol concentration exceeds the LLOQ (5 ug/L) over time for
three dosing scenarios. The Results paragraph states that 10 mg
single gives an 80% probability of being above LLOQ for 7 h, 5 mg
single gives ~5 h, and 10 mg + 10 mg at 15 min gives 10 h.

```{r figure-2}
LLOQ <- 5

prob_above <- sim |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(
    p_above = mean(Cc > LLOQ, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(prob_above, aes(time, p_above, colour = treatment)) +
  geom_hline(yintercept = 0.80, linetype = "dashed", colour = "grey50") +
  geom_line(linewidth = 0.9) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent_format()) +
  labs(
    x = "Time after first dose (h)",
    y = "P(Cc > 5 ug/L)",
    title = "Replicates Figure 2 of Foo 2016",
    caption = paste(
      "Probability of droperidol concentration > LLOQ (5 ug/L) by dosing",
      "scenario over 12 h (n = 1000 per arm). Dashed line at 80%."
    )
  )

duration_80 <- prob_above |>
  dplyr::group_by(treatment) |>
  dplyr::filter(p_above >= 0.80) |>
  dplyr::summarise(time_above_80 = max(time), .groups = "drop")
duration_80 |>
  dplyr::rename(
    "Treatment"                         = treatment,
    "Last time P(Cc > LLOQ) >= 80% (h)" = time_above_80
  ) |>
  knitr::kable(
  digits = 1,
  caption = paste(
    "Duration above the 80% probability of exceeding LLOQ.",
    "Paper Results: 5 mg ~5 h, 10 mg ~7 h, 10+10 mg ~10 h."
  )
)
```

## Initial (alpha) and terminal (beta) half-lives

Foo 2016 reports a median initial-phase half-life of 0.32 h
(IQR 0.26-0.37 h) and a median terminal-phase half-life of 3.0 h
(IQR 2.5-3.6 h). For a 2-compartment model with parameters at the
typical-value point estimates, the eigenvalues of the disposition
matrix give the alpha and beta rate constants directly.

```{r half-lives}
cl_t <- 41.9; vc_t <- 73.6; q_t <- 71.5; vp_t <- 79.8

kel_t <- cl_t / vc_t
k12_t <- q_t  / vc_t
k21_t <- q_t  / vp_t

a <- 1
b <- -(kel_t + k12_t + k21_t)
c <- kel_t * k21_t

lambda_alpha <- (-b + sqrt(b^2 - 4 * a * c)) / (2 * a)
lambda_beta  <- (-b - sqrt(b^2 - 4 * a * c)) / (2 * a)

t_alpha <- log(2) / lambda_alpha
t_beta  <- log(2) / lambda_beta

data.frame(
  Phase           = c("alpha", "beta"),
  rate_per_h      = c(lambda_alpha, lambda_beta),
  half_life_h     = c(t_alpha,      t_beta),
  paper_median_h  = c(0.32,         3.0)
) |>
  knitr::kable(
    digits  = c(0, 3, 2, 2),
    caption = paste(
      "Disposition half-lives from typical-value CL, Vc, Q, and Vp.",
      "Paper medians: alpha 0.32 h, beta 3.0 h."
    )
  )
```

# PKNCA validation

PKNCA-based NCA on the single-dose 5 mg and 10 mg arms.

```{r pknca}
sim_single <- sim |>
  dplyr::filter(treatment %in% c("5 mg single", "10 mg single")) |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

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

dose_df <- events |>
  dplyr::filter(treatment %in% c("5 mg single", "10 mg single"),
                evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(
  sim_single, Cc ~ time | treatment + id,
  concu = "ug/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,
  aucinf.obs  = TRUE,
  half.life   = TRUE
)

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

## Comparison against published quantities

Foo 2016 does not report a Cmax / AUC table per se, but the
Discussion enumerates simulated peak ranges and median terminal
half-lives that we can cross-check. The cells below summarise the
model-derived NCA values from the same simulations we used to
replicate Figures 1 and 2.

```{r nca_summary}
nca_tbl <- as.data.frame(nca_res$result) |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(
    median = stats::median(PPORRES, na.rm = TRUE),
    q05    = stats::quantile(PPORRES, 0.05, na.rm = TRUE),
    q95    = stats::quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  )
nca_tbl |>
  dplyr::rename(
    "Treatment"     = treatment,
    "NCA parameter" = PPTESTCD,
    "Median"        = median,
    "P05"           = q05,
    "P95"           = q95
  ) |>
  knitr::kable(
  digits  = c(0, 0, 2, 2, 2),
  caption = "Simulated NCA values across the two single-dose arms (per-subject medians and 5-95% range)."
)
```

The qualitative checkpoints from the paper are:

| Quantity (paper) | Paper value | Model value | Source |
|------------------|-------------|-------------|--------|
| 95th-percentile Cmax range, 5 mg | 20 to 110 ug/L | see `nca_tbl` above for `cmax` row, 5 mg single | Discussion paragraph "The model predicted maximum values" |
| 95th-percentile Cmax range, 10 mg | 25 to 220 ug/L | see `nca_tbl` above for `cmax` row, 10 mg single | Discussion paragraph "The model predicted maximum values" |
| Median terminal (beta) half-life | 3.0 h (IQR 2.5-3.6) | see `nca_tbl` above for `half.life` row | Results paragraph |
| Duration P(Cc > LLOQ) >= 80%, 5 mg | ~5 h | see `duration_80` table above | Results / Discussion |
| Duration P(Cc > LLOQ) >= 80%, 10 mg | ~7 h | see `duration_80` table above | Results / Discussion |
| Duration P(Cc > LLOQ) >= 80%, 10+10 mg | ~10 h | see `duration_80` table above | Results / Discussion |

Discrepancies > 20% should be investigated, not tuned. The
discussion-paragraph Cmax ranges are themselves simulation outputs of
the paper's own VPC, so an exact numeric match is the right reference.

# Assumptions and deviations

- **Table 2 typographic duplication of "Vp".** Table 2 of Foo 2016
  lists two rows labelled "Vp": the first reports 73.6 L (and shares
  its random effect with CL), and the second reports 79.8 L. The
  abstract resolves this directly ("volume of distribution of the
  central compartment of, 73.6 l"), so the 73.6 L row is Vc. The
  model file labels them as Vc and Vp accordingly.
- **Shared random effect on CL and Vc.** Table 2 footnote a states
  that the same random effect was used for both Vc and CL after the
  authors observed a 100% correlation between independent etas. The
  packaged model encodes this directly: a single `etalcl` (variance
  log(1 + 0.51^2) = 0.231) is referenced from both `cl <- exp(lcl +
  etalcl)` and `vc <- exp(lvc + etalcl)`. This is the
  mathematically-honest representation of "one shared eta" and
  avoids the singular-OMEGA / chol() failure mode that a rank-1
  block matrix encoding would create.
- **Absorption rate constant ka fixed at 10 1/h (extremely rapid
  absorption).** Foo 2016 Results: the individual `ka` estimates
  varied widely (1-140 h^-1) without a pattern; a variety of
  absorption models were unstable; fixing `ka` to 10 h^-1 with a
  fixed between-subject variance of 1 (~100% CV) gave the lowest OFV
  and a stable fit. The packaged model preserves both fixes.
- **Additive residual error fixed at 0.0001 ug/L.** The paper kept
  the additive term in the combined error model to retain numerical
  stability while letting the proportional component carry the
  observed variability; with `addSd` essentially zero, the residual
  is effectively proportional in this model.
- **No covariates retained.** Coingestion of alcohol was tested
  visually against CL and Vc and showed no association (Results).
  Patient weight was not available (Methods), so allometric or
  linear weight scaling could not be evaluated. The packaged model
  therefore contains no covariate effects, in line with Table 2.
- **BLQ handling.** Four observations from four different subjects
  were below the LLOQ of 5 ug/L and were processed by the Beal M6
  method (first BLQ in each series set to LLOQ/2, subsequent BLQ
  values commented out). The packaged model does not re-implement
  BLQ handling; the published point estimates already reflect the
  M6 fit.
- **One 16-year-old in the cohort.** Per Methods, a 16-year-old was
  inadvertently recruited because their age was unknown at the time
  of sedation. The fit therefore is for adults plus one adolescent
  rather than strictly adults; the model is otherwise applied as
  fit.
