---
title: "Midazolam (Franken 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Midazolam (Franken 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("Franken_2017_midazolam"))$reference`
* Description: `r rxode2::rxode(readModelDb("Franken_2017_midazolam"))$description`
* Article: <https://doi.org/10.1111/bcp.13259>

# Population

The published analysis included 45 terminally ill adult palliative-care patients
admitted to the Laurens Cadenza hospice (Rotterdam, The Netherlands) over a
two-year period. Median age 71 years (range 43-93), 51.1% female, 91.1%
Caucasian / 6.7% Afro-Caribbean / 2.2% Unknown; 97.8% had advanced malignancy as
the primary diagnosis (one patient with respiratory disease). Median duration of
admittance was 29 days (range 7-457). Patients received midazolam for insomnia
or palliative sedation per Dutch national guidelines, either orally up to four
times daily at 7.5 mg or subcutaneously (bolus or infusion) at 2.5-180 mg/day;
the intravenous route was not used. Sparse plasma sampling (median 2 samples
per subject, range 1-10) yielded 139 quantifiable concentrations across
midazolam, 1-OH-midazolam (1-OH-M), and 1-OH-midazolam-glucuronide (1-OH-MG).
Baseline blood-chemistry medians at admission were albumin 25 g/L (range 13-39),
creatinine 67 umol/L, eGFR by standard four-variable MDRD 104 mL/min/1.73 m^2
(range 6-328), CRP 92 mg/L. 37.8% used dexamethasone and 2.2% used phenytoin
during the observation window (concomitant CYP3A modulators). See Franken 2017
Table 1 for the full baseline summary; the same information is available
programmatically via `readModelDb("Franken_2017_midazolam")$population`.

# Source trace

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

| Equation / parameter | Final value | Source location |
|----------------------|-------------|-----------------|
| `lka_oral` (Ka oral) | 5.5 1/h (FIXED) | Methods 'Structural model' (literature [16, 22, 25]) |
| `lka_sc` (Ka SC) | 10 1/h (FIXED) | Methods 'Structural model' (literature [16, 22, 25]) |
| `lfdepot` (F oral) | 0.279 | Table 2 Final |
| `lcl` (midazolam CL at population median ALB) | 8.42 L/h | Table 2 Final |
| `lvc` (midazolam V; also V_1OH-M) | 113 L | Table 2 Final (V_1OH-M assumed equal per Methods 'Base model development') |
| `lcl_1ohm` (1-OH-M CL) | 45.4 L/h | Table 2 Final |
| `lcl_1ohmg` (1-OH-MG CL at median eGFR) | 5.10 L/h | Table 2 Final |
| `lvc_1ohmg` (1-OH-MG V) | 2.98 L | Table 2 Final (high RSE 71.5%) |
| `e_alb_cl` (albumin power on midazolam CL) | 1.08 | Table 2 Final |
| `e_crcl_cl_1ohmg` (eGFR power on 1-OH-MG CL) | 0.53 | Table 2 Final |
| `etalfdepot` + `etalcl` block (rho = 1) | 50.6% / 49.0% CV, off-diagonal sqrt(var_F * var_CL) | Table 2 Final + Methods 'Base model development' (correlation 0.93 fixed to unity) |
| `etalvc` (midazolam V) | 70.9% CV | Table 2 Final |
| `etalcl_1ohm` (1-OH-M CL) | 60.5% CV | Table 2 Final |
| `etalcl_1ohmg` (1-OH-MG CL) | 49.9% CV | Table 2 Final |
| `expSd` (midazolam residual) | sqrt(log(1 + 0.268^2)) | Table 2 Final (additive on log scale; CV% reported, converted to log-scale SD) |
| `expSd_1ohm` (1-OH-M residual) | sqrt(log(1 + 0.423^2)) | Table 2 Final |
| `expSd_1ohmg` (1-OH-MG residual) | sqrt(log(1 + 0.464^2)) | Table 2 Final |
| Equation: midazolam 1-cmt + 1-cmt 1-OH-M + 1-cmt 1-OH-MG cascade | n/a | Figure 2 schematic |
| Equation: power covariate forms theta = theta_pop * (cov / cov_m)^theta_cov | n/a | Methods Eq. 1 |
| Equation: ALB on midazolam CL, eGFR on 1-OH-MG CL | n/a | Methods Eqs. 4 and 5 |

# Virtual cohort

