---
title: "Morphine (Elkomy 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Morphine (Elkomy 2015)}
  %\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("Elkomy_2015_morphine"))$reference`
* Description: `r rxode2::rxode(readModelDb("Elkomy_2015_morphine"))$description`
* Article: <https://doi.org/10.1208/s12248-015-9826-5>

# Population

Elkomy 2015 enrolled 20 infants and young children aged 3 days to 5.4 years
(median 1.4 years), body weight 3.1 to 18.5 kg (median 7.8 kg, mean 7.8 +/- 4.3
kg), 10 male / 10 female, admitted to the cardiovascular intensive care unit
after congenital heart surgery at two US centres (Lucile Packard Children's
Hospital, Palo Alto, CA; Children's Hospital Colorado, Aurora, CO). Cardiac
diagnoses were tetralogy of Fallot (12), atrioventricular septal defect (3),
ventricular septal defect (1), and other (4). Renal function (Schwartz eGFR)
ranged 19 to 158 mL/min/1.73 m^2 (mean 79). Single-ventricle physiology and
significant renal / hepatic / neurologic impairment were excluded. Patients
received an IV bolus loading dose of 0.15 mg/kg morphine sulfate and
subsequent nurse-controlled analgesia (NCA) IV boluses averaging 0.06 mg/kg
(range 0.02 to 0.21). The model used 1080 concentration observations across
the 20 subjects. See Elkomy 2015 Table I for the full baseline summary; the
same information is available programmatically via
`readModelDb("Elkomy_2015_morphine")$population`.

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (theta_CL) | 8.7 L/h | Table II Fixed effects |
| `lvc` (theta_VC) | 18.7 L | Table II Fixed effects |
| `lq`  (theta_CLD) | 0.88 L/h | Table II Fixed effects |
| `lvp` (theta_VP) | 25.2 L | Table II Fixed effects |
| `e_wt_cl_q` (allometric exponent on CL, CLD) | 0.75 FIXED | Methods para 4 / Results para 2 |
| `e_wt_vc_vp` (allometric exponent on VC, VP) | 1.0 FIXED | Methods para 4 / Results para 2 |
| Reference WT (allometric normalisation) | 6 kg (study median) | Methods para 4 |
| `lkeo_m3g` (theta_Kint,M3G at GFR = 70) | 0.30 1/h | Table II Fixed effects |
| `lkeo_m6g` (theta_Kint,M6G at GFR = 70) | 0.24 1/h | Table II Fixed effects |
| GFR exponent on Kint | +1 FIXED | Results, final paragraph |
| `lemax_m3g` (theta_Mmax,M3G at GFR = 70) | 1257 nM | Table II Fixed effects |
| `lemax_m6g` (theta_Mmax,M6G at GFR = 70) | 109 nM | Table II Fixed effects |
| GFR exponent on Mmax | -1 FIXED | Results, final paragraph |
| `lec50_m3g` (Cint,50,M3G) | 40.4 nM | Table II Fixed effects |
| `lec50_m6g` (Cint,50,M6G) | 21.9 nM | Table II Fixed effects |
| Reference GFR (renal normalisation) | 70 mL/min/1.73 m^2 | Methods Simulation; Schwartz formula reference value |
| `etalcl`  (IIV variance on CL)  | 0.20 | Table II IIV column |
| `etalvc`  (IIV variance on VC)  | 0.67 | Table II IIV column |
| `etalq`   (IIV variance on CLD) | 0.17 | Table II IIV column |
| `etalvp`  (IIV variance on VP)  | 0.42 | Table II IIV column |
| `etalkeo`  (shared IIV on Kint, M3G + M6G)  | 0.27 | Table II IIV column |
| `etalemax` (shared IIV on Mmax, M3G + M6G)  | 0.23 | Table II IIV column |
| `etalec50` (shared IIV on Cint50, M3G + M6G) | 1.26 | Table II IIV column |
| `propSd`     (morphine prop residual SD) | sqrt(0.32) | Table II sigma^2(morphine) |
| `propSd_m3g` (M3G prop residual SD)      | sqrt(0.32) | Table II sigma^2(M3G) |
| `propSd_m6g` (M6G prop residual SD)      | sqrt(0.39) | Table II sigma^2(M6G) |
| Equation: morphine 2-cmt IV PK           | n/a | Methods para 1; Fig. 1 schematic |
| Equation: dCint/dt = Kint * (Cc - Cint)  | Eq. 1 | Methods para 2 |
| Equation: M = Mmax * Cint / (Cint50 + Cint) | Eq. 2 | Methods para 2 |
| Dose conversion: 1 mg morphine sulfate -> 2 * 1e6 / 758.83 nmol free morphine | n/a | Methods para 1 ("morphine sulfate molecule comprises two morphine molecules") |

# Virtual cohort

Original observed data are not publicly available. The simulation below mirrors
the published Simulation section: three age groups for morphine exposure
(group I, 1-month-old, 4.5 kg; group II, 6-month-old, 6 kg; group III,
2-year-old, 12 kg) and four GFR groups for glucuronide exposure in a 6-month-
old child (20, 50, 70, 100 mL/min/1.73 m^2). Each group has 100 subjects in
the paper; here we mirror that sampling size for the box-plot replicate.

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

# Morphine-exposure age groups (paper Fig. 7 upper row). GFR is set to the
# typical pediatric value (70 mL/min/1.73 m^2) so the morphine-channel
# comparison isolates the WT effect.
morphine_groups <- tibble::tibble(
  group = factor(c("1 mo (4.5 kg)", "6 mo (6 kg)", "2 yr (12 kg)"),
                 levels = c("1 mo (4.5 kg)", "6 mo (6 kg)", "2 yr (12 kg)")),
  WT    = c(4.5, 6, 12),
  CRCL  = c(70, 70, 70)
)

# Glucuronide-exposure GFR groups (paper Fig. 7 lower row): 6-month-old
# typical 6-kg child across four renal-function strata.
gfr_groups <- tibble::tibble(
  group = factor(paste0("GFR ", c(20, 50, 70, 100)),
                 levels = paste0("GFR ", c(20, 50, 70, 100))),
  WT    = c(6, 6, 6, 6),
  CRCL  = c(20, 50, 70, 100)
)

# Helper: build an event table for one cohort with n simulated subjects,
# dosing 0.1 mg/kg morphine sulfate every 3 h for 24 h (paper's simulation
# regimen). Doses are converted from mg morphine sulfate to nmol of free
# morphine via the paper's molar-conversion: 1 mg sulfate -> 2 * 1e6 / 758.83
# nmol of free morphine (one sulfate molecule contains two morphine moieties).
mg_sulfate_to_nmol_morphine <- function(mg) mg * 2 * 1e6 / 758.83

make_cohort <- function(group, WT, CRCL, n, id_offset = 0L) {
  dose_amt_nmol  <- mg_sulfate_to_nmol_morphine(0.1 * WT)  # 0.1 mg/kg morphine sulfate
  dose_times     <- seq(0, 21, by = 3)                       # Q3h for 24 h => 8 doses
  obs_times      <- sort(unique(c(seq(0, 24, by = 0.25))))  # 15-min observation grid

  per_id <- function(id) {
    dose_rows <- data.frame(
      id   = id,
      time = dose_times,
      amt  = dose_amt_nmol,
      evid = 1L,
      cmt  = "central",
      WT   = WT,
      CRCL = CRCL,
      group = as.character(group),
      stringsAsFactors = FALSE
    )
    obs_rows <- data.frame(
      id   = id,
      time = obs_times,
      amt  = NA_real_,
      evid = 0L,
      cmt  = "Cc",
      WT   = WT,
      CRCL = CRCL,
      group = as.character(group),
      stringsAsFactors = FALSE
    )
    dplyr::bind_rows(dose_rows, obs_rows)
  }

  dplyr::bind_rows(lapply(id_offset + seq_len(n), per_id)) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

n_per_group <- 100L

events_morphine <- dplyr::bind_rows(
  make_cohort(morphine_groups$group[1], morphine_groups$WT[1], morphine_groups$CRCL[1],
              n_per_group, id_offset =   0L),
  make_cohort(morphine_groups$group[2], morphine_groups$WT[2], morphine_groups$CRCL[2],
              n_per_group, id_offset = 100L),
  make_cohort(morphine_groups$group[3], morphine_groups$WT[3], morphine_groups$CRCL[3],
              n_per_group, id_offset = 200L)
)

events_gfr <- dplyr::bind_rows(
  make_cohort(gfr_groups$group[1], gfr_groups$WT[1], gfr_groups$CRCL[1],
              n_per_group, id_offset =   0L),
  make_cohort(gfr_groups$group[2], gfr_groups$WT[2], gfr_groups$CRCL[2],
              n_per_group, id_offset = 100L),
  make_cohort(gfr_groups$group[3], gfr_groups$WT[3], gfr_groups$CRCL[3],
              n_per_group, id_offset = 200L),
  make_cohort(gfr_groups$group[4], gfr_groups$WT[4], gfr_groups$CRCL[4],
              n_per_group, id_offset = 300L)
)

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

# Simulation

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

sim_morphine <- rxode2::rxSolve(
  mod, events = events_morphine, keep = c("group", "WT", "CRCL")
) |> as.data.frame()

sim_gfr <- rxode2::rxSolve(
  mod, events = events_gfr, keep = c("group", "WT", "CRCL")
) |> as.data.frame()
```

