---
title: "Remifentanil (Yang 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Remifentanil (Yang 2017)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{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("Yang_2017_remifentanil"))$reference`
* Description: `r rxode2::rxode(readModelDb("Yang_2017_remifentanil"))$description`
* Article: <https://doi.org/10.1038/s41598-017-16358-6> (open access in Scientific Reports 2017;7:16275)

# Population

Yang et al. enrolled 15 critically ill adults (median age 57 years, IQR
45-69; 67% male; median weight 65.4 kg, IQR 54.5-70.0; median BMI
23.8 kg/m^2) receiving continuous-infusion remifentanil sedation during
venoarterial extracorporeal membrane oxygenation (VA-ECMO) at Severance
Cardiovascular Hospital, Yonsei University (Seoul, South Korea), between
January 2015 and December 2016 (ClinicalTrials.gov NCT02581280;
single-center prospective cohort). Indications for VA-ECMO were
predominantly acute coronary syndromes (acute MI n = 5, NSTEMI n = 4,
STEMI n = 1, ischemic cardiomyopathy n = 1) with several other
cardiopulmonary failures (pulmonary embolism, coronary artery occlusive
disease, myocarditis, atrial fibrillation with bronchiolitis, angina
pectoris). Median ECMO support duration was 143 h (IQR 96-250). Median
centrifugal-pump rotational speed was 2350 RPM (IQR 2302-2532). Ten of
the 15 patients (67%) received concomitant CRRT during ECMO support;
CRRT was tested as a covariate during model development but not retained
in the final model. Baseline characteristics are listed in Table 1 of
the source paper. The same information is available programmatically via
`readModelDb("Yang_2017_remifentanil")()$meta$population`.

# Source trace

Per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Yang_2017_remifentanil.R`. The table below
collects the equations and parameter values in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| One-compartment model with zero-order input and first-order elimination | n/a | Results, paragraph 2; Methods, Structural and model development |
| Final CL model: CL (L/h) = 366 x 0.502^SEX x (ECMO pump speed / 2350)^2.04, with SEX = 0 (female), 1 (male) | n/a | Results, Population PK analysis (final-model equation) |
| `lcl` (typical female CL at 2350 RPM) | 366 L/h | Table 2 (theta_CL) |
| `lvc` (typical V) | 41 L | Table 2 (theta_V) |
| `e_sex_cl` (log(0.502); applied via (1 - SEXF) to preserve female reference) | -0.6892 | Table 2 (theta_sex on CL = 0.502) |
| `e_pumpspeed_cl` (power exponent on ECMO pump speed / 2350) | 2.04 | Table 2 (theta_ECMOpumpspeed on CL) |
| `etalcl` (IIV variance on log CL; stored as SD^2 = 0.124^2) | 0.01538 | Table 2 (omega_CL = 0.124, SD on log scale) |
| `propSd` (proportional residual SD, fraction) | 0.387 | Table 2 (sigma_proportional = 0.387) |
| `addSd` (additive residual SD, ng/mL) | 0.111 ng/mL | Table 2 (sigma_additive = 0.111 ng/mL) |
| Reference pump speed (median of cohort) | 2350 RPM | Results, Patient characteristics; Table 1 |

# Virtual cohort

Original observed concentrations are not publicly available. The
figures below use a virtual cohort that mirrors the simulation grid
used in the source paper (Methods, Simulations; Results, Predicted
concentration profiles; Figure 4):

* Sex: female (`SEXF = 1`) and male (`SEXF = 0`).
* ECMO pump speed: 1700, 2000, 2300, 2600, and 2900 RPM (Yang 2017
  Methods, Simulations).
* Remifentanil infusion rates: 0.105, 0.21, 0.42, 0.63, and 0.84 mg/h
  (= 105, 210, 420, 630, 840 ug/h). The infusion is continuous over
  48 h, matching the paper's two-day simulation horizon.

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

pump_speeds_rpm <- c(1700, 2000, 2300, 2600, 2900)
infusion_rates_mg_h <- c(0.105, 0.21, 0.42, 0.63, 0.84)
sexes <- c(Female = 1L, Male = 0L)

n_per_group <- 25L             # virtual subjects per (sex, pump speed, dose) group
inf_duration_h <- 48           # 2-day continuous infusion (Methods, Simulations)
obs_times <- seq(0, 48, by = 2)

make_cohort <- function(n, sex_label, sexf_val, pump_rpm, dose_mg_h, id_offset) {
  rate_ug_h <- dose_mg_h * 1000  # mg/h -> ug/h (model dose units are ug)
  amt_ug <- rate_ug_h * inf_duration_h

  per_subject <- function(sid) {
    dplyr::bind_rows(
      data.frame(id = sid, time = 0, amt = amt_ug, rate = rate_ug_h,
                 evid = 1L, cmt = "central",
                 SEXF = sexf_val,
                 ECMO_PUMP_SPEED = pump_rpm,
                 sex_label = sex_label,
                 pump_rpm = pump_rpm,
                 dose_mg_h = dose_mg_h),
      data.frame(id = sid, time = obs_times, amt = NA_real_, rate = NA_real_,
                 evid = 0L, cmt = NA_character_,
                 SEXF = sexf_val,
                 ECMO_PUMP_SPEED = pump_rpm,
                 sex_label = sex_label,
                 pump_rpm = pump_rpm,
                 dose_mg_h = dose_mg_h)
    )
  }

  ids <- id_offset + seq_len(n)
  do.call(rbind, lapply(ids, per_subject))
}

# Build a grid: 2 sex x 5 pump speeds x 5 doses = 50 groups
grid <- tidyr::expand_grid(
  sex_label = names(sexes),
  pump_rpm = pump_speeds_rpm,
  dose_mg_h = infusion_rates_mg_h
)
grid$id_offset <- (seq_len(nrow(grid)) - 1L) * n_per_group

events <- do.call(rbind, lapply(seq_len(nrow(grid)), function(i) {
  row <- grid[i, ]
  make_cohort(
    n = n_per_group,
    sex_label = row$sex_label,
    sexf_val = sexes[[row$sex_label]],
    pump_rpm = row$pump_rpm,
    dose_mg_h = row$dose_mg_h,
    id_offset = row$id_offset
  )
}))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
events$cohort <- paste0(events$sex_label,
                        " / ", events$pump_rpm, " RPM / ",
                        events$dose_mg_h, " mg/h")
```

