---
title: "Nefopam (Djerada 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Nefopam (Djerada 2014)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Djerada_2014_nefopam"))$reference`
* Article: <https://doi.org/10.1111/bcp.12291>

The packaged model is a joint parent + metabolite population PK model
for nefopam and its N-demethyl metabolite desmethyl-nefopam. Nefopam
is fit as a 2-compartment structural model with an IV zero-order
infusion delivered directly into the central compartment (no depot).
Desmethyl-nefopam is fit as a 1-compartment structural model that
receives an apparent first-order metabolic flux `cl_form_dnef *
central / vc` from the nefopam central compartment and eliminates
with total clearance `cl_dnef` (paper symbol K13; nefopam central to
desmethyl-nefopam formation clearance -- Figure 1B and Results
p 1031). The "apparent" qualifier on K13 and V_dnef absorbs the
fraction metabolised and the metabolite bioavailability into single
positive coefficients; the paper does not decompose these separately.
No covariates are retained. Between-subject variability is described
by exponential IIV on all seven structural parameters with a diagonal
Omega; residual variability is a combined additive + proportional
error model on each of the nefopam and desmethyl-nefopam
concentration outputs.

# Population

The study enrolled 48 elderly patients (65-99 years; mean 83 +/- 8)
scheduled for postoperative repair of a fractured hip at Reims
University Hospital (France) between June 2006 and August 2011.
Baseline demographics (Djerada 2014 Table 1): global weight
63 +/- 13 kg (range 40-100), height 161 +/- 7 cm (144-183),
BMI 23.9 +/- 4.7 kg/m^2, sex ratio male:female = 0.45 (15 males,
33 females = 68.75% female). Renal function was classified a priori
using the Cockroft-Gault GFR into three strata: normal (GFR > 60
mL/min per 1.73 m^2; n = 10), moderate impairment (GFR < 60; n = 28),
and severe impairment (GFR <= 30; n = 10). Iohexol-clearance GFR
(GFR_IO, a gold-standard direct measure) global mean 48.36 +/- 24.19
mL/min (range 10.77-131.50). Exclusion criteria included nefopam or
iodinated-contrast hypersensitivity, hepatic insufficiency, epilepsy,
glaucoma, prostate adenoma, congestive heart failure, unstable
angina, cardiac-infarction sequelae, uncontrolled arrhythmia, and
haemodialysis.

Dosing: a single 20 mg nefopam hydrochloride IV infusion (diluted in
0.9% saline) over 30 min, delivered postoperatively. Iohexol 5 mL of
Omnipaque 180 was administered as a 1 min IV infusion starting 29 min
after nefopam initiation to estimate GFR. Plasma nefopam,
desmethyl-nefopam, and iohexol concentrations were sampled at
0, 5, 10, 15, 30, 60, 120, 180, 240, 480, and 1440 min after nefopam
administration.

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

# Source trace

Every `ini()` parameter carries a source-location comment in the
model file (`inst/modeldb/specificDrugs/Djerada_2014_nefopam.R`); the
table below collects the audit trail in one place. All values come
from the "Nefopam pharmacokinetic model" and "Desmethyl-nefopam
pharmacokinetic model" paragraphs on p 1031 of Djerada 2014. BSVs are
reported as CV percent; the model file converts to variance via
`omega^2 = log(CV^2 + 1)` with CV expressed as a fraction.

| Equation / parameter    | Value (typical)     | BSV (CV %) | Source location |
|-------------------------|---------------------|------------|-----------------|
| `lcl`                   | log(17.3) (L/h)     | 53         | Results p 1031, "CL_nef = 17.3 +/- 1.9 L/h" |
| `lvc`                   | log(114) (L)        | 121        | Results p 1031, "V_1nef = 114 +/- 21 L" |
| `lq`                    | log(80.7) (L/h)     | 79         | Results p 1031, "Q_nef = 80.7 +/- 15.0 L/h" |
| `lvp`                   | log(208) (L)        | 64         | Results p 1031, "V_2nef = 208 +/- 31 L" |
| `lcl_dnef`              | log(17.81) (L/h)    | 144        | Results p 1031, "CL_dnef = 17.81 +/- 0.6 L/h" |
| `lvc_dnef`              | log(12.5) (L)       | 1          | Results p 1031, "V_dnef = 12.5 +/- 0.01 L" |
| `lcl_form_dnef` (= K13) | log(0.35) (L/h)     | 120        | Figure 1B legend + Results p 1031, "K13 = 0.35 +/- 0.022 L/h" |
| `addSd` (nefopam)       | 1.33 ug/L           | -          | Results p 1031, "additive (a = 1.33 ug/L)" |
| `propSd` (nefopam)      | 0.228 (fraction)    | -          | Results p 1031, "proportional (coefficient = 0.228)" |
| `addSd_dnef`            | 0.01 ug/L           | -          | Results p 1031, "additive (a = 0.01 ug/L)" (dnef error model) |
| `propSd_dnef`           | 0.68 (fraction)     | -          | Results p 1031, "proportional (coefficient = 0.68)" (dnef error model) |
| Omega structure         | diagonal (7 IIVs)   | -          | Methods p 1030, "parameterized as a diagonal matrix. Correlations between random effects were tested" (none retained in the final model) |
| `d/dt(central)`         | q * peripheral1/vp - (cl + q) * central/vc | - | Figure 1A (nefopam 2-cmt with zero-order input); apparent CL_nef captures all elimination |
| `d/dt(peripheral1)`     | q * central/vc - q * peripheral1/vp | - | Figure 1A |
| `d/dt(central_dnef)`    | cl_form_dnef * central/vc - cl_dnef * central_dnef/vc_dnef | - | Figure 1B (dnef 1-cmt with first-order metabolic input from nefopam) |
| Covariate effects       | none retained       | -          | Results p 1031, "None of the covariates met the predefined inclusion criteria" (age, sex, weight, height, BMI, LBW, FFM, IBW, BSA, GFR_IO, MDRD, Cockroft-Gault all screened) |

# Virtual cohort

The original observed data are not publicly available. The
simulations below use a virtual cohort of 200 elderly subjects whose
covariate distributions approximate the published Table 1
demographics. Because the final population model has no retained
covariates, the covariate values do not enter the simulation; they
are carried through only for documentation and for the post-hoc
severe-vs-non-severe stratification of NCA statistics.

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

n_subj <- 200L
mod <- readModelDb("Djerada_2014_nefopam")

# Table 1 global demographics: age 83 +/- 8, weight 63 +/- 13,
# height 161 +/- 7, BMI 23.9 +/- 4.7, 68.75% female, GFR_IO
# 48.36 +/- 24.19 (truncated at 0 to avoid negative values).
cohort <- tibble(
  id      = seq_len(n_subj),
  AGE     = pmin(pmax(rnorm(n_subj, mean = 83, sd = 8), 65), 99),
  SEXF    = as.integer(runif(n_subj) < 0.6875),
  WT      = pmin(pmax(rnorm(n_subj, mean = 63, sd = 13), 40), 100),
  HT      = pmin(pmax(rnorm(n_subj, mean = 161, sd = 7), 144), 183),
  GFR_IO  = pmax(rnorm(n_subj, mean = 48.36, sd = 24.19), 5)
) |>
  mutate(BMI      = WT / (HT / 100)^2,
         renal_gp = case_when(GFR_IO >  60 ~ "normal",
                              GFR_IO <= 30 ~ "severe",
                              TRUE         ~ "moderate"))

# Observation grid matches the paper's sampling schedule
# (converted from minutes to hours since the model's time unit is
# hours): 0, 5, 10, 15, 30, 60, 120, 180, 240, 480, 1440 min.
obs_times_h <- c(0, 5, 10, 15, 30, 60, 120, 180, 240, 480, 1440) / 60

# Denser grid for the Figure 5 replication (fine profile to catch
# Cmax and its rate of increase RCmax accurately).
dense_times_h <- sort(unique(c(obs_times_h,
                               seq(0, 1.5, by = 0.02),
                               seq(1.5, 24, by = 0.25))))

# Event helper: 20 mg IV infusion, delivered as a zero-order input on
# `central` over `infusion_h` hours. rxode2's `dur` column carries the
# infusion duration; `amt` is the total dose in mg. `treatment`
# labels the infusion regimen so it can ride through rxSolve via keep.
build_events <- function(cohort, infusion_h, obs_grid, treatment_label) {
  dose_rows <- cohort |>
    transmute(id,
              time = 0,
              amt  = 20,
              dur  = infusion_h,
              evid = 1L,
              cmt  = "central",
              treatment = treatment_label,
              AGE, SEXF, WT, HT, BMI, GFR_IO, renal_gp)
  obs_rows <- cohort |>
    select(id, AGE, SEXF, WT, HT, BMI, GFR_IO, renal_gp) |>
    tidyr::crossing(time = obs_grid) |>
    mutate(amt = NA_real_, dur = NA_real_,
           evid = 0L, cmt = "Cc",
           treatment = treatment_label)
  bind_rows(dose_rows, obs_rows) |>
    arrange(id, time, desc(evid)) |>
    as.data.frame()
}

events_30min <- build_events(cohort, infusion_h = 0.5,
                             obs_grid = dense_times_h,
                             treatment_label = "30 min infusion")
```

