---
title: "LNG-IUS 52 mg (Jensen 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{LNG-IUS 52 mg (Jensen 2023)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r readModelDb("Jensen_2023_lngIus52mg")()$meta$reference`
* Description: `r readModelDb("Jensen_2023_lngIus52mg")()$meta$description`
* Article: <https://doi.org/10.1016/j.contraception.2023.109954> (open access, CC BY-NC-ND 4.0)
* ClinicalTrials.gov (Mirena Extension Trial, MET): <https://clinicaltrials.gov/study/NCT02985541>

Jensen et al. (2023) extend the LNG-IUS 52 mg (Mirena) popPK and release
model to 8 years of continuous use by combining Mirena Extension Trial (MET)
data (years 5 to 8) with two earlier LNG-IUS 52 mg trials (Phase 2 Study
308901 and Phase 3 Study 89532). The packaged model is the 8-year popPK
model whose NONMEM code appears in the paper's Supplementary Appendix. The
sibling 8-year release-only model (fit to residual-content data only) uses
different C12 / C13 / T1 estimates and is documented but not packaged; the
popPK model is preferred because it jointly explains LNG plasma, SHBG serum,
and residual-content observations. Time is expressed in hours throughout.

# Population

Pooled 920 premenopausal women from three studies: 361 in the MET (US
Phase 3, `r "NCT02985541"`, 2016-2021), 239 in a Phase 2 study of LNG-IUS
52 mg (NCT00185380 / Study 308901), and 320 in a Phase 3 five-year study
(Study 89532). Pooled median body weight was 68 kg (individual weights
39-164 kg); ages 18-40 years (medians 29-33). MET baseline race
composition was 75.4% White, 14.1% Black or African American, 2.5% Asian,
0.6% American Indian or Alaska Native, 3.9% multiple ethnicities, 3.6%
not reported (Supplementary Section 2.1); pooled race across all three
studies was not reported. All subjects received a single LNG-IUS 52 mg
device with intended use up to 8 years.

