---
title: "Lamotrigine + N-2-glucuronide (Milosheska 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Lamotrigine + N-2-glucuronide (Milosheska 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

```{r model-meta, echo = FALSE}
.mod_meta <- rxode2::rxode(readModelDb("Milosheska_2016_lamotrigine"))
```

* Citation: `r .mod_meta$meta$reference`
* Description: `r .mod_meta$meta$description`
* Article: <https://doi.org/10.1111/bcp.12984>

The packaged model implements the Milosheska 2016 parent + metabolite joint
population-PK model for oral lamotrigine (LTG) and its N-2-glucuronide
(LTG-glu) in adult epilepsy patients. The parent side is a one-compartment
open-kinetic structure with first-order absorption; the metabolite side is
a one-compartment linear-elimination compartment fed by a
complete-conversion mass-balance arrow from the parent's total elimination
flux. Parent apparent oral clearance (`cl` = CL/F) is modified
multiplicatively by nine covariate effects (Table 4 footnote formula);
parent apparent volume (`vc` = V/F) is modified linearly by body-weight
deviation from 70 kg. Metabolite apparent clearance (`cl_gluc` = CL_LTG-glu
/ F_metab) carries a power weight and linear Cockcroft-Gault CLcr effect;
metabolite apparent volume (`vc_gluc` = V_LTG-glu / F_metab) has no
between-subject variability supported by the sparse metabolite data. The
paper's Table 4 is the parent-only final model; the joint parent +
metabolite fit reported "very similar" LTG parameters (Results paragraph 1
of the LTG-glu section) without a separate tabulation, so the packaged
model uses Table 4 values for the LTG side and the prose-reported values
for the LTG-glu side.

# Population

Milosheska 2016 enrolled 100 adult epilepsy patients on stable lamotrigine
therapy for at least two months at the Department of Neurology of the
University Medical Centre Ljubljana, Slovenia (Milosheska 2016 Table 1).
Median age was 39.8 years (range 20.7 to 80.4), median total body weight
was 70 kg (range 50 to 124), 70 of 100 were female, and median
Cockcroft-Gault CLcr was 107.6 mL/min (range 40.84 to 246). Fifty-four
patients were on lamotrigine monotherapy and 46 were on adjunctive therapy
with other antiepileptic drugs. Enzyme-inducing AED (EIAED) coadministration
was pooled by the authors into a single indicator covering carbamazepine
(n = 4), phenobarbital (n = 2), and phenytoin (n = 2), with oral
contraceptives (n = 2) excluded due to small N. UGT-inhibitor coadministration
was pooled as a single indicator covering valproic acid (n = 13) and
sertraline (n = 2). Two blood samples per patient (approximating trough
and peak at steady state) were drawn for a total of 195 lamotrigine and
195 lamotrigine-N-2-glucuronide plasma concentrations. Assay range for
both analytes was 0.25 mg/L to about 20 mg/L (Milosheska 2016 Methods
'Drug assay'). Prior population-pharmacokinetic estimates from a
meta-analysis of eleven published lamotrigine studies were incorporated
via NONMEM PRIOR functionality to stabilise Ka and V estimation, as
detailed in Milosheska 2016 Methods 'Meta-analysis of previous population
pharmacokinetic studies with LTG'. Genotype frequencies for all seven
tested variants are reported in Milosheska 2016 Table 2. The same
metadata is available programmatically via
`readModelDb("Milosheska_2016_lamotrigine")$meta$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Milosheska_2016_lamotrigine.R`. The table below
collects them in one place for review. Final point estimates for the
parent side are from Milosheska 2016 Table 4 (the paper's parent-only
final model); metabolite parameters are from the Results
'Population pharmacokinetic analysis of the lamotrigine-N-2-glucuronide'
section (prose only, no separate table).