# Simulation

The typical dosing regimen in Djerada 2014 is 20 mg of nefopam over
30 min. Simulate the full cohort at this regimen with between-subject
variability enabled; carry the covariate columns through with `keep`
so downstream figures / NCA can stratify by renal-function group.

```{r simulate}
sim_30min <- rxode2::rxSolve(
  mod,
  events = events_30min,
  keep   = c("treatment", "renal_gp", "GFR_IO", "AGE", "WT", "SEXF")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2 C, D: prediction-corrected VPC for nefopam and desmethyl-nefopam

The paper's Figure 2 C / D shows the prediction-corrected VPC across
1000 simulated data sets, with the observed 5th / 50th / 95th
percentiles compared against the 90% confidence intervals of the
simulated percentiles. The replication below plots the simulated 5th
/ 50th / 95th percentiles from the virtual cohort at each nominal
sampling time; the shape should be the same as Figure 2C / D even
though the observations are not available to overlay.

```{r figure-2cd, fig.height = 5}
vpc_data <- sim_30min |>
  dplyr::filter(time %in% obs_times_h, time > 0) |>
  dplyr::select(time, Cc, Cc_dnef) |>
  tidyr::pivot_longer(c(Cc, Cc_dnef),
                      names_to = "analyte", values_to = "conc") |>
  dplyr::mutate(analyte = dplyr::recode(analyte,
                                        Cc = "Nefopam",
                                        Cc_dnef = "Desmethyl-nefopam")) |>
  dplyr::group_by(analyte, time) |>
  dplyr::summarise(
    Q05 = quantile(conc, 0.05, na.rm = TRUE),
    Q50 = quantile(conc, 0.50, na.rm = TRUE),
    Q95 = quantile(conc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_data, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~analyte, scales = "free_y", ncol = 2) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Plasma concentration (ug/L)",
       title = "Figure 2 C / D -- prediction-corrected VPC",
       caption = paste("Replicates Figure 2 C (nefopam) and 2 D",
                       "(desmethyl-nefopam) of Djerada 2014.",
                       "Ribbon: simulated 5th-95th percentile;",
                       "line: simulated median.")) +
  theme_bw()