```{r metadata}
mod_fn <- readModelDb("Jensen_2023_lngIus52mg")
mod    <- mod_fn()
str(mod$meta$population)
str(mod$meta$covariateData)
```

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lc12  = log(7.98)`        | c12 = 7.98e-6 /hr        | Supplementary Appendix, `$THETA` TH1 |
| `lc13  = fixed(log(303))`  | c13 = 303e-6 mg/hr (FIXED) | Supplementary Appendix, `$THETA` TH11 FIX |
| `lt1   = log(58.8)`        | t1 = 58.8 hr             | Supplementary Appendix, `$THETA` TH12 |
| `lvc   = fixed(log(20.7))` | vc (V2) = 20.7 L (FIXED) | Supplementary Appendix, `$THETA` TH2 FIX |
| `lcl   = log(242)`         | cl (CL/F, apparent free) = 242 L/hr | Supplementary Appendix, `$THETA` TH3 |
| `lvp   = fixed(log(4690))` | vp (V3) = 4690 L (FIXED) | Supplementary Appendix, `$THETA` TH4 FIX |
| `lq    = fixed(log(600))`  | q (Q3) = 600 L/hr (FIXED) | Supplementary Appendix, `$THETA` TH5 FIX |
| `ltau  = fixed(log(13.7))` | tau = 13.7 hr (FIXED) | Supplementary Appendix, `$THETA` TH6 FIX |
| `lri   = fixed(log(0.232))` | ri = 0.232 L/nmol (FIXED) | Supplementary Appendix, `$THETA` TH7 FIX |
| `lrbase_shbg = log(51.8)`  | SBL = 51.8 nmol/L        | Supplementary Appendix, `$THETA` TH8 |
| `lkout_shbg  = fixed(log(0.00313))` | kout = 0.00313 /hr (FIXED) | Supplementary Appendix, `$THETA` TH9 FIX |
| `logitfdepot = fixed(3.51)` | F1 = plogis(3.51) = 0.971 (FIXED) | Supplementary Appendix, `$THETA` TH10 FIX |
| `e_wt_cl         = 0.728`  | WT exponent on CL        | Supplementary Appendix, `$THETA` TH13 |
| `e_wt_rbase_shbg = -0.987` | WT exponent on SBL       | Supplementary Appendix, `$THETA` TH14 |
| Reference weight medWT = 68 kg | scaling anchor for CO1 and CO2 | Supplementary Appendix, `$PK` block (`medWT = 68`) |
| Physicochemical constants MWLNG=312.5 g/mol, KDS=1.82 nmol/L, KDA=18209 nmol/L, ALB=700000 nmol/L | free-fraction of LNG closed-form | Supplementary Section 1.5 table ("Free fraction of LNG") and `$PK` block |
| `etalcl + etalrbase_shbg ~ c(0.0345, -0.0505, 0.173)` | OMEGA block (CL var, cov, SBL var) | Supplementary Appendix, `$OMEGA BLOCK(2)` |
| `propSd = sqrt(0.0329)`         | sqrt(0.0329) ~ 0.181 (proportional LNG) | Supplementary Appendix, `$SIGMA` row 1 |
| `propSd_shbg = sqrt(0.0344)`    | sqrt(0.0344) ~ 0.186 (proportional SHBG) | Supplementary Appendix, `$SIGMA` row 2 |
| `addSd_iusResidual = sqrt(1.52)` | sqrt(1.52) ~ 1.233 mg (additive residual content) | Supplementary Appendix, `$SIGMA` row 3 |
| `d/dt(depot) = -input1 - input2` with `input1 = c12*depot`, `input2 = c13*(1 + depot/(t1+t))` | mixed zero + first + time-dependent release | Supplementary Section 1.5, "LNG, release from IUS" row |
| `d/dt(central) = input1 + input2 - (k20+k23)*fuLNG*central + k32*peripheral1` | 2-compartment LNG with free-fraction on eliminative + distributive flows | Supplementary Section 1.5, "LNG, central compartment" row |
| `d/dt(peripheral1) = k23*fuLNG*central - k32*peripheral1` | LNG peripheral compartment | Supplementary Section 1.5, "LNG, peripheral compartment" row |
| `d/dt(effect) = (1/tau) * (A3nM - effect)` | delayed unbound-LNG signal driving SHBG inhibition | Supplementary Section 1.5, "LNG, delay compartment" row |
| `d/dt(shbg) = kin*(1-inh) - kout*shbg`; `kin = SBL*kout`; `inh = min(ri*effect, 1)` | SHBG indirect-response turnover | Supplementary Section 1.5, "SHBG" row |
| `shbg(0) = SBL * (WT/68)^e_wt_rbase_shbg` | SHBG initial condition at insertion | Supplementary Appendix, `$PK` block (`A_INITIAL(6) = SBL`) |

# Virtual cohort

Original observed data are not publicly available. The vignette uses a
typical-value simulation at the reference weight of 68 kg and a
body-weight subgroup analysis (`<=55` kg vs `>55` kg) matching
Supplementary Figure 1. The stochastic-VPC cohort uses 200 subjects per
arm at WT = 68 kg drawn from log-normal individual-parameter
distributions.

```{r cohort}
set.seed(20260724L)

# Time grid: dense early (release dynamics + fast SHBG turnover), then
# monthly through year 8. Times are in HOURS throughout.
sample_days <- c(seq(0.5, 14, by = 0.5),
                 seq(15, 60, by = 3),
                 seq(70, 365, by = 10),
                 seq(400, 8 * 365, by = 30))
sample_hr   <- sample_days * 24

# Typical-value trajectory at reference WT = 68 kg. Observation rows use
# dvid = 1 (the internal rxode2 routing for the Cc observation output);
# this avoids referencing the algebraic observable name in the event
# table's cmt column, per the vignette lint. The rxode2 solve emits all
# algebraic observables (Cc, CcUnbound, shbg, iusResidual) as output
# columns regardless of the dvid choice on the observation rows.
events_typical <- data.frame(
  id   = 1L,
  time = c(0, sample_hr),
  amt  = c(52, rep(0, length(sample_hr))),
  evid = c(1L, rep(0L, length(sample_hr))),
  cmt  = c("depot", rep(NA_character_, length(sample_hr))),
  dvid = c(NA_integer_, rep(1L, length(sample_hr)))
)

# Body-weight subgroup: sample n_arm = 24 (<=55 kg) and n_ref = 200 (>55 kg)
# to match Supp Fig 1 sizes; log-normal individual etas.
n_arm  <- 24L
n_ref  <- 200L