# Replicate published figures

## Figure 7 (upper row): morphine AUC0-24h across age groups

```{r figure-7-morphine}
trapz <- function(x, y) {
  ord <- order(x); x <- x[ord]; y <- y[ord]
  sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2)
}

auc_morphine <- sim_morphine |>
  dplyr::group_by(id, group) |>
  dplyr::summarise(
    AUC_nM_h = trapz(time, Cc),
    .groups  = "drop"
  ) |>
  # Convert nM*h to ng*h/mL (paper unit): 1 nM = 1 nmol/L = (285.34 ng/mol) ng/L
  # => 1 nM*h = 285.34 / 1000 ng/mL*h = 0.28534 ng*h/mL.
  dplyr::mutate(AUC_ng_h_mL = AUC_nM_h * 285.34 / 1000)

ggplot(auc_morphine, aes(group, AUC_ng_h_mL)) +
  geom_boxplot(fill = "#a3c4f3") +
  labs(x = NULL,
       y = "Morphine AUC0-24h (ng*h/mL)",
       title = "Replicates Figure 7 (upper row) of Elkomy 2015",
       caption = paste(
         "Morphine exposure across three age groups under 0.1 mg/kg morphine",
         "sulfate Q3h x 24 h. Median AUC0-24h (paper Results, final paragraph):",
         "331 (1 mo), 353 (6 mo), 408 (2 yr) ng*h/mL."))
```