```

## Figure 5: probability of Cmax / RCmax / AUC crossing cut-values vs infusion duration

Djerada 2014 Figure 5 shows how the probabilities of `Cmax_nef` being
below its tachycardia / PONV cut-value (389 ug/L), `RCmax_nef` being
below its cut-value (764 ug/L/h), and `AUC0-infinity_nef` being above
the morphine-sparing cut-value (950 ug*h/L) evolve as a function of
the infusion duration. Reproduce this by simulating the same cohort
at four different infusion durations (15, 30, 45, and 60 min) and
tallying the per-subject exceedances of each cut-value.

```{r figure-5, fig.height = 4.4}
durations_min <- c(15, 30, 45, 60)

infusion_sims <- lapply(durations_min, function(dur_min) {
  ev <- build_events(cohort, infusion_h = dur_min / 60,
                     obs_grid = dense_times_h,
                     treatment_label = paste0(dur_min, " min infusion"))
  rxode2::rxSolve(mod, events = ev,
                  keep = c("treatment", "renal_gp")) |>
    as.data.frame() |>
    dplyr::mutate(dur_min = dur_min)
})

# For each subject and each infusion duration, compute:
#   * Cmax_nef  = max Cc over 0-24 h
#   * RCmax_nef = Cmax / t_Cmax                 (paper definition: rate of increase
#                                                from t=0 to t=Cmax)
#   * AUC_nef   = trapezoidal AUC over 0-24 h
per_subject <- dplyr::bind_rows(infusion_sims) |>
  dplyr::filter(time <= 24) |>
  dplyr::group_by(dur_min, id) |>
  dplyr::summarise(
    Cmax_nef  = max(Cc, na.rm = TRUE),
    tmax_nef  = time[which.max(Cc)],
    AUC_nef   = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2,
                    na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(RCmax_nef = ifelse(tmax_nef > 0, Cmax_nef / tmax_nef, NA_real_))

probs <- per_subject |>
  dplyr::group_by(dur_min) |>
  dplyr::summarise(
    `P(Cmax < 389 ug/L)`         = mean(Cmax_nef  < 389),
    `P(RCmax < 764 ug/L/h)`      = mean(RCmax_nef < 764),
    `P(AUC > 950 ug*h/L)`        = mean(AUC_nef   > 950),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(-dur_min, names_to = "criterion", values_to = "prob")

ggplot(probs, aes(dur_min, prob, colour = criterion, shape = criterion)) +
  geom_line() + geom_point(size = 2) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(breaks = durations_min) +
  labs(x = "Infusion duration (min)", y = "Probability",
       colour = NULL, shape = NULL,
       title = "Figure 5 -- effect of infusion duration on exposure metrics",
       caption = paste("Replicates Figure 5 of Djerada 2014. Cut-values",
                       "from Table 3 (tachycardia PMxax cut 389 ug/L,",
                       "RCmax cut 764 ug/L/h) and Results p 1032",
                       "(morphine-sparing AUC cut 950 ug*h/L).")) +
  theme_bw() + theme(legend.position = "bottom")

knitr::kable(
  probs |> tidyr::pivot_wider(names_from = criterion, values_from = prob),
  digits = 3,
  caption = paste("Simulated probabilities of the tachycardia,",
                  "PONV, and morphine-sparing cut-values as a function",
                  "of infusion duration. Djerada 2014 Figure 5 reports",
                  "P(Cmax<389) and P(RCmax<764) both ~0.94 at 45 min",
                  "and ~0.999 at 60 min; P(AUC>950) stays ~0.72 across",
                  "all durations (approximately 0.28 fall below the",
                  "morphine-sparing cut-value regardless of duration)."))
```

# PKNCA validation -- nefopam and desmethyl-nefopam

Compute noncompartmental parameters for the 30 min infusion regimen
using PKNCA, stratifying by simulated renal-function group. Nefopam
and desmethyl-nefopam are handled in separate PKNCA blocks
(multi-output convention).

```{r pknca-nef}
# Nefopam concentrations
sim_nca_nef <- sim_30min |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, renal_gp)

# Guarantee a time=0 row per (id, renal_gp).
sim_nca_nef <- dplyr::bind_rows(
  sim_nca_nef,
  sim_nca_nef |> dplyr::distinct(id, renal_gp) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, renal_gp, time, .keep_all = TRUE) |>
  dplyr::arrange(id, renal_gp, time)

conc_nef <- PKNCA::PKNCAconc(sim_nca_nef, Cc ~ time | renal_gp + id)

dose_nef <- events_30min |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt, dur, renal_gp)
dose_obj_nef <- PKNCA::PKNCAdose(
  dose_nef, amt ~ time | renal_gp + id,
  duration = "dur", route = "intravascular")

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