make_bw_cohort <- function(wt_vec, arm_label, id_offset) {
  ids <- id_offset + seq_along(wt_vec)
  dose_rows <- data.frame(
    id = ids, time = 0, amt = 52, cmt = "depot",
    evid = 1L, dvid = NA_integer_, WT = wt_vec, arm = arm_label
  )
  obs_grid <- expand.grid(id = ids, time = sample_hr)
  obs_rows <- obs_grid |>
    dplyr::mutate(amt = 0, cmt = NA_character_, evid = 0L, dvid = 1L) |>
    dplyr::left_join(data.frame(id = ids, WT = wt_vec, arm = arm_label),
                     by = "id")
  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

# Draw weight distributions consistent with Supp Fig 1 groupings (LNG-IUS 52 mg
# MET cohort mean BMI 27.9 kg/m^2). The vignette only needs cohort-mean weight
# behaviour, not the exact study distribution; using deterministic weight
# vectors avoids run-to-run variability while still exercising the WT
# covariate.
wt_ref <- rep(70, n_ref)                                # >55 kg reference
wt_low <- seq(42, 55, length.out = n_arm)               # <=55 kg subgroup

events_wt <- dplyr::bind_rows(
  make_bw_cohort(wt_ref, arm_label = ">55 kg (typical)", id_offset =   0L),
  make_bw_cohort(wt_low, arm_label = "<=55 kg (subgroup)", id_offset = n_ref)
)
```

# Simulation

```{r simulate}
# readModelDb returns the model function; call it to build the rxUi.
mod         <- readModelDb("Jensen_2023_lngIus52mg")()
mod_typical <- rxode2::zeroRe(mod)

# Typical-value trajectory (reference WT = 68 kg).
sim_typ <- rxode2::rxSolve(mod_typical, events_typical, params = c(WT = 68))
sim_typ$day    <- sim_typ$time / 24
sim_typ$year   <- sim_typ$day / 365.25

# Compute the instantaneous release rate as input1 + input2 (mg/hr) using the
# packaged typical values, then convert to ug/day for the Table 1 comparison.
th <- mod$iniDf
theta <- setNames(th$est, th$name)
c12_typ <- exp(theta[["lc12"]]) * 1e-6
c13_typ <- exp(theta[["lc13"]]) * 1e-6
t1_typ  <- exp(theta[["lt1"]])

sim_typ$release_rate_ug_per_d <- with(sim_typ,
  1e3 * 24 * (c12_typ * depot + c13_typ * (1 + depot / (t1_typ + time)))
)

# Stochastic-VPC cohort with IIV, WT subgroup.
sim_wt <- rxode2::rxSolve(mod, events_wt, keep = c("arm", "WT"))
sim_wt$day  <- sim_wt$time / 24
sim_wt$year <- sim_wt$day / 365.25
```

# Replicate published typical values (Table 1 and Figure 2)

The typical model prediction is compared against the paper's Table 1
values for LNG-IUS 52 mg at nine time points from 24 days through 8
years.

```{r table1_compare}
paper_table1 <- tibble::tribble(
  ~day,   ~label,    ~release_rate_ug_d, ~Cc_ng_L, ~CcUnbound_ng_L, ~shbg_nmol_L,
  24,     "24 d",    21.4,               175,      2.76,            41.9,
  60,     "2 mo",    21.0,               169,      2.68,            41.2,
  90,     "3 mo",    20.8,               167,      2.65,            41.2,
  180,    "6 mo",    20.1,               164,      2.60,            41.4,
  270,    "9 mo",    19.4,               161,      2.55,            41.5,
  365,    "1 y",     18.7,               159,      2.50,            41.6,
  730,    "2 y",     16.2,               149,      2.34,            41.9,
  1825,   "5 y",     10.7,               123,      1.91,            42.9,
  2920,   "8 y",     7.04,               100,      1.56,            43.2
)

# Interpolate the typical simulation at the requested times.
interp_at <- function(target_day, col) {
  approx(x = sim_typ$day, y = sim_typ[[col]], xout = target_day)$y
}

comparison <- paper_table1 |>
  dplyr::mutate(
    release_pred = interp_at(day, "release_rate_ug_per_d"),
    Cc_pred      = interp_at(day, "Cc"),
    CcUnbound_pred = interp_at(day, "CcUnbound"),
    shbg_pred    = interp_at(day, "shbg")
  ) |>
  dplyr::mutate(
    release_pct_diff  = 100 * (release_pred - release_rate_ug_d) / release_rate_ug_d,
    Cc_pct_diff       = 100 * (Cc_pred - Cc_ng_L) / Cc_ng_L,
    CcUnbound_pct_diff = 100 * (CcUnbound_pred - CcUnbound_ng_L) / CcUnbound_ng_L,
    shbg_pct_diff     = 100 * (shbg_pred - shbg_nmol_L) / shbg_nmol_L
  )

comparison |>
  dplyr::transmute(
    "Time point"           = label,
    "Release paper (ug/d)" = round(release_rate_ug_d, 2),
    "Release model (ug/d)" = round(release_pred, 2),
    "Delta Release (%)"    = sprintf("%+.1f", release_pct_diff),
    "LNG paper (ng/L)"     = round(Cc_ng_L, 1),
    "LNG model (ng/L)"     = round(Cc_pred, 1),
    "Delta LNG (%)"        = sprintf("%+.1f", Cc_pct_diff),
    "Unbound paper (ng/L)" = round(CcUnbound_ng_L, 2),
    "Unbound model (ng/L)" = round(CcUnbound_pred, 2),
    "Delta Unbound (%)"    = sprintf("%+.1f", CcUnbound_pct_diff),
    "SHBG paper (nmol/L)"  = round(shbg_nmol_L, 1),
    "SHBG model (nmol/L)"  = round(shbg_pred, 1),
    "Delta SHBG (%)"       = sprintf("%+.1f", shbg_pct_diff)
  ) |>
  knitr::kable(
    caption = "Replicates Table 1 of Jensen 2023 for LNG-IUS 52 mg. Signed percent differences show a systematic +10 to +14% overprediction of the LNG plasma concentration by the extracted popPK model at reference WT = 68 kg (see Assumptions and deviations)."
  )
```

## Figure 2 - release-rate curve over 8 years

```{r figure2}
# Replicates Figure 2 of Jensen 2023 for LNG-IUS 52 mg release rate.
paper_pts <- paper_table1 |>
  dplyr::transmute(year = day / 365.25,
                   release_rate_ug_per_d = release_rate_ug_d,
                   source = "Paper Table 1")

ggplot() +
  geom_line(data = dplyr::filter(sim_typ, year <= 8.1),
            aes(year, release_rate_ug_per_d, colour = "Extracted popPK model"),
            linewidth = 0.8) +
  geom_point(data = paper_pts,
             aes(year, release_rate_ug_per_d, colour = source), size = 2) +
  scale_colour_manual(values = c("Extracted popPK model" = "steelblue",
                                 "Paper Table 1" = "black")) +
  scale_y_continuous(limits = c(0, 25)) +
  labs(x = "Years after insertion",
       y = "Typical LNG release rate (ug/day)",
       colour = NULL,
       title = "Figure 2 - LNG-IUS 52 mg release rate over 8 years",
       caption = "Replicates Figure 2 of Jensen 2023 (LNG-IUS 52 mg trace).")
```

## Figure 3 - LNG, unbound LNG, and SHBG trajectories

```{r figure3}
# Replicates Figure 3 of Jensen 2023 (panels A, B, C for LNG-IUS 52 mg).
long_typ <- sim_typ |>
  dplyr::filter(year <= 8.1) |>
  dplyr::select(year, Cc, CcUnbound, shbg) |>
  tidyr::pivot_longer(-year, names_to = "series", values_to = "value") |>
  dplyr::mutate(
    panel = dplyr::case_when(
      series == "Cc"        ~ "A: Total LNG (ng/L)",
      series == "CcUnbound" ~ "B: Unbound LNG (ng/L)",
      series == "shbg"      ~ "C: SHBG (nmol/L)"
    )
  )

paper_long <- paper_table1 |>
  dplyr::transmute(
    year = day / 365.25,
    "A: Total LNG (ng/L)"   = Cc_ng_L,
    "B: Unbound LNG (ng/L)" = CcUnbound_ng_L,
    "C: SHBG (nmol/L)"      = shbg_nmol_L
  ) |>
  tidyr::pivot_longer(-year, names_to = "panel", values_to = "value")

ggplot(long_typ, aes(year, value)) +
  geom_line(aes(colour = "Extracted popPK model"), linewidth = 0.8) +
  geom_point(data = paper_long,
             aes(year, value, colour = "Paper Table 1"), size = 2) +
  facet_wrap(~ panel, scales = "free_y", ncol = 1) +
  scale_colour_manual(values = c("Extracted popPK model" = "steelblue",
                                 "Paper Table 1" = "black")) +
  labs(x = "Years after insertion", y = NULL, colour = NULL,
       title = "Figure 3 - LNG-IUS 52 mg total LNG, unbound LNG, and SHBG",
       caption = "Replicates Figure 3 of Jensen 2023 (LNG-IUS 52 mg traces).")
```

## Supplementary Figure 1 - body-weight subgroup

The paper's subgroup analysis notes that women weighing `<=55` kg had LNG
concentrations 27-43% higher than women weighing `>55` kg during years
5-8. The packaged model reproduces this qualitatively via the
(WT / 68)^0.728 scaling on apparent clearance.

```{r suppfig1}
sim_wt_summary <- sim_wt |>
  dplyr::filter(year >= 5, year <= 8, !is.na(Cc)) |>
  dplyr::group_by(arm, year) |>
  dplyr::summarise(Cc_geo = exp(mean(log(pmax(Cc, 1e-6)))),
                   .groups = "drop")

sim_wt_summary |>
  ggplot(aes(year, Cc_geo, colour = arm)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Years after insertion",
       y = "Typical LNG (ng/L, geometric mean)",
       colour = "Body weight group",
       title = "Supp Fig 1 - LNG-IUS 52 mg by body weight subgroup",
       caption = "Replicates Supplementary Figure 1 of Jensen 2023 (LNG-IUS 52 mg by body-weight group).")
```

# PKNCA validation - year-8 window average concentration

For sustained-release intrauterine devices the natural PKNCA parameter
is the average concentration `Cavg` over a defined window (year of use).
This vignette computes `Cavg = AUC / (t_end - t_start)` over year 8
(days 2555 to 2920) from the typical-value simulation, and compares
against the paper's Table 1 8-year geometric-mean value of 100 ng/L. The
PKNCA formula includes a `treatment` grouping variable per convention
(only one group here, `"typical"`, but the group column is required).

```{r pknca}
# Concentrations. Keep the column named Cc.
sim_nca_input <- sim_typ |>
  dplyr::mutate(id = 1L, treatment = "typical") |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# Guarantee a time = 0 row (see pknca-recipes.md, "Time-zero records
# (mandatory)"). For LNG-IUS the depot is inserted at t = 0 with no
# preloaded plasma exposure, so Cc(0) = 0.
sim_nca_input <- dplyr::bind_rows(
  sim_nca_input,
  data.frame(id = 1L, time = 0, Cc = 0, treatment = "typical")
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca_input, Cc ~ time | treatment + id)

# One dose per subject; PKNCA needs the dose row for calculating Cavg
# over user-defined intervals via `AUC / (t_end - t_start)`.
dose_df <- data.frame(
  id = 1L, treatment = "typical",
  time = 0, amt = 52
)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

# Interval: year 8 in hours (5-6 year window, 6-7, and 7-8 for context).
year_hours <- function(y_start, y_end) c(y_start * 365 * 24, y_end * 365 * 24)

intervals <- data.frame(
  start   = c(5, 6, 7) * 365 * 24,
  end     = c(6, 7, 8) * 365 * 24,
  cav     = TRUE,
  auclast = TRUE
)

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

nca_wide <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD %in% c("cav", "auclast")) |>
  dplyr::mutate(year_window = sprintf("year %d-%d", start / (365 * 24),
                                                     end / (365 * 24))) |>
  dplyr::select(year_window, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

nca_wide |>
  dplyr::rename(
    "Window"          = year_window,
    "Cavg (ng/L)"     = cav,
    "AUC(win) (ng*h/L)" = auclast
  ) |>
  knitr::kable(
    digits  = c(0, 1, 0),
    caption = "PKNCA Cavg and AUC of typical LNG plasma concentration over years 5-6, 6-7, and 7-8 (h scale). Year-8 window Cavg from the model can be compared against the paper's Table 1 8-year geometric-mean LNG of 100 ng/L; the model overpredicts by roughly 10% at typical WT = 68 kg (see Assumptions and deviations)."
  )
```

## Comparison against published NCA

The paper does not report a strict `Cmax` / `Tmax` / `AUC0-inf` NCA
because the LNG-IUS device is a multi-year sustained-release
formulation. Instead the paper reports geometric-mean LNG concentration
at fixed calendar time points (Table 1) and the average LNG-IUS 52 mg
release rate over selected windows (Table 1 "Average over ..." rows).
The comparison for those metrics is in the *Replicate published typical
values* table above. Cavg over year 8 from the packaged model is close
to the paper's year-8 geometric mean of 100 ng/L, with the systematic
+10% bias discussed below.

# Assumptions and deviations

- **Typical-value bias of +10% to +14% on total LNG plasma
  concentration** (see the Replicate published typical values table).
  My extracted 8-year popPK model consistently overpredicts the paper's
  Table 1 typical LNG concentrations by roughly 10-14% across all time
  points from 24 days through 8 years, while matching the SHBG
  trajectory, free-fraction fu = 1.5-1.6% (paper reports ~1.6%), and
  8-year residual content within ~7% (12.5 mg model vs ~13.5 mg paper).
  Candidate causes considered and ruled out during extraction:
  (a) numerical ODE tolerance (tightening `atol = 1e-10, rtol = 1e-9`
  in `rxode2::rxSolve` did not change the trajectory); (b) initial
  condition setup for the SHBG compartment (both the paper's `A_INITIAL(6) =
  SBL` and the packaged `shbg(0) <- rbase_shbg` set SHBG to its
  weight-adjusted baseline at insertion); (c) misparsed 10^-6 scale on
  the C12 / C13 release parameters (release-rate mass balance over 8
  years matches the paper's 13 ug/day average exactly). The most likely
  remaining cause is a small (~5%) difference between the NONMEM $THETA
  values printed in the Supplementary Appendix and the values actually
  used to render Figures 2-3 and Table 1 of the main paper. Bayesian
  point estimates and confidence intervals were reported with 3
  significant figures throughout, so parameter-rounding effects
  compound. Downstream users needing exact paper-typical trajectories
  should be aware of this ~10% offset; the packaged model reproduces
  the model *structure* faithfully.

- **Time unit is hours.** The NONMEM control stream in the
  Supplementary Appendix does not explicitly declare a time unit; hours
  is inferred from the SHBG elimination-rate kout = 0.00313, which
  gives a half-life of log(2) / 0.00313 = 221 h (~9 days) if kout is
  per hour but 221 days (~7 months) if per day. The 9-day SHBG
  half-life is consistent with published SHBG turnover; the 7-month
  half-life is not physiological. All rate constants are therefore
  expressed per hour in the packaged model.

- **Reference weight fixed at 68 kg.** The NONMEM control stream fixes
  `medWT = 68` unconditionally and imputes missing / zero WT to 68 kg.
  Pooled population median is closer to 66 kg (weighted average of
  study medians 70 / 67 / 62 kg from Supplementary Table 1); the
  packaged model uses 68 kg exactly as coded in the Appendix.

- **Weight covariate distributions in the vignette body-weight
  subgroup analysis are deterministic** (WT = 70 kg for the reference
  arm and WT = 42-55 kg spaced linearly for the low-weight arm),
  rather than sampled from an empirical distribution. The paper does
  not publish the individual weight distribution, so any specific
  sampling scheme would be a model of a model; the deterministic
  spacing gives a reproducible cohort-mean behaviour without
  fabricating subject-level data.

- **The 8-year release-only sibling model is not packaged.** Jensen
  2023 fits both a joint popPK model (LNG + SHBG + residual content,
  packaged here) and a release-only model (residual content only, with
  different C12 / C13 / T1 estimates). Table 1 typical release rates
  are computed from the release-only model, whereas Table 1 typical
  concentrations are computed from the popPK model. The packaged
  model's release rates therefore differ from Table 1 at each
  instantaneous time point (~15-18% lower at 15 days, matching by 8
  years) even though the depot-content trajectory matches within ~7%.
  The release-only model is documented in Supplementary Section 1.6
  and its Appendix code, but not exported to `modellib()` because it
  duplicates the release process of the popPK model.

- **The paper's proprietary NM.Mirena.E.v17.csv dataset is not
  available.** All simulations are typical-value or synthetic-cohort;
  Visual Predictive Checks against measured concentrations require the
  proprietary MET dataset.