```{r morphine-medians}
morphine_medians <- auc_morphine |>
  dplyr::group_by(group) |>
  dplyr::summarise(`median AUC0-24h (ng*h/mL)` = median(AUC_ng_h_mL),
                   .groups = "drop")
knitr::kable(morphine_medians,
             caption = "Simulated median morphine AUC0-24h by age group. Compare to paper-reported 331, 353, 408 ng*h/mL (Results).")
```

## Figure 7 (lower row): M3G and M6G AUC0-24h across GFR groups

```{r figure-7-metabolites}
auc_metab <- sim_gfr |>
  dplyr::group_by(id, group) |>
  dplyr::summarise(
    AUC_M3G_nM_h = trapz(time, Cc_m3g),
    AUC_M6G_nM_h = trapz(time, Cc_m6g),
    .groups = "drop"
  ) |>
  # Convert nM*h to ng*h/mL using metabolite MW 461.46.
  dplyr::mutate(
    AUC_M3G_ng_h_mL = AUC_M3G_nM_h * 461.46 / 1000,
    AUC_M6G_ng_h_mL = AUC_M6G_nM_h * 461.46 / 1000
  ) |>
  tidyr::pivot_longer(
    cols = c(AUC_M3G_ng_h_mL, AUC_M6G_ng_h_mL),
    names_to = "metabolite", values_to = "AUC"
  ) |>
  dplyr::mutate(metabolite = ifelse(metabolite == "AUC_M3G_ng_h_mL", "M3G", "M6G"))

ggplot(auc_metab, aes(group, AUC, fill = metabolite)) +
  geom_boxplot() +
  facet_wrap(~metabolite, scales = "free_y") +
  labs(x = NULL,
       y = "AUC0-24h (ng*h/mL)",
       title = "Replicates Figure 7 (lower row) of Elkomy 2015",
       caption = paste(
         "M3G and M6G exposure in a 6-month-old (6 kg) child across four",
         "renal-function strata (GFR in mL/min/1.73 m^2) under 0.1 mg/kg",
         "morphine sulfate Q3h x 24 h. Paper Discussion: GFR=20 vs GFR=70",
         "metabolite AUC is ~2.7-fold higher."))
```