nca_nef <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_nef, dose_obj_nef, intervals = intervals_nef)
)
```

```{r pknca-dnef}
# Desmethyl-nefopam concentrations (no exogenous dose; formation
# clearance from nefopam is the only input).
sim_nca_dnef <- sim_30min |>
  dplyr::filter(!is.na(Cc_dnef)) |>
  dplyr::select(id, time, Cc_dnef, renal_gp)

sim_nca_dnef <- dplyr::bind_rows(
  sim_nca_dnef,
  sim_nca_dnef |> dplyr::distinct(id, renal_gp) |>
    dplyr::mutate(time = 0, Cc_dnef = 0)
) |>
  dplyr::distinct(id, renal_gp, time, .keep_all = TRUE) |>
  dplyr::arrange(id, renal_gp, time)

conc_dnef <- PKNCA::PKNCAconc(sim_nca_dnef, Cc_dnef ~ time | renal_gp + id)

# Nefopam dose drives desmethyl-nefopam kinetics; register the same
# nefopam dose here so PKNCA has a dosing anchor even though the
# metabolite has no independent exogenous input.
dose_obj_dnef <- PKNCA::PKNCAdose(
  dose_nef, amt ~ time | renal_gp + id,
  duration = "dur", route = "intravascular")

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

nca_dnef <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_dnef, dose_obj_dnef, intervals = intervals_dnef)
)
```

## Comparison against published NCA

Djerada 2014 Table 2 reports individual PK parameters stratified by
severe (GFR_IO <= 30 mL/min; n = 10) vs non-severe (n = 38) renal
impairment. Note that the population model has NO retained
renal-function covariate; the Table 2 stratification is a post-hoc
split of EBE-derived individual parameters. Because the simulation
uses the same model equations for every subject regardless of GFR,
the simulated distributions in each renal-function group will
NOT reproduce the published stratified values -- the simulated
severe-vs-non-severe split is expected to be much smaller than the
published one, and the comparison flags this cleanly as an
above-tolerance discrepancy on `cmax`, `auc`, and half-life.

The reference below is a weighted-average "global" estimate derived
from Djerada 2014 Table 2 (n = 10 severe + n = 38 non-severe = 48
patients):

* Global mean CL_nef      = (10 x 13.04 + 38 x 24.61) / 48 = 22.20 L/h
* Global mean Cmax_nef    = (10 x 156.2 + 38 x 148.2) / 48 = 149.87 ug/L
* Global mean AUC0inf_nef = (10 x 1617.0 + 38 x 1184.0) / 48 = 1274.2 ug*h/L
* Global mean Cmax_dnef   = (10 x 7.0  + 38 x  4.9) / 48 =   5.34 ug/L
* Global mean AUC0inf_dnef= (10 x 141.7 + 38 x 95.1) / 48 = 104.8 ug*h/L

These global means are shown in a single "global" row of the
reference table below. Simulated per-renal-group values are shown
alongside; because the model has no renal-function effect, the
severe and non-severe simulated rows should agree with each other
(and with the global reference) rather than reproduce the
published Table 2 stratification.

```{r nca_compare_nef}
published_nef <- tibble::tribble(
  ~cmax,   ~tmax,   ~aucinf.obs,
  149.87,  0.5,     1274.2
)

