---
title: "Pregabalin (Shoji 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pregabalin (Shoji 2011)}
  %\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

```{r meta}
mod <- readModelDb("Shoji_2011_pregabalin")
mod_meta <- nlmixr2est::nlmixr(mod)$meta
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article (DOI): <https://doi.org/10.1111/j.1365-2125.2011.03932.x>

This vignette validates the packaged `Shoji_2011_pregabalin` model -- a
one-compartment oral PK model for pregabalin with first-order absorption,
absorption lag-time, and creatinine-clearance-driven elimination, fit to
data from 14 clinical trials in 616 subjects (195 healthy volunteers,
267 post-herpetic neuralgia patients, 154 painful diabetic peripheral
neuropathy patients). Validation focuses on Table 4 of the source paper,
which reports mean steady-state AUC by dose and renal-function stratum
in PT05 (the diabetic-peripheral-neuropathy efficacy study), and on the
deterministic CL/F vs CLcr relationship reproduced in the paper's Figure 1.

# Population

The Shoji 2011 analysis pooled data from 14 pregabalin clinical trials
(Table 1). Nine phase-1 studies (HV01-HV09) enrolled 195 healthy adults
with dense PK sampling, including substudies of subjects with impaired
renal function (HV05) and elderly subjects (HV07). Four phase-3 trials
in post-herpetic neuralgia (PT01-PT04, n=267) and one in painful
diabetic peripheral neuropathy (PT05, n=154) used sparse outpatient
sampling (1-2 samples per patient). Mean (range) age was 59.1 (19-101)
years, mean total body weight 71.0 (31-142) kg, and mean Cockcroft-Gault
CrCl 86.0 (10.0-230) mL/min (Table 1, All Total row). The cohort was
37.2% female (Table 2). Ethnicity composition: 52.9% White, 41.6% Asian
(including 252 Japanese subjects), 4.55% Other, 0.97% Black. Pregabalin
doses ranged from a single 1 mg dose to 600 mg/day BID or 900 mg/day
TID in phase-1 studies; phase-3 doses were 75-300 mg BID or 25-200 mg
TID up to 600 mg/day.

The same information is available programmatically via the model's
`population` metadata:

```{r population}
str(mod_meta$population)
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in `inst/modeldb/specificDrugs/Shoji_2011_pregabalin.R`.
The table below collects them in one place; all final-model values come
from Shoji 2011 Table 3 final-model column unless otherwise noted.

| Parameter / equation                                  | Value                | Source location |
|-------------------------------------------------------|----------------------|-----------------|
| `lcl` (typical CL/F at reference CRCL=86 mL/min)      | `log(0.0462 * 86)`   | Table 3 q_CL/F final model (slope 0.0462 L/h per mL/min) x Table 1 mean CRCL 86 mL/min |
| `th_bp` (CRCL break point above which CL/F saturates) | 107 mL/min           | Table 3 q_BP final model |
| `e_ibw_cl` (IBW power exponent on CL/F)               | 0.354                | Table 3 q_IBW on CL/F final model |
| `lvc` (typical V/F at reference covariates)           | `log(35.6)`          | Table 3 q_V/F final model |
| `e_ibw_vc` (IBW power exponent on V/F)                | 0.819                | Table 3 q_IBW on V/F final model |
| `e_bmi_vc` (BMI power exponent on V/F)                | 0.525                | Table 3 q_BMI on V/F final model |
| `e_age_vc` (AGE power exponent on V/F)                | -0.125               | Table 3 q_Age on V/F final model |
| `e_sexf_vc` (female-vs-male V/F ratio)                | 0.906                | Table 3 q_Gender on V/F final model |
| `lka` (ka at FED = 0)                                 | `log(7.99)`          | Table 3 q_ka final model |
| `e_food_ka` (food effect on ka)                       | -0.930               | Table 3 q_Food on ka final model |
| `ltlag` (tlag at FED = 0)                             | `log(0.243)`         | Table 3 q_t_lag final model |
| `e_food_tlag` (food effect on tlag)                   | 0.811                | Table 3 q_Food on tlag final model |
| `etalcl ~ 0.02168`                                    | `log(1 + 0.148^2)`   | Table 3 CV%(CL/F) final model 14.8% |
| `etalvc ~ 0.00727`                                    | `log(1 + 0.0854^2)`  | Table 3 CV%(V/F) final model 8.54% |
| `etalka ~ 0.62526`                                    | `log(1 + 0.932^2)`   | Table 3 CV%(ka) final model 93.2% |
| `propSdHealthy <- 0.220`                              | 0.220                | Table 3 CV%(healthy) final model 22.0% |
| `addSdHealthy  <- 0.0239`                             | 0.0239 ug/mL         | Table 3 SD(healthy) final model |
| `propSdPatient <- 0.285`                              | 0.285                | Table 3 CV%(patient) final model 28.5% |
| `addSdPatient  <- 0.236`                              | 0.236 ug/mL          | Table 3 SD(patient) final model |
| CL/F = `exp(lcl) * (min(CRCL, th_bp) / 86) * (IBW / 62)^e_ibw_cl * exp(etalcl)` | n/a | Final-model equation, reformulated from the paper's slope form (see model file header) |
| V/F = `exp(lvc) * (IBW/62)^e_ibw_vc * (BMI/25)^e_bmi_vc * (AGE/59)^e_age_vc * e_sexf_vc^SEXF * exp(etalvc)` | n/a | Final-model equation, Results "Final model development" |
| ka = `exp(lka + etalka) * (1 + e_food_ka * FED)`      | n/a                  | Final-model equation, paper's NONMEM (1 + theta * FED) form |
| tlag = `exp(ltlag) * (1 + e_food_tlag * FED)`         | n/a                  | Same form on lag-time |
| `Cc ~ prop(propSd) + add(addSd)`                      | n/a                  | Methods paragraph 3: Y_ij = C_ij + C_ij*eps1 + eps2 |
| Residual error stratification by `DIS_HEALTHY`        | n/a                  | Table 3 separate (healthy)/(patient) residual rows |