Original observed concentrations are not publicly available. The simulations
below use typical-value virtual subjects whose covariates match the scenarios
shown in Franken 2017 Figures 4 (albumin sweep) and 5 (eGFR sweep), plus the
combined effect verification quoted in the Discussion / Simulations section.

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

# Figure 4 cohort: vary albumin at constant eGFR = 90 mL/min/1.73 m^2.
fig4 <- data.frame(
  id    = 1:3,
  ALB   = c(15, 25, 35),
  CRCL  = 90,
  scenario = factor(
    paste0("albumin ", c(15, 25, 35), " g/L"),
    levels = paste0("albumin ", c(15, 25, 35), " g/L")
  ),
  stringsAsFactors = FALSE
)

# Figure 5 cohort: vary eGFR at constant albumin = 25 g/L.
fig5 <- data.frame(
  id    = 4:7,
  CRCL  = c(10, 30, 60, 90),
  ALB   = 25,
  scenario = factor(
    paste0("eGFR ", c(10, 30, 60, 90), " mL/min"),
    levels = paste0("eGFR ", c(10, 30, 60, 90), " mL/min")
  ),
  stringsAsFactors = FALSE
)
```

# Simulation

Both figures use the published simulation regimen from Methods 'Simulations':
a 10 mg subcutaneous midazolam loading dose followed by 5 mg subcutaneously
every 4 hours over 72 hours (six times daily). All simulations zero out the
between-subject variability and residual error via `rxode2::zeroRe()` so that
the resulting profiles are typical-value replications of the published figures.

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

build_events <- function(cov_df, obs_times, dose_events) {
  events_list <- lapply(seq_len(nrow(cov_df)), function(i) {
    row <- cov_df[i, , drop = FALSE]
    dose_rows <- dose_events |>
      dplyr::mutate(
        id       = row$id,
        ALB      = row$ALB,
        CRCL     = row$CRCL,
        scenario = as.character(row$scenario)
      )
    # One observation row per output per time point. Use the algebraic
    # observable name in cmt (rxode2 / rxUi auto-injects cmt(<observable>)
    # slots so the observable name is accepted on observation rows for
    # multi-output models).
    obs_rows <- tidyr::expand_grid(
      time = obs_times,
      cmt  = c("Cc", "Cc_1ohm", "Cc_1ohmg")
    ) |>
      dplyr::mutate(
        id       = row$id,
        amt      = NA_real_,
        rate     = NA_real_,
        evid     = 0L,
        ALB      = row$ALB,
        CRCL     = row$CRCL,
        scenario = as.character(row$scenario)
      )
    dplyr::bind_rows(dose_rows, obs_rows)
  })
  dplyr::bind_rows(events_list) |>
    dplyr::arrange(id, time, evid)
}

# 10 mg SC loading dose at t = 0 followed by 5 mg SC every 4 h through 72 h.
dose_regimen <- dplyr::bind_rows(
  data.frame(time = 0,                 amt = 10, rate = 0, evid = 1L, cmt = "depot2"),
  data.frame(time = seq(4, 72, by = 4), amt =  5, rate = 0, evid = 1L, cmt = "depot2")
)

obs_times_72 <- c(seq(0, 8, by = 0.25), seq(8.5, 72, by = 0.5))

events_fig4 <- build_events(fig4, obs_times_72, dose_regimen)
events_fig5 <- build_events(fig5, obs_times_72, dose_regimen)

sim_fig4 <- rxode2::rxSolve(
  mod_typical, events = events_fig4,
  keep = c("scenario", "ALB", "CRCL")
) |> as.data.frame()

sim_fig5 <- rxode2::rxSolve(
  mod_typical, events = events_fig5,
  keep = c("scenario", "ALB", "CRCL")
) |> as.data.frame()
```

# Replicate published figures

The paper expresses the combined "total effective concentration" as the sum
of parent and metabolite-equivalents, with 1-OH-M weighted at 0.8 and 1-OH-MG
at 0.1 to reflect their relative midazolam potencies (Figure 4 / Figure 5
captions). The plots below show the three species separately so the
underlying disposition is visible, plus the combined effective concentration
as a fourth panel.

## Figure 4: effect of serum albumin on midazolam and metabolite concentrations

