Skip to contents

LNG-IUS 52 mg (Jensen 2023)

Source: vignettes/articles/Jensen_2023_lngIus52mg.Rmd
Jensen_2023_lngIus52mg.Rmd

Model and source

  • Citation: Jensen JT, Reinecke I, Post TM, Lukkari-Lax E, Hofmann BM. Extended use of levonorgestrel-releasing intrauterine system (LNG-IUS) 52 mg: A population pharmacokinetic approach to estimate in vivo levonorgestrel release rates and systemic exposure including comparison with two other LNG-IUSs. Contraception. 2023 May;121:109954. doi:10.1016/j.contraception.2023.109954. Upstream 5-year integrated LNG contraceptive popPK meta-analysis: Reinecke I, Hofmann B, Mesic E, Drenth HJ, Garmann D. J Clin Pharmacol. 2018 Dec;58(12):1639-1654. doi:10.1002/jcph.1288.
  • Description: Eight-year population PK / SHBG-turnover model for levonorgestrel (LNG) released from the LNG-releasing intrauterine system (LNG-IUS) 52 mg (Mirena). In-vivo release from the reservoir is described as a mixture of a first-order term (rate coefficient c12), a constant zero-order term (c13), and a time-dependent first-order term with time-decay constant t1, giving DADT(depot) = -c12 * depot - c13 * (1 + depot / (t1 + t)). Two- compartment LNG disposition treats only unbound drug as distributable / eliminable via K20 * fuLNG * central and K23 * fuLNG * central, where fuLNG is the closed-form free-fraction solution to reversible LNG binding to albumin (constant, KDA = 18209 nmol/L, ALB = 700000 nmol/L) and to SHBG (KDS = 1.82 nmol/L). SHBG serum concentration is modelled with an indirect-response turnover (zero-order kin, first-order kout) whose synthesis is linearly inhibited by a delay-compartment-smoothed unbound LNG signal (delay time-constant tau). Body weight enters as an allometric-style power on apparent (free) LNG clearance and on the SHBG baseline. Bioavailability of the loaded LNG reservoir is fixed via a logit anchor (F1 ~ 0.971 of the 52 mg reservoir). All rate constants are expressed per hour.
  • 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, 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.

mod_fn <- readModelDb("Jensen_2023_lngIus52mg")
mod    <- mod_fn()
str(mod$meta$population)
#> List of 12
#>  $ species       : chr "human"
#>  $ n_subjects    : int 920
#>  $ n_studies     : int 3
#>  $ studies       : chr [1:3] "NCT02985541 (Mirena Extension Trial, MET, Phase 3, US 54 centers, 2016-2021)" "NCT00185380 (Phase 2 Study 308901, LNG-IUS 52 mg, up to ~3 y)" "Leiras Study Report 02-89532-07 (Phase 3 Study 89532, LNG-IUS 52 mg, up to ~5 y)"
#>  $ age_range     : chr "18-40 years (medians 29-33 across the three studies)"
#>  $ weight_range  : chr "39-164 kg (medians 62 / 67 / 70 kg across studies; pooled median 68 kg)"
#>  $ sex_female_pct: num 100
#>  $ race_ethnicity: Named num [1:6] 75.4 14.1 2.5 0.6 3.9 3.6
#>   ..- attr(*, "names")= chr [1:6] "White" "Black" "Asian" "AmericanIndianAlaskaNative" ...
#>  $ disease_state : chr "Premenopausal fertile women using LNG-IUS 52 mg (Mirena) for contraception, with or without heavy menstrual bleeding."
#>  $ dose_range    : chr "Single insertion of a 52 mg LNG intrauterine reservoir (F1 ~ 0.971 of the loaded amount enters the release pool"| __truncated__
#>  $ regions       : chr "United States (MET); global for the pooled Phase 2 / 3 upstream studies."
#>  $ notes         : chr "Race percentages are baseline demographics for the 362 MET subjects starting extended treatment in year 6 (Supp"| __truncated__
str(mod$meta$covariateData)
#> List of 1
#>  $ WT:List of 6
#>   ..$ description       : chr "Body weight (kg). Continuous covariate on apparent clearance of LNG and on SHBG baseline via power-law scaling."
#>   ..$ units             : chr "kg"
#>   ..$ type              : chr "continuous"
#>   ..$ reference_category: NULL
#>   ..$ notes             : chr "Reference weight is 68 kg (medWT), the pooled median across the 8-year popPK dataset (median weight column in S"| __truncated__
#>   ..$ source_name       : chr "WGHTC"

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.

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 

# 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))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalrbase_shbg'
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.

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)."
  )
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).
Time point Release paper (ug/d) Release model (ug/d) Delta Release (%) LNG paper (ng/L) LNG model (ng/L) Delta LNG (%) Unbound paper (ng/L) Unbound model (ng/L) Delta Unbound (%) SHBG paper (nmol/L) SHBG model (nmol/L) Delta SHBG (%)
24 d 21.40 17.43 -18.5 175 193.3 +10.5 2.76 3.01 +8.9 41.9 45.5 +8.5
2 mo 21.00 16.98 -19.1 169 186.7 +10.5 2.68 2.92 +9.1 41.2 44.6 +8.3
3 mo 20.80 16.80 -19.2 167 184.8 +10.7 2.65 2.89 +9.2 41.2 44.7 +8.4
6 mo 20.10 16.43 -18.2 164 181.0 +10.4 2.60 2.83 +8.9 41.4 44.8 +8.3
9 mo 19.40 16.13 -16.9 161 177.8 +10.4 2.55 2.78 +8.9 41.5 44.9 +8.3
1 y 18.70 15.82 -15.4 159 174.6 +9.8 2.50 2.72 +9.0 41.6 45.1 +8.3
2 y 16.20 14.73 -9.1 149 163.2 +9.6 2.34 2.54 +8.4 41.9 45.5 +8.6
5 y 10.70 11.94 +11.5 123 133.6 +8.6 1.91 2.06 +7.6 42.9 46.7 +8.7
8 y 7.04 9.68 +37.4 100 109.2 +9.2 1.56 1.67 +6.8 43.2 47.6 +10.2

Figure 2 - release-rate curve over 8 years 

# 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 

# 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.

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).

# 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)
#> Warning: Requesting an AUC range starting (0) before the first measurement
#> (360) is not allowed
#> Warning: Requesting an AUC range starting (0) before the first measurement
#> (240) is not allowed
#> Warning: Requesting an AUC range starting (0) before the first measurement
#> (120) is not allowed

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)."
  )
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).
Window AUC(win) (ng*h/L) Cavg (ng/L)
year 5-6 NA NA
year 6-7 NA NA
year 7-8 NA NA

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:

    1. 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.