# Virtual cohort

The original observed pregabalin concentrations are not publicly
available. The virtual cohort below mirrors the three groups in
Shoji 2011 Table 4 -- PT05 diabetic-peripheral-neuropathy patients
allocated to (i) 300 mg BID with CrCl >= 60 mL/min (n=31), (ii)
150 mg BID with CrCl >= 60 mL/min (n=109), and (iii) 150 mg BID with
30 <= CrCl < 60 mL/min (n=14) -- so that simulated steady-state AUC
can be compared 1:1 against the paper's reported values. Covariate
distributions match the PT05 phase-3 demographics summary in
Shoji 2011 Table 1 (mean age 60.9 years; mean total body weight
65.6 kg; mean CrCl 99.3 mL/min; ranges in the helpers below).

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

# Helper: per-cohort tibble of covariates. id_offset shifts subject IDs
# so multi-cohort bind_rows() produces disjoint id ranges (rxSolve treats
# id as the subject key; duplicate ids across cohorts silently collapse).
make_cohort_covs <- function(n, dose_mg, treatment,
                             crcl_mean, crcl_sd, crcl_range,
                             id_offset) {
  # Truncated log-normal CRCL so the cohort sits in its stratum (>=60 or
  # 30-60); SD is small so most subjects stay inside the band.
  crcl <- exp(rnorm(n, mean = log(crcl_mean), sd = crcl_sd))
  crcl <- pmin(pmax(crcl, crcl_range[1]), crcl_range[2])

  # Total body weight: log-normal around PT05 mean 65.6 kg, range 31-113.
  wt <- exp(rnorm(n, mean = log(65.6), sd = log(113 / 31) / 4))
  wt <- pmin(pmax(wt, 31), 113)

  # Age: normal around PT05 mean 60.9 y, range 35-85.
  age <- pmin(pmax(rnorm(n, mean = 60.9, sd = 10), 35), 85)

  # Sex: PT05-typical breakdown roughly 50/50 (paper-summary value not
  # reported per cohort; use overall 37.2% female from Table 2).
  sexf <- rbinom(n, 1, prob = 0.372)

  # Height: log-normal around 167 cm (population-mean imputation in
  # Methods, Demographic data); SD chosen to land BMI near 25.
  ht_cm <- exp(rnorm(n, mean = log(167), sd = 0.06))

  # IBW (Devine-family): men 50 + 0.91 * (ht_cm - 152.4);
  # women 45.5 + 0.91 * (ht_cm - 152.4). The Methods reference [12]
  # is the Devine variant; we use the metric-units form here.
  ibw_men   <- 50.0 + 0.91 * (ht_cm - 152.4)
  ibw_women <- 45.5 + 0.91 * (ht_cm - 152.4)
  ibw <- ifelse(sexf == 1, ibw_women, ibw_men)
  ibw <- pmax(ibw, 30)  # guard for very short stature

  # BMI (kg/m^2)
  bmi <- wt / (ht_cm / 100)^2

  tibble::tibble(
    id          = id_offset + seq_len(n),
    treatment   = treatment,
    dose_mg     = dose_mg,
    CRCL        = crcl,
    WT          = wt,
    AGE         = age,
    SEXF        = sexf,
    HT          = ht_cm,
    IBW         = ibw,
    BMI         = bmi,
    FED         = 0L,           # PT05 protocol does not stratify by food
    DIS_HEALTHY = 0L            # PT05 patients are not healthy volunteers
  )
}