```{r gfr-ratios}
gfr_medians <- auc_metab |>
  dplyr::group_by(metabolite, group) |>
  dplyr::summarise(`median AUC0-24h (ng*h/mL)` = median(AUC),
                   .groups = "drop") |>
  tidyr::pivot_wider(names_from = metabolite,
                     values_from = `median AUC0-24h (ng*h/mL)`)
knitr::kable(gfr_medians,
             caption = "Simulated median metabolite AUC0-24h by GFR group (6-month-old, 6 kg).")

# Severe (GFR=20) vs typical (GFR=70) ratio per metabolite -- paper reports
# ~2.7-fold.
sev_vs_norm <- auc_metab |>
  dplyr::group_by(metabolite, group) |>
  dplyr::summarise(median_AUC = median(AUC), .groups = "drop") |>
  tidyr::pivot_wider(names_from = group, values_from = median_AUC)
sev_vs_norm$ratio_GFR20_vs_GFR70 <- sev_vs_norm$`GFR 20` / sev_vs_norm$`GFR 70`
knitr::kable(sev_vs_norm,
             caption = "GFR=20 vs GFR=70 median AUC ratio per metabolite. Paper reports ~2.7-fold (Discussion).")
```

# PKNCA validation

PKNCA-derived single-dose NCA after a single 0.1 mg/kg IV bolus of morphine
sulfate, computed separately for morphine and each glucuronide metabolite.
Treatment grouping is by `group` (the three age cohorts at typical GFR = 70).

```{r pknca-setup}
# Reuse the morphine cohort, but with a single dose at t = 0 and sampling
# from 0 to 24 h.
make_single_dose_cohort <- function(group, WT, CRCL, n, id_offset = 0L) {
  dose_amt_nmol  <- mg_sulfate_to_nmol_morphine(0.1 * WT)
  obs_times      <- sort(unique(c(seq(0, 24, by = 0.25))))

  per_id <- function(id) {
    dose_rows <- data.frame(
      id   = id,
      time = 0,
      amt  = dose_amt_nmol,
      evid = 1L,
      cmt  = "central",
      WT   = WT,
      CRCL = CRCL,
      group = as.character(group),
      stringsAsFactors = FALSE
    )
    obs_rows <- data.frame(
      id   = id,
      time = obs_times,
      amt  = NA_real_,
      evid = 0L,
      cmt  = "Cc",
      WT   = WT,
      CRCL = CRCL,
      group = as.character(group),
      stringsAsFactors = FALSE
    )
    dplyr::bind_rows(dose_rows, obs_rows)
  }

  dplyr::bind_rows(lapply(id_offset + seq_len(n), per_id)) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events_sd <- dplyr::bind_rows(
  make_single_dose_cohort(morphine_groups$group[1], morphine_groups$WT[1], morphine_groups$CRCL[1],
                          n_per_group, id_offset =   0L),
  make_single_dose_cohort(morphine_groups$group[2], morphine_groups$WT[2], morphine_groups$CRCL[2],
                          n_per_group, id_offset = 100L),
  make_single_dose_cohort(morphine_groups$group[3], morphine_groups$WT[3], morphine_groups$CRCL[3],
                          n_per_group, id_offset = 200L)
)
stopifnot(!anyDuplicated(unique(events_sd[, c("id", "time", "evid")])))

sim_sd <- rxode2::rxSolve(mod, events = events_sd,
                          keep = c("group", "WT", "CRCL")) |>
  as.data.frame()

dose_df <- events_sd |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, group)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | group + id,
                             doseu = "nmol")

intervals <- data.frame(
  start      = 0,
  end        = 24,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  aucinf.obs = TRUE,
  half.life  = TRUE
)
```

## Morphine