cmp_nef <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_nef,
  reference     = published_nef,
  by            = NULL,
  units         = c(cmax = "ug/L", aucinf.obs = "ug*h/L",
                    tmax = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_nef,
  digits  = 2,
  caption = paste("Nefopam: simulated pooled median across the",
                  "virtual cohort vs. published Djerada 2014 Table 2",
                  "weighted-global means (Cmax = 149.9 ug/L,",
                  "Tmax = 0.5 h at end of 30 min infusion,",
                  "AUC0-inf = 1274 ug*h/L).",
                  "* differs from reference by >20%."))
```

```{r nca_compare_dnef}
published_dnef <- tibble::tribble(
  ~cmax,  ~tmax,   ~aucinf.obs,
  5.34,   3.0,     104.8
)

cmp_dnef <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_dnef,
  reference     = published_dnef,
  by            = NULL,
  units         = c(cmax = "ug/L", aucinf.obs = "ug*h/L",
                    tmax = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_dnef,
  digits  = 2,
  caption = paste("Desmethyl-nefopam: simulated pooled median across",
                  "the virtual cohort vs. published Djerada 2014",
                  "Table 2 weighted-global means (Cmax = 5.3 ug/L,",
                  "AUC0-inf = 105 ug*h/L).",
                  "* differs from reference by >20%."))