# Three cohorts matching Shoji 2011 Table 4.
covs_a <- make_cohort_covs(
  n = 31L, dose_mg = 300, treatment = "300 mg BID, CrCl>=60",
  crcl_mean = 95, crcl_sd = 0.20, crcl_range = c(60, 230),
  id_offset = 0L
)
covs_b <- make_cohort_covs(
  n = 109L, dose_mg = 150, treatment = "150 mg BID, CrCl>=60",
  crcl_mean = 95, crcl_sd = 0.20, crcl_range = c(60, 230),
  id_offset = 1000L
)
covs_c <- make_cohort_covs(
  n = 14L, dose_mg = 150, treatment = "150 mg BID, 30<=CrCl<60",
  crcl_mean = 45, crcl_sd = 0.10, crcl_range = c(30, 60),
  id_offset = 2000L
)
covs_all <- bind_rows(covs_a, covs_b, covs_c)
stopifnot(!anyDuplicated(covs_all$id))

# Build event table per subject: BID dosing for 20 doses (10 days) so
# the system is at steady-state in the final dosing interval. Sample
# the final tau-interval [228, 240] every 0.25 h.
make_events <- function(cov_row) {
  amt  <- cov_row$dose_mg
  dose_times <- seq(0, 228, by = 12)        # 20 BID doses, last at 228 h
  sample_times <- c(0, seq(228, 240, by = 0.25))

  doses <- tibble::tibble(
    id   = cov_row$id,
    time = dose_times,
    evid = 1L,
    amt  = amt,
    cmt  = "depot",
    dv   = NA_real_
  )
  obs <- tibble::tibble(
    id   = cov_row$id,
    time = sample_times,
    evid = 0L,
    amt  = NA_real_,
    cmt  = NA_character_,
    dv   = NA_real_
  )
  bind_rows(doses, obs) |>
    mutate(
      treatment   = cov_row$treatment,
      CRCL        = cov_row$CRCL,
      WT          = cov_row$WT,
      AGE         = cov_row$AGE,
      SEXF        = cov_row$SEXF,
      HT          = cov_row$HT,
      IBW         = cov_row$IBW,
      BMI         = cov_row$BMI,
      FED         = cov_row$FED,
      DIS_HEALTHY = cov_row$DIS_HEALTHY
    ) |>
    arrange(time, desc(evid))
}

events <- bind_rows(lapply(seq_len(nrow(covs_all)), function(i) {
  make_events(covs_all[i, ])
}))

# Cheap regression guard for the disjoint-id property; see
# references/vignette-template.md Notes "Multi-cohort simulations".
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
sim_stoch <- rxode2::rxSolve(
  object = mod, events = events,
  keep   = c("treatment", "CRCL", "IBW", "BMI", "AGE", "SEXF", "FED", "DIS_HEALTHY")
) |>
  as.data.frame()