```{r pknca-morphine}
sim_nca_morphine <- sim_sd |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::select(id, time, Cc, group)

conc_obj_morphine <- PKNCA::PKNCAconc(
  sim_nca_morphine, Cc ~ time | group + id,
  concu = "nmol/L", timeu = "hour"
)

nca_morphine <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_morphine, dose_obj, intervals = intervals)
)
knitr::kable(as.data.frame(summary(nca_morphine)),
             caption = "Morphine NCA after single 0.1 mg/kg morphine sulfate IV bolus.")
```

## M3G

```{r pknca-m3g}
sim_nca_m3g <- sim_sd |>
  dplyr::filter(!is.na(Cc_m3g), time > 0) |>
  dplyr::select(id, time, Cc_m3g, group) |>
  dplyr::rename(Cc = Cc_m3g)

conc_obj_m3g <- PKNCA::PKNCAconc(
  sim_nca_m3g, Cc ~ time | group + id,
  concu = "nmol/L", timeu = "hour"
)

nca_m3g <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_m3g, dose_obj, intervals = intervals)
)
knitr::kable(as.data.frame(summary(nca_m3g)),
             caption = "M3G NCA after single 0.1 mg/kg morphine sulfate IV bolus.")
```

## M6G

```{r pknca-m6g}
sim_nca_m6g <- sim_sd |>
  dplyr::filter(!is.na(Cc_m6g), time > 0) |>
  dplyr::select(id, time, Cc_m6g, group) |>
  dplyr::rename(Cc = Cc_m6g)

conc_obj_m6g <- PKNCA::PKNCAconc(
  sim_nca_m6g, Cc ~ time | group + id,
  concu = "nmol/L", timeu = "hour"
)

nca_m6g <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_m6g, dose_obj, intervals = intervals)
)
knitr::kable(as.data.frame(summary(nca_m6g)),
             caption = "M6G NCA after single 0.1 mg/kg morphine sulfate IV bolus.")
```

## Comparison against published values

The paper reports the median morphine AUC0-24h at the multiple-dose
0.1 mg/kg-Q3h regimen as 331 ng*h/mL (1 mo, 4.5 kg), 353 ng*h/mL (6 mo, 6 kg),
and 408 ng*h/mL (2 yr, 12 kg) (Results, final paragraph). The simulated
medians from the multiple-dose simulation above are shown in the
`morphine_medians` table; the same monotonically-increasing pattern with age
emerges, with the 1-month and 6-month values close to one another and the
2-year-old value substantially higher, matching the paper's pairwise
comparison (p > 0.05 between 1 mo and 6 mo, p < 0.05 between 6 mo and 2 yr).

For metabolite exposure under reduced renal function, the paper reports an
approximately 2.7-fold higher metabolite AUC at GFR = 20 versus GFR = 70 in
a 6-month-old (Results, final paragraph). The `sev_vs_norm` table above shows
the simulated GFR=20 vs GFR=70 ratio for M3G and M6G; the value is identically
determined by the model's GFR^(+1) on Kint and GFR^(-1) on Mmax under
typical-value scaling, giving an asymptotic GFR=20/70 = 3.5-fold ratio at
near-saturation (`M_max * Cint / (Cint50 + Cint)` saturates well above
Cint50 for the paper's morphine dosing range) and approaching the typical
linear `(70/20)` = 3.5 ratio in the deep-tail at lower exposures; the paper's
2.7-fold value reflects the AUC integrated over the Q3h dosing-period and
includes between-subject variability not present in the deterministic
asymptote.

# Assumptions and deviations

- **Dosing unit conversion.** The packaged model accepts dose in `nmol` of
  free morphine equivalents, matching the molar units the paper used for
  parameter estimation (Methods para 1). Clinical doses are typically
  reported in mg of morphine sulfate; the vignette includes a helper
  `mg_sulfate_to_nmol_morphine()` that applies the paper's conversion
  (1 mg sulfate -> 2 * 1e6 / 758.83 nmol of free morphine) before passing
  to `rxSolve`. The same conversion would apply in any downstream
  user-facing simulation.
- **Reference body weight 6 kg.** The Methods para 4 statement
  "allometric body weight-based model scaled to a 6-kg child (median body
  weight in this study)" identifies 6 kg as the explicit allometric
  reference. The simulated 6-month cohort therefore matches the reference
  exactly; the 1-month (4.5 kg) and 2-year (12 kg) cohorts inherit the
  paper's allometric scaling.