```{r figure-4}
plot_fig4 <- sim_fig4 |>
  dplyr::mutate(
    Effective = Cc + 0.8 * Cc_1ohm + 0.1 * Cc_1ohmg
  ) |>
  tidyr::pivot_longer(
    cols      = c(Cc, Cc_1ohm, Cc_1ohmg, Effective),
    names_to  = "species",
    values_to = "conc"
  ) |>
  dplyr::mutate(species = factor(
    dplyr::recode(species,
                  "Cc"        = "Midazolam",
                  "Cc_1ohm"   = "1-OH-M",
                  "Cc_1ohmg"  = "1-OH-MG",
                  "Effective" = "Effective (Mid + 0.8*1-OH-M + 0.1*1-OH-MG)"),
    levels = c("Midazolam", "1-OH-M", "1-OH-MG",
               "Effective (Mid + 0.8*1-OH-M + 0.1*1-OH-MG)")
  ))

ggplot(plot_fig4, aes(x = time, y = conc, color = scenario)) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~species, scales = "free_y", ncol = 2) +
  labs(x = "Time (h)",
       y = "Concentration (ug/L; midazolam equivalents)",
       color = "Albumin scenario",
       title = "Replicates Figure 4 of Franken 2017",
       caption = paste(
         "10 mg subcutaneous midazolam loading dose then 5 mg SC every 4 h x 72 h.",
         "eGFR held at 90 mL/min/1.73 m^2.",
         "Typical-value simulation (no IIV / no RUV)."
       ))
```

## Figure 5: effect of eGFR on midazolam and metabolite concentrations

```{r figure-5}
plot_fig5 <- sim_fig5 |>
  dplyr::mutate(
    Effective = Cc + 0.8 * Cc_1ohm + 0.1 * Cc_1ohmg
  ) |>
  tidyr::pivot_longer(
    cols      = c(Cc, Cc_1ohm, Cc_1ohmg, Effective),
    names_to  = "species",
    values_to = "conc"
  ) |>
  dplyr::mutate(species = factor(
    dplyr::recode(species,
                  "Cc"        = "Midazolam",
                  "Cc_1ohm"   = "1-OH-M",
                  "Cc_1ohmg"  = "1-OH-MG",
                  "Effective" = "Effective (Mid + 0.8*1-OH-M + 0.1*1-OH-MG)"),
    levels = c("Midazolam", "1-OH-M", "1-OH-MG",
               "Effective (Mid + 0.8*1-OH-M + 0.1*1-OH-MG)")
  ))

ggplot(plot_fig5, aes(x = time, y = conc, color = scenario)) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~species, scales = "free_y", ncol = 2) +
  labs(x = "Time (h)",
       y = "Concentration (ug/L; midazolam equivalents)",
       color = "eGFR scenario",
       title = "Replicates Figure 5 of Franken 2017",
       caption = paste(
         "10 mg subcutaneous midazolam loading dose then 5 mg SC every 4 h x 72 h.",
         "Albumin held at 25 g/L.",
         "Typical-value simulation (no IIV / no RUV)."
       ))
```

# Numeric checks against the Discussion / Simulations section

Franken 2017 'Simulations' quotes four numeric anchor results from the final
model. The model file's `ini()` reproduces these closed-form, so the table
below also serves as a direct provenance audit for the two covariate
coefficients.

```{r numeric-checks}
typical_CL <- function(theta_pop, cov, cov_ref, exponent) {
  theta_pop * (cov / cov_ref)^exponent
}

anchor <- tibble::tribble(
  ~quantity,                       ~scenario,                    ~paper,  ~model,
  "Midazolam CL (L/h)",            "ALB = 35 g/L (ref 25, exp 1.08)",
                                                                  12.1,
                                                                  typical_CL(8.42, 35, 25, 1.08),
  "Midazolam CL (L/h)",            "ALB = 25 g/L (reference)",     8.4,
                                                                  typical_CL(8.42, 25, 25, 1.08),
  "Midazolam CL (L/h)",            "ALB = 15 g/L",                 4.8,
                                                                  typical_CL(8.42, 15, 25, 1.08),
  "Midazolam CL (L/h)",            "ALB = 45 g/L (healthy)",       15.3,
                                                                  typical_CL(8.42, 45, 25, 1.08),
  "1-OH-MG CL (L/h)",              "eGFR = 90 (ref 104, exp 0.53)",
                                                                   4.7,
                                                                   typical_CL(5.10, 90, 104, 0.53),
  "1-OH-MG CL (L/h)",              "eGFR = 50",                    3.5,
                                                                   typical_CL(5.10, 50, 104, 0.53),
  "1-OH-MG CL (L/h)",              "eGFR = 30",                    2.7,
                                                                   typical_CL(5.10, 30, 104, 0.53)
) |>
  dplyr::mutate(
    diff_pct = round(100 * (model - paper) / paper, 1)
  )

knitr::kable(
  anchor,
  caption = paste(
    "Typical-value covariate-anchored CL predictions vs Franken 2017",
    "Discussion / Simulations text. Sub-3% differences trace to rounding",
    "of the published coefficients (1.08 and 0.53) and reference values",
    "(median ALB 25 g/L, median eGFR 104 mL/min/1.73 m^2)."
  ),
  align = c("l", "l", "r", "r", "r")
)
```