| Equation / parameter    | Value                     | Source location |
|-------------------------|---------------------------|-----------------|
| `lka` (Ka)              | log(1.96) -> 1.96 1/h     | Table 4 row 'Ka': 1.96 (95% CI 1.72 - 2.24) |
| `lcl` (CL/F)            | log(2.40) -> 2.40 L/h     | Table 4 row 'CL': 2.40 (95% CI 2.30 - 2.50) |
| `lvc` (V/F)             | log(76.2) -> 76.2 L       | Table 4 row 'V': 76.2 (95% CI 66.9 - 84.2) |
| `lcl_gluc`              | log(3.16) -> 3.16 L/h     | Results LTG-glu paragraph 1: 3.16 L/h |
| `lvc_gluc`              | log(110) -> 110 L         | Results LTG-glu paragraph 1: 110 L |
| `e_conmed_ugt_inh_cl`   | -0.579                    | Table 4 row 'Co-treatment with inhibitors': -0.579 |
| `e_wt_cl` (power)       | 0.938                     | Table 4 row 'Body weight' on CL: 0.938 |
| `e_smoke_cl`            | +0.340                    | Table 4 row 'Cigarette smoking': +0.340 |
| `e_conmed_eiaed_cl`     | +0.546                    | Table 4 row 'Co-treatment with inducers': +0.546 |
| `e_ugt2b7_m161ct_cl`    | -0.0358                   | Table 4 row 'UGT2B7 -161 CT vs CC': -0.0358 |
| `e_ugt2b7_m161tt_cl`    | -0.204                    | Table 4 row 'UGT2B7 -161 TT vs CC': -0.204 |
| `e_ugt2b7_372ag_cl`     | +0.194                    | Table 4 row 'UGT2B7 372 AG vs AA': +0.194 |
| `e_ugt2b7_372gg_cl`     | +1.17                     | Table 4 row 'UGT2B7 372 GG vs AA': +1.17 |
| `e_crcl_cl`             | +0.00328 per mL/min       | Table 4 row 'CLcr': +0.00328 per 1 mL/min deviation from 110 |
| `e_wt_vc`               | +0.0181 per kg            | Table 4 row 'Body weight' on V: +0.0181 per 1 kg deviation from 70 |
| `e_wt_cl_gluc` (power)  | 1.01                      | Results LTG-glu paragraph 1: power exponent 1.01 on weight |
| `e_crcl_cl_gluc`        | +0.00759 per mL/min       | Results LTG-glu paragraph 1: -0.759% per 1 mL/min decrease in CLcr |
| `etalka` (omega^2)      | log(1 + 0.711^2) = 0.4090 | Table 4 row 'IIV Ka (%)': 71.1 |
| `etalcl` (omega^2)      | log(1 + 0.331^2) = 0.1040 | Table 4 row 'IIV CL (%)': 33.1 |
| `etalvc` (omega^2)      | log(1 + 0.301^2) = 0.0867 | Table 4 row 'IIV V (%)': 30.1 |
| `etalcl_gluc` (omega^2) | log(1 + 0.417^2) = 0.1603 | Results LTG-glu paragraph 1: 'Unexplained IIV on CL_LTG-glu was 41.7%' |
| `propSd` (parent)       | 0.180                     | Table 4 row 'Residual variability (%)': 18.0 |
| `propSd_gluc`           | 0.138                     | Results LTG-glu paragraph 1: 'residual ... variability of LTG-glu concentration was 13.8%' |
| ODE `d/dt(depot)`       | `-ka * depot`             | Methods 'Model development' page 4 (ADVAN 6) |
| ODE `d/dt(central)`     | `ka * depot - cl * Cc`    | Methods (parent one-compartment linear elimination) |
| ODE `d/dt(central_gluc)`| `cl * Cc - cl_gluc * Cc_gluc` | Methods (complete parent-to-metabolite conversion, ADVAN 6) |

# Virtual cohort

The Milosheska 2016 individual-patient data are not publicly available.
The virtual cohort below reproduces the paper's demographic distribution
and stratifies simulations across four illustrative subgroups so all
covariate effects can be inspected on the same axes:

* `Reference`: 70 kg non-smoker, no comedications, CLcr = 110 mL/min,
  UGT2B7 -161 CC / 372 AA (the paper's reference stratum).
* `+Inducer`: same as Reference but with `CONMED_EIAED = 1`
  (Table 4 row 'Co-treatment with inducers'; +54.6% CL/F).
* `+Inhibitor`: same as Reference but with `CONMED_UGT_INH = 1`
  (Table 4 row 'Co-treatment with inhibitors'; -57.9% CL/F).
* `+UGT2B7 372 GG`: same as Reference but with `UGT2B7_372GG = 1`
  (Table 4 row 'UGT2B7 372 GG vs AA'; +117% CL/F; the largest single-variant
  effect in the paper).

Each arm carries 200 subjects sampled from a log-normal weight distribution
centred on 70 kg with SD chosen to reproduce the paper's 50 to 124 kg range
approximately. Every subject receives 100 mg oral lamotrigine BID
(200 mg/day, the paper's median daily dose per Table 1), with dosing over
15 half-lives to reach steady state and a fine 0.5 h observation grid over
the final 12 h interval for NCA.

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

mod <- rxode2::rxode(readModelDb("Milosheska_2016_lamotrigine"))

n_per_arm <- 200L
wt_ref    <- 70
crcl_ref  <- 110
amt_mg    <- 100
tau_h     <- 12
n_doses   <- 30L
last_dose_time <- (n_doses - 1L) * tau_h

# 15 half-lives at Ka = 1.96 1/h is ~5 h; the CL/vc time constant at the
# reference stratum is 76.2/2.40 = 31.75 h, giving t1/2 ~ 22 h. Thirty BID
# doses spans 15 days, well beyond 5 t1/2 for the parent (110 h) and about
# 15 t1/2 for the metabolite (110/3.16 * ln(2) ~ 24 h -> 15 days = 15
# doses' worth of accumulation).
make_arm <- function(n, arm_label,
                     conmed_ugt_inh = 0L, smoke = 0L,
                     conmed_eiaed   = 0L,
                     ugt2b7_m161ct  = 0L, ugt2b7_m161tt = 0L,
                     ugt2b7_372ag   = 0L, ugt2b7_372gg  = 0L,
                     crcl_center    = crcl_ref,
                     id_offset      = 0L) {
  wt   <- pmin(pmax(exp(stats::rnorm(n, mean = log(wt_ref), sd = 0.20)),
                    50), 124)
  crcl <- pmin(pmax(stats::rnorm(n, mean = crcl_center, sd = 30),
                    40), 246)
  base <- tibble::tibble(
    id             = id_offset + seq_len(n),
    WT             = wt,
    CRCL           = crcl,
    SMOKE          = smoke,
    CONMED_EIAED   = conmed_eiaed,
    CONMED_UGT_INH = conmed_ugt_inh,
    UGT2B7_M161CT  = ugt2b7_m161ct,
    UGT2B7_M161TT  = ugt2b7_m161tt,
    UGT2B7_372AG   = ugt2b7_372ag,
    UGT2B7_372GG   = ugt2b7_372gg,
    arm            = arm_label,
    amt_per_dose   = amt_mg
  )
  doses <- tidyr::crossing(id = base$id, dose_idx = seq_len(n_doses)) |>
    dplyr::mutate(time = (dose_idx - 1L) * tau_h,
                  evid = 1L,
                  cmt  = "depot") |>
    dplyr::left_join(base, by = "id") |>
    dplyr::mutate(amt = amt_per_dose)
  obs_grid <- c(0, seq(last_dose_time + 0.25, last_dose_time + tau_h,
                       by = 0.5))
  # Setting `cmt = "Cc"` on observation rows is enough for this two-output
  # model: rxode2 computes both Cc and Cc_gluc at every solve step and
  # returns both columns in the simulation output. See the Abduljalil 2009
  # clarithromycin vignette for the same pattern.
  obs <- tidyr::crossing(id = base$id, time = obs_grid) |>
    dplyr::mutate(evid = 0L,
                  amt  = 0,
                  cmt  = "Cc") |>
    dplyr::left_join(base, by = "id")
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

arms <- dplyr::bind_rows(
  make_arm(n_per_arm, "Reference",         id_offset = 0L * n_per_arm),
  make_arm(n_per_arm, "+Inducer",          conmed_eiaed   = 1L,
           id_offset = 1L * n_per_arm),
  make_arm(n_per_arm, "+Inhibitor",        conmed_ugt_inh = 1L,
           id_offset = 2L * n_per_arm),
  make_arm(n_per_arm, "+UGT2B7 372 GG",    ugt2b7_372gg   = 1L,
           id_offset = 3L * n_per_arm)
)
stopifnot(!anyDuplicated(unique(arms[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
sim_pop <- rxode2::rxSolve(
  object     = mod,
  events     = arms,
  keep       = c("arm", "WT", "CRCL", "SMOKE", "CONMED_EIAED",
                 "CONMED_UGT_INH", "UGT2B7_M161CT", "UGT2B7_M161TT",
                 "UGT2B7_372AG", "UGT2B7_372GG", "amt_per_dose"),
  returnType = "data.frame"
)
```

For deterministic typical-value replication of the mean profile across
each arm, zero out between-subject variability:

```{r simulate-typical}
sim_typ <- rxode2::rxSolve(
  object     = rxode2::zeroRe(mod),
  events     = arms,
  keep       = c("arm", "WT", "CRCL", "SMOKE", "CONMED_EIAED",
                 "CONMED_UGT_INH", "UGT2B7_M161CT", "UGT2B7_M161TT",
                 "UGT2B7_372AG", "UGT2B7_372GG", "amt_per_dose"),
  returnType = "data.frame"
)
```

# Replicate published figures

## Steady-state Cc profile by arm

Milosheska 2016 does not publish a per-subject or per-covariate visual
predictive check (the paper's Figure 3 shows a whole-population VPC).
The chunk below plots the steady-state parent lamotrigine concentration
across the last dosing interval, stratified by covariate arm, to make
the multiplicative CL effects of the Table 4 footnote formula directly
visible.

```{r css-vpc, fig.width = 7.5, fig.height = 4.2}
last_iv <- sim_pop |>
  dplyr::filter(time >= last_dose_time,
                time <= last_dose_time + tau_h) |>
  dplyr::mutate(tau_time = time - last_dose_time)

last_iv |>
  dplyr::group_by(arm, tau_time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(tau_time, Q50, colour = arm, fill = arm)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), colour = NA, alpha = 0.20) +
  geom_line(linewidth = 0.7) +
  labs(x = "Time in dosing interval (h)",
       y = "Parent LTG Cc (mg/L)",
       colour = "Arm", fill = "Arm",
       title  = "Steady-state parent lamotrigine Cc by covariate arm",
       caption = "100 mg BID; 200 subjects per arm; ribbon = 5th-95th percentile.")
```

## Steady-state LTG-glu profile by arm

```{r css-vpc-gluc, fig.width = 7.5, fig.height = 4.2}
last_iv |>
  dplyr::group_by(arm, tau_time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc_gluc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc_gluc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc_gluc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(tau_time, Q50, colour = arm, fill = arm)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), colour = NA, alpha = 0.20) +
  geom_line(linewidth = 0.7) +
  labs(x = "Time in dosing interval (h)",
       y = "Metabolite LTG-glu Cc_gluc (mg/L)",
       colour = "Arm", fill = "Arm",
       title  = "Steady-state N-2-glucuronide Cc_gluc by covariate arm",
       caption = "Same simulation as the parent panel above; the metabolite compartment is fed by the complete-conversion mass-balance arrow.")
```

# PKNCA validation

Steady-state NCA is computed on the typical-value simulation over the last
BID dosing interval, per arm. The paper does not report per-subject NCA
metrics, so the check below is a self-consistency audit: PKNCA output
should reproduce the Table 4 footnote covariate multipliers when Cmax and
Cavg across arms are compared against the reference arm.

```{r pknca-parent}
# Re-anchor the last interval to time = 0 for PKNCA.
pk_conc_parent <- sim_typ |>
  dplyr::filter(time >= last_dose_time,
                time <= last_dose_time + tau_h) |>
  dplyr::mutate(time = time - last_dose_time) |>
  dplyr::select(id, time, Cc, arm, WT, CRCL) |>
  dplyr::filter(!is.na(Cc))

# Guarantee a time = 0 row per (id, arm) so PKNCA anchors AUC at t = 0.
pk_conc_parent <- dplyr::bind_rows(
  pk_conc_parent,
  pk_conc_parent |>
    dplyr::distinct(id, arm) |>
    dplyr::mutate(time = 0)
) |>
  dplyr::distinct(id, arm, time, .keep_all = TRUE) |>
  dplyr::arrange(id, arm, time) |>
  dplyr::mutate(Cc = dplyr::coalesce(Cc, 0))

pk_dose_parent <- pk_conc_parent |>
  dplyr::distinct(id, arm) |>
  dplyr::mutate(time = 0, amt = amt_mg)

conc_obj_parent <- PKNCA::PKNCAconc(
  pk_conc_parent, Cc ~ time | arm + id,
  concu = "mg/L", timeu = "hr"
)
dose_obj_parent <- PKNCA::PKNCAdose(
  pk_dose_parent, amt ~ time | arm + id, doseu = "mg"
)

intervals_ss <- data.frame(
  start   = 0,
  end     = tau_h,
  cmax    = TRUE,
  tmax    = TRUE,
  auclast = TRUE,
  cmin    = TRUE,
  cav     = TRUE
)

nca_parent <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_parent, dose_obj_parent, intervals = intervals_ss)
)
knitr::kable(
  summary(nca_parent),
  caption = paste("Steady-state parent LTG NCA over the last 12-hour",
                  "dosing interval, by covariate arm. Typical-value",
                  "simulation; heterogeneity comes only from the",
                  "sampled body-weight / CLcr distribution within each arm.")
)
```

```{r pknca-metabolite}
pk_conc_gluc <- sim_typ |>
  dplyr::filter(time >= last_dose_time,
                time <= last_dose_time + tau_h) |>
  dplyr::mutate(time = time - last_dose_time) |>
  dplyr::select(id, time, Cc_gluc, arm, WT, CRCL) |>
  dplyr::filter(!is.na(Cc_gluc)) |>
  dplyr::rename(Cc = Cc_gluc)

pk_conc_gluc <- dplyr::bind_rows(
  pk_conc_gluc,
  pk_conc_gluc |>
    dplyr::distinct(id, arm) |>
    dplyr::mutate(time = 0)
) |>
  dplyr::distinct(id, arm, time, .keep_all = TRUE) |>
  dplyr::arrange(id, arm, time) |>
  dplyr::mutate(Cc = dplyr::coalesce(Cc, 0))

conc_obj_gluc <- PKNCA::PKNCAconc(
  pk_conc_gluc, Cc ~ time | arm + id,
  concu = "mg/L", timeu = "hr"
)

nca_gluc <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_gluc, dose_obj_parent, intervals = intervals_ss)
)
knitr::kable(
  summary(nca_gluc),
  caption = paste("Steady-state metabolite LTG-glu NCA over the last",
                  "12-hour dosing interval, by covariate arm. Same",
                  "typical-value simulation as the parent table above.")
)
```

## Comparison against published typical CL values

Milosheska 2016 does not publish a per-arm NCA table, but Table 4 reports
the typical CL/F point estimate (2.40 L/h at the reference stratum) and
gives closed-form covariate multipliers in the footnote formula. The
chunk below computes the paper's typical CL prediction for each arm from
the Table 4 formula and compares it to the model's simulated Cavg-derived
CL (dose / (Cavg * tau)). Any material disagreement here points at either
a covariate-encoding bug in the model file or a bug in the simulation
event table, not at a paper-vs-model discrepancy.

```{r cl-comparison}
# Local trapezoidal integrator (avoids the extra `pracma` dependency).
pracma_trapz <- function(x, y) {
  o <- order(x)
  x <- x[o]; y <- y[o]
  sum(0.5 * (y[-1] + y[-length(y)]) * diff(x))
}

# Paper's typical CL at each arm's covariate values (Milosheska 2016
# Table 4 footnote formula).
paper_cl <- function(smoke = 0, eiaed = 0, ugt_inh = 0,
                     m161ct = 0, m161tt = 0,
                     s372ag = 0, s372gg = 0,
                     wt = 70, crcl = 110) {
  2.40 *
    (1 - 0.579 * ugt_inh) *
    (wt / 70)^0.938 *
    (1 + 0.340 * smoke) *
    (1 + 0.546 * eiaed) *
    (1 - 0.0358 * m161ct) *
    (1 - 0.204  * m161tt) *
    (1 + 0.194  * s372ag) *
    (1 + 1.17   * s372gg) *
    (1 + 0.00328 * (crcl - 110))
}

paper_cl_by_arm <- tibble::tibble(
  arm      = c("Reference", "+Inducer", "+Inhibitor", "+UGT2B7 372 GG"),
  paper_CL = c(
    paper_cl(),
    paper_cl(eiaed   = 1),
    paper_cl(ugt_inh = 1),
    paper_cl(s372gg  = 1)
  )
)

# Simulated Cavg-derived CL from the typical-value simulation, restricted
# to subjects whose WT and CRCL are within +/- 5% of the reference (so
# the check is against the covariate multiplier alone, not weight-scaled
# CL differences).
sim_cl_at_ref <- sim_typ |>
  dplyr::filter(time >= last_dose_time,
                time <= last_dose_time + tau_h) |>
  dplyr::mutate(tau_time = time - last_dose_time) |>
  dplyr::filter(abs(WT - 70) < 3.5,       # ~+/- 5% of 70 kg
                abs(CRCL - 110) < 5.5) |> # ~+/- 5% of 110 mL/min
  dplyr::group_by(arm, id) |>
  dplyr::summarise(
    Cav_sim = pracma_trapz(tau_time, Cc) / tau_h,
    .groups = "drop"
  ) |>
  dplyr::group_by(arm) |>
  dplyr::summarise(
    sim_CL = mean(amt_mg / (Cav_sim * tau_h)),
    n_ref  = dplyr::n(),
    .groups = "drop"
  )

cmp <- paper_cl_by_arm |>
  dplyr::left_join(sim_cl_at_ref, by = "arm") |>
  dplyr::mutate(pct_diff = round(100 * (sim_CL - paper_CL) / paper_CL, 1))

knitr::kable(
  cmp,
  digits = c(NA, 3, 3, 0, 1),
  col.names = c("Arm", "Paper CL (L/h)",
                "Simulated CL (L/h)", "Ref subjects", "% diff"),
  caption = paste("Simulated typical CL/F at the reference WT / CRCL",
                  "vs the Milosheska 2016 Table 4 footnote formula",
                  "prediction. Absolute % diff > 20% would indicate a",
                  "covariate-encoding or event-table bug.")
)
```

# Assumptions and deviations

- **Table 4 for the LTG side, prose for the LTG-glu side.**
  Milosheska 2016 fits two sequential models: a parent-only model
  (Table 4) and a joint parent + metabolite model (Results
  'Population pharmacokinetic analysis of the lamotrigine-N-2-glucuronide'
  section, prose only). The joint-model paper text explicitly states
  that the LTG-side parameters "were very similar to those obtained
  with the parent drug model" but the paper does not tabulate them
  separately. The packaged model uses Table 4 values for the LTG side
  and the prose-reported values for the LTG-glu side, consistent with
  the paper's headline "parent + metabolite" narrative.
- **Reference for the metabolite CLcr effect.**
  Milosheska 2016 states the metabolite CLcr effect as "a decrease of
  0.759% per 1 mL/min decrease in CLcr" without explicitly re-declaring
  the reference CLcr. The packaged model uses the same reference as
  the parent equation (110 mL/min) because Milosheska 2016 Discussion
  paragraph on renal function explicitly uses that same 110 mL/min
  anchor when discussing metabolite CL: "The influence of renal
  function on CL_LTG-glu is more pronounced (a decrease by 0.759% per
  1 mL deviation from a standard CLcr of 110 mL/min)".
- **Complete parent-to-metabolite conversion.**
  The model equations follow the paper's stated assumption that
  lamotrigine is completely converted to N-2-glucuronide (Results
  paragraph 1 of the LTG-glu section: "We assumed that LTG was
  completely converted to LTG-glu ... renal elimination of unchanged
  LTG is insignificant ... conversion to LTG-N-5-glucuronide is a
  minor route of LTG metabolism"). Mass balance is written in the
  same mass units on both sides of the parent -> metabolite arrow,
  without an explicit molecular-weight (256 -> 430 g/mol) correction,
  following the Abduljalil 2009 clarithromycin + 14-OH-CLA precedent
  in nlmixr2lib. The reported `CL_LTG-glu / F_metab = 3.16 L/h` and
  `V_LTG-glu / F_metab = 110 L` are therefore "apparent" values that
  absorb the MW ratio and any deviation of the true metabolic
  conversion fraction from 1.
- **No IIV on metabolite volume.**
  Milosheska 2016 explicitly reports (Results LTG-glu paragraph 1):
  "Due to sparse concentration measurements and no prior data on
  pharmacokinetics of the metabolite we were not able to estimate the
  IIV of V_LTG-glu." The packaged model has no `etalvc_gluc` term.
- **CV%-to-omega^2 conversion for IIV.**
  Milosheska 2016 Table 4 reports IIV as CV% (Ka 71.1%, CL 33.1%,
  V 30.1%, LTG-glu CL 41.7%). The packaged `ini()` translates these
  to the internal log-normal variance via omega^2 = log(1 + CV^2)
  so the simulated between-subject geometric CV reproduces the
  reported values. The corresponding numeric omega^2 values
  (0.4090 / 0.1040 / 0.0867 / 0.1603) fall inside the bootstrap
  95% CIs the paper reports on the omega^2 scale (Table 4 columns
  '95% CI' for the IIV rows).
- **Cockcroft-Gault CLcr, not BSA-normalised eGFR.**
  Milosheska 2016 Methods 'Patients and blood sampling' uses raw
  Cockcroft-Gault CLcr (mL/min) throughout. The canonical `CRCL`
  register entry accepts raw Cockcroft-Gault CLcr as a documented
  variant (Delattre 2010 amikacin precedent). The vignette samples
  CRCL from a truncated normal centered at 110 with SD 30 mL/min to
  approximate the paper's Table 1 range 40.84 - 246.
- **Reference stratum indicators.**
  The paper reports its UGT2B7 -161C>T and 372A>G effects as two
  non-reference indicators each (CT vs CC, TT vs CC; AG vs AA, GG
  vs AA). The packaged `covariateData` lists only the non-reference
  indicators (`UGT2B7_M161CT` + `UGT2B7_M161TT` and `UGT2B7_372AG`
  + `UGT2B7_372GG`); the reference stratum (`UGT2B7_M161CC` and
  `UGT2B7_372AA`) is implicit when the paired non-reference
  indicators are both 0, matching the canonical `UGT2B7_M161CC` /
  `UGT2B7_372AA` register entries. `UGT2B7_M161CC` and
  `UGT2B7_372AA` remain registered as covariate canonicals for
  potential future reuse.
- **Prior-based Ka and V stabilisation.**
  Milosheska 2016 Methods 'Meta-analysis of previous population
  pharmacokinetic studies with LTG' incorporates prior information
  from eleven previously published lamotrigine popPK studies via
  NONMEM PRIOR functionality (NWPRI / INFV / INFTHETA) to stabilise
  Ka and V estimation. The paper's Table 4 point estimates are
  posterior-mode values under those informative priors. The packaged
  model treats them as the paper's final point estimates; users
  running a re-fit will need to supply the priors themselves.
- **Screened-but-not-retained covariates.**
  Age, sex (SEXF), Devine ideal body weight (IBW), aspartate and
  alanine transaminases (AST, ALT), UGT1A4 70C>A, SLC22A1 1222G>A,
  ABCB1 2677G>T/A, ABCB1 3435C>T, ABCB1 1236C>T, ABCB1 T-T-T
  haplotype, oral contraceptives, and LTG daily dose were tested
  but not retained in the final model. They are documented in the
  packaged `covariatesDataExcluded` list for provenance.
- **Errata.**
  No erratum or corrigendum to Milosheska 2016 was located on disk
  for this extraction. A search of the Wiley / BJCP corrections feed
  for "Milosheska 2016 erratum" and "bcp.12984 erratum" returned no
  hits (searched 2026-06-20); operators should reconfirm against
  the journal's current corrections listing if a re-extraction is
  undertaken.
- **Race / ethnicity.**
  The study cohort was 100% Slovenian adult epilepsy patients
  recruited from a single Ljubljana centre. Race and ethnicity were
  not tested as covariates. The packaged `population$race_ethnicity`
  is recorded as `c(Slovenian = 100)` and the virtual cohort does
  not stratify on race.