```

# Assumptions and deviations

- **Concentration units and the 1000 x scaling in `model()`.** Doses
  are entered in milligrams (matching the paper's 20 mg dose
  language) and volumes are in litres, so the intrinsic `central /
  vc` ratio is in mg / L. The paper reports concentrations in
  ug / L (equivalent to ng / mL). The model applies a `1000 *`
  factor when computing `Cc` and `Cc_dnef` to convert mg/L to ug/L
  so simulated concentrations agree with the paper's numerical
  scale.
- **No covariate stratification in the simulation.** The paper's
  final model has no retained covariates (Results p 1031); the
  Table 2 severe / non-severe stratification is a post-hoc split of
  individual EBE parameters, not a covariate effect. The virtual
  cohort therefore samples all patients from the same
  covariate-independent distribution regardless of GFR / age / etc.
  Assigning a `renal_gp` label from the simulated `GFR_IO` value is
  documentation-only. The `NCA comparison` section flags the
  expected disagreement between simulated per-group values and
  published Table 2 stratified values.
- **`V_dnef` IIV of 1% is preserved at face value.** Djerada 2014
  reports BSV of `V_dnef` as 1 +/- 0.9 % (Results p 1031), a
  near-zero value that likely reflects near-unidentifiability of
  `V_dnef` separately from the apparent formation clearance `K13`
  in the apparent-parameters parameterisation. The packaged model
  preserves the reported value (`omega^2 ~= 1e-4`) so a downstream
  user can see it explicitly rather than fixing IIV to zero.
- **Apparent-clearance parameterisation of K13 and V_dnef.** Both
  the metabolic formation clearance K13 and the desmethyl-nefopam
  distribution volume V_dnef carry the "apparent" qualifier in
  Djerada 2014 (implicit through the K13 label and the low V_dnef
  value). Under this parameterisation the fraction metabolised and
  the metabolite bioavailability are absorbed into the apparent
  coefficients; the parent's total CL_nef captures all elimination
  and is NOT reduced by the K13 flux to desmethyl-nefopam in the
  parent ODE. This is the standard Monolix parent-metabolite
  formulation.
- **Diagonal Omega structure.** The paper tested correlations
  between random effects (Methods p 1030) but reports the diagonal
  Omega as the final structure; the packaged model preserves the
  diagonal Omega with no off-diagonal covariance terms.
- **Table 2 CL comparison uses the weighted-global mean.** Djerada
  2014 does not report a global CL_nef mean in Table 2. The
  reference value used in the comparison table
  ((10 x 13.04 + 38 x 24.61) / 48 = 22.20 L/h) is the arithmetic
  weighted mean of the severe and non-severe subgroup means; this
  agrees within about 15% with the typical value 17.3 L/h once the
  lognormal mean-vs-median correction is applied
  (17.3 x exp(0.248 / 2) = 19.58 L/h).
- **Figure 5 replication uses per-subject empirical AUC and Cmax.**
  The paper's Figure 5 is drawn from 1000 Monte Carlo simulations
  in Monolix. The replication here uses a 200-subject virtual
  cohort with the same model and reads Cmax, tmax, and trapezoidal
  AUC directly from the simulated concentration-time curves; the
  same qualitative shape (Cmax and RCmax cut-value probabilities
  increasing sharply with duration, AUC cut-value probability
  roughly constant across durations) is reproduced.
- **Sampling grid.** The virtual cohort observation grid matches
  the paper's protocol (0, 5, 10, 15, 30, 60, 120, 180, 240, 480,
  1440 min) plus a denser 0.02-h / 0.25-h grid over the first 24 h
  for accurate Cmax / AUC extraction in the Figure 5 replication.