# PKNCA validation

Single-dose NCA after the first 10 mg subcutaneous midazolam loading dose
(window 0-4 h, before the second dose), computed separately for each analyte
and grouped by the eGFR scenarios of Figure 5. This exercises the parent +
both metabolites within a single NONMEM-style cohort. The paper itself does
not report a baseline NCA table (sparse-sampling design with no IV reference
arm), so this section validates that PKNCA runs cleanly on the simulated
profiles and surfaces the dose-proportional response in midazolam Cmax with
no dependence on eGFR (as expected -- eGFR only affects 1-OH-MG CL).

```{r pknca-setup}
sim_nca_window <- sim_fig5 |> dplyr::filter(time <= 4)

dose_df_fig5 <- events_fig5 |>
  dplyr::filter(evid == 1L, time == 0) |>
  dplyr::select(id, time, amt, scenario)

intervals_sd <- data.frame(
  start     = 0,
  end       = 4,
  cmax      = TRUE,
  tmax      = TRUE,
  auclast   = TRUE,
  half.life = FALSE
)
```

## Midazolam

```{r pknca-midazolam}
# Guarantee a time = 0 anchor row per (id, scenario) for PKNCA AUC0-*
add_t0 <- function(df, value_col) {
  df_t0 <- df |>
    dplyr::distinct(id, scenario) |>
    dplyr::mutate(time = 0, !!value_col := 0)
  dplyr::bind_rows(df, df_t0) |>
    dplyr::distinct(id, scenario, time, .keep_all = TRUE) |>
    dplyr::arrange(id, scenario, time)
}

sim_nca_midaz <- sim_nca_window |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, scenario) |>
  add_t0("Cc")

conc_obj_midaz <- PKNCA::PKNCAconc(
  sim_nca_midaz, Cc ~ time | scenario + id,
  concu = "ug/L", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df_fig5, amt ~ time | scenario + id, doseu = "mg"
)

nca_midaz <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_midaz, dose_obj, intervals = intervals_sd)
)
knitr::kable(as.data.frame(summary(nca_midaz)),
             caption = "Midazolam NCA over 4 h after the first 10 mg SC bolus.")
```

## 1-OH-M

```{r pknca-1ohm}
sim_nca_1ohm <- sim_nca_window |>
  dplyr::filter(!is.na(Cc_1ohm)) |>
  dplyr::select(id, time, Cc_1ohm, scenario) |>
  dplyr::rename(Cc = Cc_1ohm) |>
  add_t0("Cc")

conc_obj_1ohm <- PKNCA::PKNCAconc(
  sim_nca_1ohm, Cc ~ time | scenario + id,
  concu = "ug/L", timeu = "h"
)
nca_1ohm <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_1ohm, dose_obj, intervals = intervals_sd)
)
knitr::kable(as.data.frame(summary(nca_1ohm)),
             caption = "1-OH-M NCA over 4 h after the first 10 mg SC midazolam bolus.")
```

## 1-OH-MG

```{r pknca-1ohmg}
sim_nca_1ohmg <- sim_nca_window |>
  dplyr::filter(!is.na(Cc_1ohmg)) |>
  dplyr::select(id, time, Cc_1ohmg, scenario) |>
  dplyr::rename(Cc = Cc_1ohmg) |>
  add_t0("Cc")

conc_obj_1ohmg <- PKNCA::PKNCAconc(
  sim_nca_1ohmg, Cc ~ time | scenario + id,
  concu = "ug/L", timeu = "h"
)
nca_1ohmg <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_1ohmg, dose_obj, intervals = intervals_sd)
)
knitr::kable(as.data.frame(summary(nca_1ohmg)),
             caption = "1-OH-MG NCA over 4 h after the first 10 mg SC midazolam bolus.")
```

# Assumptions and deviations