# Simulation

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

sim <- rxode2::rxSolve(
  mod, events = events, addCov = TRUE,
  keep = c("sex_label", "pump_rpm", "dose_mg_h", "cohort")
)
sim <- as.data.frame(sim)
```

# Replicate Figure 4 -- simulated remifentanil concentrations by sex x pump speed x dose

The paper's Figure 4 shows simulated mean remifentanil concentrations
in female and male patients across the 5 pump speeds (1700, 2000, 2300,
2600, 2900 RPM) for 5 infusion rates (0.84, 0.63, 0.42, 0.21,
0.105 mg/h). The horizontal reference line at 1.5 ng/mL is the
sedation target referenced in the paper (Discussion; "a target
concentration >= 1.5 ng/mL").

```{r figure-4}
fig4 <- sim |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::group_by(sex_label, pump_rpm, dose_mg_h, time) |>
  dplyr::summarise(
    median_Cc = median(Cc, na.rm = TRUE),
    p05 = quantile(Cc, 0.05, na.rm = TRUE),
    p95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(fig4, aes(time, median_Cc, colour = factor(pump_rpm))) +
  geom_hline(yintercept = 1.5, linetype = "dashed", colour = "grey50") +
  geom_line() +
  facet_grid(dose_mg_h ~ sex_label, scales = "free_y",
             labeller = labeller(dose_mg_h = function(x) paste0(x, " mg/h"))) +
  scale_colour_brewer("Pump speed (RPM)", palette = "Set1") +
  labs(
    x = "Time (h)",
    y = "Simulated remifentanil concentration (ng/mL)",
    title = "Median simulated Cc by sex x ECMO pump speed x infusion rate",
    caption = "Replicates Figure 4 of Yang 2017; dashed line = 1.5 ng/mL sedation target."
  ) +
  theme_bw()
```

# Comparison against the published dosing recommendation

Yang 2017 reports (Results, Predicted concentration profiles):

* For female patients: pump speed 1700-2000 RPM, >= 0.42 mg/h;
  pump speed 2000-2900 RPM, >= 0.63 mg/h achieves 95% of subjects
  maintaining Cc >= 1.5 ng/mL.
* For male patients: pump speed 1700-2000 RPM, >= 0.21 mg/h;
  pump speed 2000-2900 RPM, >= 0.42 mg/h.

Reproduce the same threshold check from the simulated cohort by
computing, per (sex, pump speed, dose) cell, the fraction of subjects
whose steady-state Cc (taken at t = 48 h) meets or exceeds 1.5 ng/mL.

```{r dose-target-check}
target_check <- sim |>
  dplyr::filter(!is.na(Cc), time == 48) |>
  dplyr::group_by(sex_label, pump_rpm, dose_mg_h) |>
  dplyr::summarise(
    pct_at_target = mean(Cc >= 1.5, na.rm = TRUE) * 100,
    median_Cc_48h = median(Cc, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(
  target_check |>
    dplyr::mutate(
      pct_at_target = sprintf("%.0f%%", pct_at_target),
      median_Cc_48h = sprintf("%.2f ng/mL", median_Cc_48h)
    ) |>
    tidyr::pivot_wider(
      names_from = pump_rpm,
      values_from = c(pct_at_target, median_Cc_48h)
    ),
  caption = "Per-cell percentage of virtual subjects at or above 1.5 ng/mL at t = 48 h, and median Cc at t = 48 h. Compare with Yang 2017 dosing recommendation (Results, Predicted concentration profiles)."
)
```

# PKNCA validation

PKNCA NCA on the 0-48 h continuous-infusion profile, grouped by sex,
pump speed, and dose. The treatment grouping comes before `id` per
project convention.

For closed-form sanity-checks at steady state (typical-value, no IIV):

* Steady-state concentration: Css = (infusion rate) / CL_i, where
  CL_i = 366 x 0.502^SEX x (ECMO_PUMP_SPEED / 2350)^2.04 L/h.
* Elimination half-life: t_(1/2) = ln(2) x V / CL_i. With V = 41 L:

```{r closed-form-css}
closed_form <- tidyr::expand_grid(
  sex_label = names(sexes),
  pump_rpm = pump_speeds_rpm,
  dose_mg_h = infusion_rates_mg_h
) |>
  dplyr::mutate(
    sex_male_indicator = ifelse(sex_label == "Male", 1, 0),
    CL_Lh = 366 * 0.502^sex_male_indicator * (pump_rpm / 2350)^2.04,
    Css_ng_mL = (dose_mg_h * 1000) / CL_Lh,
    t_half_min = log(2) * 41 / CL_Lh * 60
  )

knitr::kable(
  closed_form |>
    dplyr::filter(dose_mg_h == 0.35) |>
    dplyr::transmute(
      Sex = sex_label,
      `Pump speed (RPM)` = pump_rpm,
      `CL (L/h)` = round(CL_Lh, 1),
      `Css at 0.35 mg/h (ng/mL)` = round(Css_ng_mL, 2),
      `t_(1/2) (min)` = round(t_half_min, 2)
    ),
  caption = "Closed-form typical-value CL, Css, and t_(1/2) for the cohort-median infusion rate (0.35 mg/h)."
)
```

PKNCA on a manageable subset of the simulation -- female / male at the
cohort-median pump speed (2300 RPM, nearest grid point to the cohort
median of 2350 RPM) and a representative subset of doses:

```{r pknca}
sim_for_nca <- sim |>
  dplyr::filter(pump_rpm == 2300, dose_mg_h %in% c(0.21, 0.42, 0.84), !is.na(Cc)) |>
  dplyr::mutate(
    treatment = paste0(sex_label, "/", dose_mg_h, "mg/h")
  )

sim_nca <- sim_for_nca |>
  dplyr::select(id, time, Cc, treatment)

dose_df <- events |>
  dplyr::filter(evid == 1, pump_rpm == 2300,
                dose_mg_h %in% c(0.21, 0.42, 0.84)) |>
  dplyr::transmute(
    id = id, time = time, amt = amt,
    treatment = paste0(sex_label, "/", dose_mg_h, "mg/h"),
    rate = rate, duration = amt / rate
  )

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)
dose_obj <- PKNCA::PKNCAdose(
  dose_df,
  amt ~ time | treatment + id,
  duration = "duration"
)

intervals <- data.frame(
  start    = 0,
  end      = 48,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE,
  cav      = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)
knitr::kable(summary(nca_res),
  caption = "NCA over the 0-48 h continuous-infusion window for female and male subjects at 2300 RPM and three infusion rates. Cav at steady state equals the closed-form Css = (infusion rate) / CL_i."
)
```

# Assumptions and deviations

* **Variance-vs-SD reporting in Table 2.** Yang 2017 Methods defines
  the inter-individual variance as omega^2 and the residual variance
  as sigma^2, while Table 2 column headers report `omega CL`,
  `sigma proportional`, and `sigma additive` without superscripts. The
  numeric values are interpreted here as the SDs (i.e. the square root
  of the variances defined in Methods), based on three pieces of
  evidence: (1) the Greek letter without a `^2` superscript is the
  standard mathematical convention for the SD; (2) Table 2 lists
  `sigma additive = 0.111 ng/mL` in linear concentration units, whereas
  a variance would carry concentration-squared units (the additive
  variance interpretation 0.111 (ng/mL)^2 would correspond to SD =
  0.333 ng/mL ~= 6.6 x the assay LLOQ of 0.05 ng/mL, larger than
  expected for a validated LC-MS/MS assay); and (3) this matches the
  convention used in the `Abboud_2009_epinephrine` model file, whose
  source paper has an analogous "variance is omega^2, table reports the
  SD" structure. The `ini()` block accordingly stores
  `etalcl ~ 0.124^2 = 0.01538` (variance, derived from SD = 0.124) and
  the residual error parameters `propSd = 0.387` and
  `addSd = 0.111 ng/mL` as SDs. If a future reading of the original
  NONMEM control stream shows the Table 2 values were the variances
  themselves, the encoding would change to
  `etalcl ~ 0.124` (variance), `propSd = sqrt(0.387) = 0.622`, and
  `addSd = sqrt(0.111) = 0.333 ng/mL`.
* **Sex encoded under canonical SEXF with an inverted reference.**
  Yang 2017 uses the male-indicator convention (SEX = 0 female, SEX =
  1 male) with female as the structural reference (the published
  equation `CL = 366 x 0.502^SEX` yields CL = 366 L/h for SEX = 0 /
  female). To store under the canonical SEXF (1 = female, 0 = male)
  while preserving the paper's female-reference CL of 366 L/h, the
  effect is applied in `model()` as `exp(e_sex_cl * (1 - SEXF))` with
  `e_sex_cl = log(0.502)`. This is the same convention as
  `Bajaj_2017_nivolumab`.
* **Sex effect should be re-interpreted cautiously.** The original
  paper concludes (Discussion) that "there is no conclusive
  information indicating that true sex differences exist in the PK of
  remifentanil due to small sample size and uneven sex ratio in the
  present study", with 10 male and 5 female patients in the cohort.
  The covariate is preserved here because it is part of the published
  final model, but downstream users should interpret simulated male
  vs female differences with the same caveat as the paper authors.
* **ECMO pump speed as a continuous covariate.** Tested as
  continuous in the source paper (Methods, Covariate model) using
  power, exponential, and linear forms, with the power form retained.
  The simulation grid here spans 1700-2900 RPM (Methods, Simulations);
  extrapolating outside this range is not supported by the source
  data.
* **CRRT tested but not retained.** The source paper tested CRRT
  presence/absence (10/15 patients on CRRT during ECMO support) as a
  covariate and did not retain it in the final model. The library
  model therefore does not include RRT_CRRT_STATUS as an effect; the
  population metadata records the cohort CRRT prevalence (67%) for
  reference.
* **No bolus dosing.** The source cohort received remifentanil
  exclusively by continuous infusion (Methods, Study procedures: "no
  patients required the bolus injection of remifentanil during the
  study"). The library model has no depot / first-order absorption
  compartment; the event table dosing convention is to drive the
  central compartment with `evid = 1`, `amt = rate x duration`, and
  `rate > 0`.
* **No IIV on V.** The source paper estimated IIV on V as near zero
  due to the narrow weight range in the 15-patient cohort (Results,
  Population PK analysis), so the final model retains IIV only on CL.
  The library model follows the same convention; `vc <- exp(lvc)` has
  no associated eta.
* **Pump speed treated as time-fixed.** Yang 2017 reports a per-
  patient ECMO pump speed (the median over the PK sampling window)
  and does not resolve session-level pump-speed changes. The
  simulation event table accordingly carries `ECMO_PUMP_SPEED` as a
  per-subject covariate. For prospective designs that resolve
  session-level pump-speed adjustments, the covariate is naturally
  time-varying and the event table would carry one row per time
  point.