- **Reference GFR 70 mL/min/1.73 m^2.** The Schwartz-formula-derived eGFR
  enters the Kint and Mmax covariate equations normalised to 70 mL/min
  /1.73 m^2 (paper Table II model column). The reference value is read
  directly from the table; the population mean GFR (79) is close but not
  identical to the reference.
- **Allometric exponents 0.75 / 1.0 FIXED.** Paper Methods para 4 and
  Results para 2: the authors tested free-exponent estimation and the OFV
  reduction was below the 3.84 threshold (p > 0.05, 1 df); they therefore
  fixed the canonical Anderson-Holford values. The packaged model wraps
  the exponents in `fixed()` in `ini()`.
- **GFR exponents +/-1 FIXED.** Paper Results, final paragraph: the
  authors fit free GFR exponents on Kint and Mmax, found them very close
  to unity, and fixed them to +1 and -1 with only a +2 OFV penalty. The
  packaged model inlines the +/-1 exponents in `model()` rather than
  carrying explicit `e_crcl_*` parameters, since the structural assumption
  is a paper-imposed equality rather than a free fitted parameter.
- **Shared M3G / M6G inter-individual variability.** Paper Results
  ("Between-subject random effects of M3G and M6G were very close.
  Therefore, inter-subject variability parameters were shared between M3G
  and M6G."): the packaged model encodes a single eta per metabolite-
  channel parameter group (`etalkeo`, `etalemax`, `etalec50`) applied to
  both `lkeo_m3g` / `lkeo_m6g`, `lemax_m3g` / `lemax_m6g`, and
  `lec50_m3g` / `lec50_m6g`. Per Table II this matches the paper's
  Final-model column (one IIV row per group rather than per metabolite).
- **Log-additive residual error -> proportional on linear scale.** Paper
  Methods para 1: concentrations were log-transformed and an additive
  residual error model was used. This is equivalent to a proportional
  error model on the linear concentration scale with `propSd ~=
  sqrt(sigma^2)` (verification-checklist.md A). The packaged model
  declares `Cc ~ prop(propSd)`, `Cc_m3g ~ prop(propSd_m3g)`, and
  `Cc_m6g ~ prop(propSd_m6g)`.
- **Effect-compartment ODE form for the "intermediate" compartment.**
  Paper Eq. 1 reads `dCint/dt = Kint * (C - Cint)`, identical to the
  Sheiner / Holford effect-compartment ODE form (canonical compartment
  name `effect` in nlmixr2lib conventions). The packaged model uses
  `d/dt(effect_m3g) <- keo_m3g * (Cc - effect_m3g)` and similarly for
  M6G, with the metabolite-suffix `_m3g` / `_m6g` per the
  registered-metabolite naming convention.
- **Algebraic Emax transduction.** Paper Eq. 2 reads
  `M = Mmax * Cint / (Cint50 + Cint)`; the packaged model declares
  `Cc_m3g` and `Cc_m6g` as algebraic outputs (not ODE states) computed
  inside `model()` from the corresponding effect-compartment state.
  Metabolite concentrations therefore have no mass-balance or
  elimination dynamics of their own -- they are deterministic Emax
  transductions of the intermediate signal, as in the source paper.
- **Initial condition for the effect compartment.** The model assumes
  `effect_m3g(0) = effect_m6g(0) = 0`; before any morphine is
  administered, neither metabolite signal exists. The paper does not
  explicitly state the initial condition but Eq. 1 with `C = 0` gives
  `Cint -> 0` as the stationary point.
- **Population baseline gaps.** The paper does not report race or
  ethnicity; `population$race_ethnicity` is left `NULL`. Postnatal age
  was not a significant covariate once weight was included; the
  vignette does not stratify by PNA beyond the three age-by-weight
  groups Figure 7 uses.
- **Cohort sizes.** The paper used 100 simulated subjects per group for
  Figure 7; the vignette matches that. The variance reported in the
  vignette's PKNCA tables therefore reflects the same Monte-Carlo
  sample size, but the random-number stream differs.