- **Residual-error interpretation.** Franken 2017 Methods 'Base model
  development' reports residual variability as "additive residual error on
  logarithmic transformed concentrations" (LTBS), with Table 2 Final values
  reported as CV percentages (26.8% / 42.3% / 46.4% for midazolam / 1-OH-M /
  1-OH-MG). The packaged model converts each CV via the canonical log-normal
  identity `omega^2 = log(1 + CV^2)` and passes the square root to rxode2's
  `lnorm()` residual model. The Franken 2015 morphine paper from the same
  group reports residual variability as raw NONMEM `$SIGMA` variance values
  (no percent sign); the difference in reporting style is paper-specific.
  Both interpretations of the midazolam Table 2 values (CV% -> log-normal
  conversion versus log-scale SD directly) give numerically indistinguishable
  predictions for CV <= 0.5.
- **Cross-output residual correlation not encoded.** Franken 2017 Methods
  'Base model development' notes that "since the parent and metabolite
  concentrations were measured in the same samples using a single assay, a
  correlation between the residual errors was incorporated in the model."
  The paper does not report the correlation magnitude in Table 2, and
  rxode2 / nlmixr2 do not support cross-output residual correlations
  directly via the `Cc ~ lnorm(expSd)` interface. The packaged model treats
  the three residuals as independent, which slightly underestimates the
  joint shrinkage of simulated parent-vs-metabolite trajectories at any
  single time point. Marginal univariate VPCs and per-output NCA summaries
  are unaffected.
- **1-OH-M volume of distribution.** The packaged model sets `vc_1ohm = vc`
  (a structural equality, not a separately estimated parameter), matching
  Methods 'Base model development': "The volume of distribution (V) of
  1-OH-M was assumed to be equal to the volume of distribution of
  midazolam."
- **Apparent vs absolute clearance.** Per Discussion: "the 4-OH-M metabolite
  was not included in the model and the true fraction metabolized is
  unknown, the clearance and distribution volumes of the other metabolites
  are apparent values." The packaged model uses the published apparent
  values directly without attempting to back out fractions; downstream
  users who want absolute clearance should rescale by their own assumed
  fraction-metabolised.
- **Absorption rate constants fixed from literature.** Methods 'Structural
  model': "the absorption constants (Ka) could not be estimated. They were
  therefore derived from literature (5.5 h^-1 for oral administration, 10
  h^-1 for subcutaneous injection)." Wrapped in `fixed()` in the model
  file; both values are structural assumptions, not estimated.
- **SC bioavailability assumed 1.** Methods 'Base model development':
  "Bioavailability of subcutaneous midazolam was assumed to be 100%." The
  packaged model leaves the SC depot's bioavailability at the rxode2
  default 1 (no `f(depot2) <- ...` assignment) and estimates oral F via
  `lfdepot`.
- **F-vs-CL correlation fixed to unity.** Methods 'Base model development':
  "The correlation between IIV of midazolam clearance and F of oral
  midazolam was high (0.93) and therefore fixed to unity." Encoded as a
  2x2 OMEGA block with off-diagonal `sqrt(var_F * var_CL)`, which is
  mathematically equivalent to the published shared-eta-with-scaling-theta
  NONMEM construct.
- **eGFR formula.** The packaged model uses the standard four-variable
  MDRD eGFR (paper's "eGFR-a", Methods Eq. 4), which is the equation
  retained in the Final model. The original six-variable MDRD formula
  ("eGFR-b", Eq. 5) was tested but not retained because the difference in
  fit was minimal.
- **Time-to-death (TTD) not retained.** TTD was tested as a covariate on
  midazolam CL (delta-OFV = -4.69), midazolam V (delta-OFV = -6.50), and
  1-OH-MG CL (delta-OFV = -15.25) per Table 3, but did not survive
  backward elimination at p < 0.001. The packaged model omits TTD
  accordingly. The Franken 2015 companion morphine model
  (`modellib("Franken_2015_morphine")`) retains TTD on parent CL because
  it did clear the morphine paper's stricter selection process.
- **Bootstrap-derived 95% confidence intervals not surfaced.** Table 2
  reports bootstrap (200-run) 95% CI bounds for every estimate. The model
  file carries only the Final-model point estimates; the bootstrap CIs
  are documented in the paper's Table 2 and the in-file comments
  reference them implicitly via the RSE column.
- **NPDE not reproduced.** The paper validated the final model via NPDE
  (global adjusted p = 0.75 / 0.20 / 0.41 for the three analytes), not
  VPC, because of the variable dosing-regimen design that precludes a
  standard VPC. The packaged vignette replicates the typical-value
  Figure 4 / Figure 5 simulations (the only quantitative figures the
  paper presents) rather than re-running NPDE.