mod_typical <- rxode2::zeroRe(mod)
sim_typical <- rxode2::rxSolve(
  object = mod_typical, events = events,
  keep   = c("treatment", "CRCL", "IBW", "BMI", "AGE", "SEXF")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 1 -- CL/F vs CrCl

Shoji 2011 Figure 1 plots individual CL/F vs CrCl from the final model
and overlays the typical-value relationship CL/F = `q_CL/F * min(CrCl, th_bp)`.
The deterministic typical-value relationship is reproduced below across
the full CrCl range observed in the pooled cohort (10-230 mL/min).

```{r figure-1}
# Replicates Figure 1 of Shoji 2011 (right panel): typical-value CL/F
# as a strict linear function of min(CrCl, th_bp), with the slope and
# break point taken directly from the model file.
crcl_grid <- seq(10, 230, by = 1)
clf_typ <- 0.0462 * pmin(crcl_grid, 107)  # at IBW = 62, by definition

# Empirical-Bayes-style scatter: per-subject CL/F from the simulated
# stochastic cohort, displayed as a Cl/F vs CRCL scatter to mirror the
# paper's Figure 1 visual.
clf_per_subject <- sim_stoch |>
  group_by(id) |>
  summarise(
    CRCL    = first(CRCL),
    IBW     = first(IBW),
    cl_indiv = first(CRCL) * 0.0462 *
                 pmin(first(CRCL), 107) / first(CRCL) *
                 (first(IBW) / 62)^0.354,
    .groups = "drop"
  )

ggplot() +
  geom_point(data = clf_per_subject, aes(CRCL, cl_indiv),
             alpha = 0.4, colour = "black", size = 1.1) +
  geom_line(data = tibble(crcl_grid = crcl_grid, clf_typ = clf_typ),
            aes(crcl_grid, clf_typ),
            colour = "red", linewidth = 0.9) +
  labs(
    x = "CrCl (mL/min)",
    y = "CL/F (L/h)",
    title    = "Figure 1 -- CL/F vs CrCl",
    subtitle = paste0("Black: per-subject typical CL/F across the virtual cohort. ",
                      "Red: typical-value line CL/F = 0.0462 * min(CrCl, 107).")
  ) +
  theme_minimal()
```

## Figure 4 -- VPC by dose and renal-function stratum

Shoji 2011 Figure 4 shows the visual predictive check of pregabalin
concentration vs time after the most recent administration, stratified
by dose and renal function. The reproduction below shows the
steady-state final tau-interval (228-240 h) by treatment group.

```{r figure-4}
# Replicates Figure 4 of Shoji 2011: predicted percentiles per dose-by-
# renal-function group. Time is mapped to time-after-most-recent-dose
# (shift t -> t - 228 so the x-axis runs 0-12 h).
sim_panel <- sim_stoch |>
  filter(time >= 228) |>
  mutate(t_tad = time - 228)

sim_panel |>
  group_by(treatment, t_tad) |>
  summarise(
    Q025 = quantile(Cc, 0.025, na.rm = TRUE),
    Q50  = quantile(Cc, 0.50,  na.rm = TRUE),
    Q975 = quantile(Cc, 0.975, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(t_tad, Q50)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975),
              fill = "gray70", alpha = 0.5) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~treatment, ncol = 3) +
  labs(
    x = "Time after most recent dose (h)",
    y = "Pregabalin plasma concentration (ug/mL)",
    title    = "Figure 4 -- VPC by treatment group at steady state",
    subtitle = paste0("Median (line) and 95% prediction interval (ribbon) ",
                      "over the final tau interval (h 228-240).")
  ) +
  theme_minimal()
```

# PKNCA validation

PKNCA is used to compute AUC over the steady-state tau interval
(228-240 h) for each treatment group. The published comparator is
Shoji 2011 Table 4, which reports mean AUC(tau,ss) in the same three
PT05 sub-groups.

```{r pknca}
# Filter to the SS tau interval and shift t -> t - 228 so PKNCA's
# interval [0, 12] aligns with the final dosing interval. Time-zero
# row is guaranteed by the simulation grid (sample includes t = 228 = 0
# in shifted coordinates).
sim_for_nca <- sim_stoch |>
  dplyr::filter(time >= 228, !is.na(Cc)) |>
  dplyr::mutate(t_tad = time - 228) |>
  dplyr::select(id, time = t_tad, Cc, treatment) |>
  as.data.frame()

# Defensive time-zero ensure (pknca-recipes.md "Time-zero records").
sim_for_nca <- dplyr::bind_rows(
  sim_for_nca,
  sim_for_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

# Dose record at t = 0 in the shifted frame; one row per subject for
# the final dose carried into the SS tau interval.
doses_for_nca <- covs_all |>
  dplyr::transmute(
    id, time = 0, amt = dose_mg, treatment
  ) |>
  as.data.frame()

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ time | treatment + id,
  concu   = "ug/mL",
  timeu   = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  data    = doses_for_nca,
  formula = amt ~ time | treatment + id,
  doseu   = "mg"
)

# Single SS tau interval [0, 12]; cmax + tmax + cmin + AUC over tau.
intervals <- data.frame(
  start    = 0,
  end      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  auclast  = TRUE,
  cav      = TRUE
)
nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)
```

## Comparison against Shoji 2011 Table 4

```{r nca_compare}
published <- tibble::tribble(
  ~treatment,                       ~cmax,  ~auclast,
  "300 mg BID, CrCl>=60",           NA,     75.5,
  "150 mg BID, CrCl>=60",           NA,     37.5,
  "150 mg BID, 30<=CrCl<60",        NA,     80.3
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = published,
  by        = "treatment",
  units     = c(auclast = "ug*h/mL", cmax = "ug/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste0("Simulated vs. Shoji 2011 Table 4 mean AUC(tau,ss). ",
                   "* differs from reference by >20%."),
  align   = c("l", "l", "r", "r", "r")
)
```

# Assumptions and deviations

- **CL/F parameterisation rewritten in canonical ratio form.** Shoji 2011
  parameterises CL/F as a strict slope on CrCl with no intercept:
  `q_CL/F = 0.0462 L/h per mL/min`. To pair the `etalcl` random effect
  with a canonical `lcl` fixed-effect, the model file anchors CL/F at
  the population-mean CrCl 86 mL/min as `exp(lcl) * (min(CRCL, 107) / 86)`,
  with `exp(lcl) = 0.0462 * 86 = 3.97 L/h`. The two forms are
  mathematically identical -- the slope is recovered as
  `exp(lcl) / 86 = 0.0462`.
- **Reference covariate values for V/F.** The paper's Discussion
  explicitly cites IBW = 62 kg as the population-mean IBW and AGE =
  59 years as the population-mean age. The BMI reference value 25 kg/m^2
  is inferred as the rounded population-mean BMI (mean TBW 71 kg + mean
  HT 167 cm -> BMI 25.4); the paper's own BMI sensitivity comment ("V/F
  decreased to 84% when BMI changed from 25 to 18 kg/m^2") confirms 25
  as the reference comparator. At the reference covariates (IBW=62,
  BMI=25, AGE=59, SEXF=0) the typical V/F is `exp(lvc) = 35.6 L`,
  matching Table 3 q_V/F directly.
- **Diagonal IIV.** Methods paragraph 3 specifies independent
  correlation structure on the IIV variance-covariance matrix W. We
  encode this as diagonal `etalcl`, `etalvc`, `etalka` with no
  off-diagonal block.
- **Residual error stratification.** Table 3 reports separate residual
  CV / SD for healthy (CV 22.0%, SD 0.0239 ug/mL) and patient
  (CV 28.5%, SD 0.236 ug/mL) cohorts; the model file switches the
  combined proportional + additive error magnitudes at runtime via
  `DIS_HEALTHY`. The simulation in this vignette uses the patient
  values (`DIS_HEALTHY = 0`) because the Table 4 validation target is
  the PT05 patient cohort.
- **IBW formula in the virtual cohort.** The source paper references a
  Devine-family IBW formula [12] but does not transcribe the formula
  in the trimmed methods text. The cohort builder uses the metric-units
  Devine variant (men 50 + 0.91 * (HT_cm - 152.4); women 45.5 + 0.91 *
  (HT_cm - 152.4)) to derive IBW from a simulated HT_cm. Per-paper IBW
  is recommended when a user has access to the source dataset.
- **Food and disease-state covariates not stratified in this validation.**
  Table 4's PT05 outpatient sampling protocol did not stratify by food
  state; we therefore set FED = 0 (fasted) for the SS-AUC validation,
  which is unaffected because AUC over a complete tau interval depends
  only on CL/F, not on ka or tlag. DIS_HEALTHY is set to 0 (patient
  cohort) to select the patient residual-error magnitudes.
- **Bootstrap-based parameter uncertainty is not propagated.** The
  packaged model uses the point estimates from Shoji 2011 Table 3
  final-model column; the nonparametric bootstrap CIs (column 4 of
  Table 3) are documentary and not part of the simulation.
