---
title: "Sumatriptan (Lee 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sumatriptan (Lee 2015)}
  %\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 rxode2::rxode(readModelDb("Lee_2015_sumatriptan"))$reference`
* Description: `r rxode2::rxode(readModelDb("Lee_2015_sumatriptan"))$description`
* Article: <https://doi.org/10.12793/tcp.2015.23.2.66>

# Population

Lee 2015 fits a population PK model to plasma sumatriptan concentrations
collected during the reference-formulation arm of a single-center,
randomized, open-label, two-period, single-dose crossover bioequivalence
study at Kyungpook National University Hospital, Daegu, Korea. Twenty-six
healthy adult Korean males between 22 and 28 years (mean 23.9), body
weight 51-84 kg (mean 66.7), height 166.2-184.8 cm (mean 174.2), received
a single 50 mg oral dose of sumatriptan succinate (50 mg expressed as
sumatriptan free base; Table 1 of the source paper). 364 plasma
concentrations were collected on a dense 0-12 h grid (predose, 0.25, 0.5,
0.75, 1, 1.5, 2, 2.5, 3, 4, 6, 8, 10, 12 h after dose) and informed the
analysis.

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

```{r population-metadata}
str(rxode2::rxode(readModelDb("Lee_2015_sumatriptan"))$population)
```

# Source trace

Per-parameter origin is recorded inline next to each `ini()` entry in
`inst/modeldb/specificDrugs/Lee_2015_sumatriptan.R`. Consolidated here for
review:

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (CL/F)                                            | log(418) = 6.035        | Table 3 CL/F = 418 L/h (RSE 4%) |
| `lvc` (V/F)                                             | log(56.9) = 4.041       | Table 3 V/F = 56.9 L (RSE 35.4%) |
| `lka1` (first-order absorption rate)                    | log(0.62) = -0.478      | Table 3 ka1 = 0.62 1/h (RSE 9.13%) |
| `lka2` (transit -> central absorption rate)             | log(0.29) = -1.238      | Table 3 ka2 = 0.29 1/h (RSE 6.89%) |
| `lmtt` (mean transit time)                              | log(1.94) = 0.663       | Table 3 MTT = 1.94 h (RSE 9.89%) |
| `ln` (number of transit compartments)                   | log(11) = 2.398         | Table 3 n = 11 (RSE 23.2%); estimated as a continuous value via the analytical input form |
| `lalag1` (first-order absorption lag)                   | log(0.24) = -1.427      | Table 3 ALAG1 = 0.24 h (RSE 1.32%) |
| `lfr` (fraction routed via transit chain)               | log(0.56) = -0.580      | Table 3 f = 0.56 (RSE 6.18%); the remaining 1 - f = 0.44 is routed via the first-order arm |
| `etalcl` (BSV CL, variance)                             | 0.0337                  | Table 3 BSV(CL/F) = 18.5% CV (omega^2 = log(1 + 0.185^2)) |
| `etalvc` (BSV V, variance)                              | 0.4072                  | Table 3 BSV(V/F) = 70.9% CV (omega^2 = log(1 + 0.709^2)) |
| `etalka2` (BSV ka2, variance)                           | 0.0587                  | Table 3 BSV(ka2) = 24.6% CV (omega^2 = log(1 + 0.246^2)) |
| `etalmtt` (BSV MTT, variance)                           | 0.1193                  | Table 3 BSV(MTT) = 35.6% CV (omega^2 = log(1 + 0.356^2)) |
| `etalfr` (BSV f, variance)                              | 0.0205                  | Table 3 BSV(f) = 14.4% CV (omega^2 = log(1 + 0.144^2)) |
| `propSd` (proportional residual SD)                     | 0.21                    | Table 3 proportional error = 0.21 (RSE 7.3%) |
| `addSd` (additive residual SD, ng/mL)                   | 0.30                    | Table 3 additive error = 0.30 ng/mL (RSE 15.2%) |
| One-compartment central + two parallel absorption arms  | n/a                     | Methods page 67-68 ("base PK model" section), Figure 2 schematic, and supplement NONMEM code page 74 ($MODEL block: DEPOT1, DEPOT2, CENTRAL) |
| Savic 2007 analytical transit-compartment input rate    | n/a                     | Methods page 68 (transit compartment equations), supplement $DES block page 74 (analytical kernel using Stirling's approximation); replaced here by rxode2's `transit(n, mtt, fr)` which uses `lgamma(n+1)` |
| Combined additive + proportional residual error         | n/a                     | Supplement $ERROR block page 74: W = sqrt(THETA(9)^2 + THETA(10)^2 * IPRED^2); Y = IPRED + W * ERR(1) |

# Virtual cohort

Original individual concentration-time data are not publicly available;
the figures below use a virtual population whose covariate distributions
match the published trial (Table 1): 26 healthy Korean adult males, ages
22-28 years, weights 51-84 kg, receiving a single 50 mg oral dose of
sumatriptan.

The Lee 2015 model has two parallel absorption arms; each administered
dose therefore generates **two** dose records (one to `depot` for the
first-order arm and one to `depot2` for the transit-chain arm). Both
records carry the same `amt`; the model performs the `(1 - fr) / fr`
split internally via bioavailability settings. This is a constraint of
rxode2's compartment-aware `transit()` function (see the
implementation notes in the model file).

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

n_sub <- 26L

obs_grid <- c(0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4, 6, 8, 10, 12)

make_cohort <- function(n, label, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  base <- tibble::tibble(id = ids, cohort = label)
  expand_one <- function(row) {
    dose_depot  <- data.frame(id = row$id, cohort = row$cohort, time = 0,
                              evid = 1L, amt = 50, cmt = "depot")
    dose_depot2 <- data.frame(id = row$id, cohort = row$cohort, time = 0,
                              evid = 1L, amt = 50, cmt = "depot2")
    obs <- data.frame(id = row$id, cohort = row$cohort, time = obs_grid,
                      evid = 0L, amt = 0, cmt = "Cc")
    rbind(dose_depot, dose_depot2, obs)
  }
  rows <- split(base, base$id)
  do.call(rbind, lapply(rows, expand_one))
}

events <- make_cohort(n_sub, "50 mg oral, single dose")
# Multiple dose events per id at the same time are intentional (two
# parallel absorption arms); only the (id, time, evid, cmt) combination
# needs to be unique, which it is.
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid", "cmt")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Lee_2015_sumatriptan")
sim <- rxode2::rxSolve(mod(), events = events, keep = c("cohort"))
sim <- as.data.frame(sim)
```

# Replicate published figures

## Figure 1 -- individual concentration profiles with multiple peaks

Lee 2015 Figure 1 shows individual plasma concentration vs time plots
for the 26 subjects, with the bold red line marking the median. Each
profile shows the characteristic multiple-peaks pattern that motivated
the combined transit-compartment + first-order absorption model. The
plot below renders the simulated individual profiles in light gray and
overlays the population median in red.

```{r figure-1-individual-profiles}
sim_obs <- sim |> dplyr::filter(time > 0)

median_sim <- sim_obs |>
  dplyr::group_by(time) |>
  dplyr::summarise(median_Cc = median(Cc, na.rm = TRUE), .groups = "drop")

ggplot(sim_obs, aes(time, Cc, group = id)) +
  geom_line(alpha = 0.25, colour = "gray40", linewidth = 0.3) +
  geom_line(data = median_sim, aes(time, median_Cc, group = NULL),
            colour = "red", linewidth = 1.0) +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 6, 8, 10, 12)) +
  labs(x = "Time after dose (h)", y = "Sumatriptan plasma conc (ng/mL)",
       title = "Individual concentration profiles, 50 mg oral sumatriptan",
       subtitle = "Replicates Lee 2015 Figure 1 (bold red line = median)") +
  theme_minimal()
```

## Figure 5 -- representative individual showing multiple peaks

Lee 2015 Figure 5 displays the multiple-peaks behavior in a single
representative subject (observed circles vs individual model
prediction). We pick a representative simulated subject (using a high
eta on the transit fraction to amplify the multi-peak shape) and plot
its typical-individual profile.

```{r figure-5-representative}
# Simulate a typical-value (no-IIV) profile to visualise the structural
# multi-peak shape clearly. Two dose records (depot + depot2) per
# administration as required by the model's parallel-absorption layout.
typ_grid <- seq(0.05, 12, by = 0.05)
typ_events <- rbind(
  data.frame(id = 1L, time = 0, evid = 1L, amt = 50, cmt = "depot"),
  data.frame(id = 1L, time = 0, evid = 1L, amt = 50, cmt = "depot2"),
  data.frame(id = 1L, time = typ_grid, evid = 0L, amt = 0, cmt = "Cc")
)

mod_typ <- rxode2::zeroRe(mod())
typ_sim <- as.data.frame(rxode2::rxSolve(mod_typ, events = typ_events))

ggplot(typ_sim |> dplyr::filter(time > 0),
       aes(time, Cc)) +
  geom_line(colour = "red", linewidth = 0.8) +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 6, 8, 10, 12)) +
  labs(x = "Time after dose (h)", y = "Sumatriptan plasma conc (ng/mL)",
       title = "Typical-value profile, 50 mg oral sumatriptan",
       subtitle = "Replicates the multi-peak shape of Lee 2015 Figure 5 (representative individual prediction)") +
  theme_minimal()
```

## Figure 4 -- visual predictive check 5th / median / 95th

Lee 2015 Figure 4 is a visual predictive check from 1,000 simulated
datasets, showing the 90% prediction interval (gray band) and observed
5th / 50th / 95th percentiles. The simulated VPC below uses the
packaged model to reproduce that figure's shape.

```{r figure-4-vpc}
vpc_grid <- seq(0.1, 12, by = 0.1)
vpc_events <- dplyr::bind_rows(lapply(seq_len(1000), function(i) {
  rbind(
    data.frame(id = i, time = 0, evid = 1L, amt = 50, cmt = "depot"),
    data.frame(id = i, time = 0, evid = 1L, amt = 50, cmt = "depot2"),
    data.frame(id = i, time = vpc_grid, evid = 0L, amt = 0, cmt = "Cc")
  )
}))
vpc_sim <- as.data.frame(rxode2::rxSolve(mod(), events = vpc_events))

vpc_summary <- vpc_sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    p05 = quantile(Cc, 0.05, na.rm = TRUE),
    p50 = quantile(Cc, 0.50, na.rm = TRUE),
    p95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_summary, aes(time, p50)) +
  geom_ribbon(aes(ymin = p05, ymax = p95), fill = "gray70", alpha = 0.6) +
  geom_line(colour = "red", linewidth = 0.8) +
  scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 6, 8, 10, 12)) +
  labs(x = "Time after dose (h)", y = "Sumatriptan plasma conc (ng/mL)",
       title = "Visual predictive check, 50 mg oral sumatriptan (1,000 simulated subjects)",
       subtitle = "Replicates Lee 2015 Figure 4 (90% prediction interval gray; median red)") +
  theme_minimal()
```

# PKNCA validation

Compute Cmax, Tmax, AUClast, AUCinf, and terminal half-life on the
26-subject simulated cohort and compare against published literature
values for 50 mg oral sumatriptan (the source paper does not tabulate
NCA in Table 3 -- it reports population PK parameters only; published
literature values are taken from the Cosson et al. 1999 and Christensen
et al. 2004 popPK comparisons that Lee 2015 cites).

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort)

dose_df <- events |>
  dplyr::filter(evid == 1, cmt == "depot") |>     # one row per user-facing 50 mg dose (the depot record stands in for the full administration)
  dplyr::select(id, time, amt, cohort)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | cohort + id,
                             concu = "ng/mL", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | cohort + id,
                             doseu = "mg")

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

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

nca_summary <- summary(nca_res)
knitr::kable(nca_summary,
             caption = "Simulated NCA parameters, 50 mg oral sumatriptan (n = 26 virtual subjects).")
```

## Comparison against published reference values

Lee 2015 does not publish a per-subject NCA table; the paper validates
its model via VPC (Figure 4) and goodness-of-fit (Figure 3). Two
internally-consistent sanity checks against the published popPK
literature are available:

* **Apparent oral clearance**. Lee 2015 reports CL/F = 418 L/h. The
  typical-value model produces an AUCinf consistent with
  dose / (CL/F) = 50 mg / 418 L/h ~= 0.120 mg.h/L = 120 ng.h/mL for the
  typical subject. The simulated median AUCinf in the NCA table above
  should fall within ~20 % of this analytical target.
* **Multi-peak signature**. The original popPK comparisons (Cosson et al.
  1999, Christensen et al. 2004) report Cmax in the 50-80 ng/mL range
  and Tmax in the 1.5-2.5 h window for 50 mg oral sumatriptan in
  healthy subjects, with the median profile peaking, dipping, and
  re-peaking around 4-6 h after dose. The packaged model's VPC (Figure
  4 above) should reproduce that double-peak shape.

```{r dose-cl-check}
analytical_aucinf_ng_h_per_mL <- 50 * 1000 / 418   # dose mg * 1000 / CL/F L/h = ng.h/mL
cat(sprintf("Analytical AUCinf at CL/F = 418 L/h: %.1f ng*h/mL\n",
            analytical_aucinf_ng_h_per_mL))
```

# Assumptions and deviations

* **Final published model is M6, not M5.** Lee 2015 selects "Model M6"
  (M5 base plus a creatinine clearance covariate on the transit
  fraction `f`, Table 2 row M6) as the final published model. The
  packaged model implements the **M5 base** structure (no covariate)
  for the following reason: the source paper reports the M6 covariate
  coefficient as `f, CrCL = -0.985` (Table 3, RSE 34.6 %) but does
  **not** state the functional form (power vs linear-additive vs
  exponential vs logit) or the reference CrCL value, and the
  supplement's NONMEM code (page 74) only contains the M5 control
  stream without the covariate. The cohort's CrCL distribution is also
  not tabulated (Table 1 lists only age, sex, weight, height). Per the
  skill's "Covariate encoding ambiguous" stop-and-ask trigger and the
  operator's directive, M5 is the faithfully-reproducible structure
  from on-disk sources; downstream simulation for healthy adults with
  normal renal function (the published cohort) is essentially
  unaffected by the M5 -> M6 simplification. Users wishing to add the
  M6 effect can wrap the typical `fr` in their preferred functional
  form, e.g.
  `fr_i = exp(lfr + etalfr) * (CRCL / CRCL_ref)^(-0.985)`, and
  document the assumed reference. **Erratum/clarification candidate:**
  the operator confirmed via sidecar that the M6 covariate form is
  unrecoverable from on-disk sources; this is documented for any
  future erratum or author correspondence that resolves the form.

* **Analytical Savic transit input via lgamma rather than Stirling.**
  The supplement's NONMEM code (page 74) evaluates the gamma kernel
  log(N!) via the Stirling approximation
  `log(sqrt(2*pi)) + (N + 0.5) * log(N) - N`. The packaged model uses
  rxode2's built-in `transit(n, mtt, bio)` which evaluates the same
  kernel using `lgamma(n + 1)` directly. At Lee 2015's estimated
  N = 11, Stirling's approximation differs from the exact gamma by
  about 0.03 % in log(N!) -- well below simulation precision.

* **No covariate columns required.** The M5 base structure has no
  covariate dependencies, so the packaged model's `covariateData` is
  an empty list. The published cohort is a single ethnic and
  demographic group (healthy Korean adult males, ages 22-28); the
  paper notes that age, weight, and BMI did not appear to affect any
  PK parameter, and CrCL was the only covariate selected in Lee
  2015's covariate search.

* **Two-depot dose routing.** Lee 2015's NONMEM control stream
  routes (1 - f) * AMT into DEPOT1 (the first-order arm) via
  F1 = 1 - FR, and the remaining f * AMT enters DEPOT2 via the
  analytical Savic input rate while F2 = 0 suppresses any bolus
  deposition into DEPOT2. The packaged rxode2 model mirrors this
  exactly: `f(depot) <- 1 - fr` (the first-order arm's bioavailability
  in rxode2 syntax), `f(depot2) <- 0` (suppressed bolus), and the
  analytical input rate is supplied through `transit(n, mtt, fr)`
  in the `d/dt(depot2)` line. `transit()` reads the raw parent dose
  amount via `podo()` and time-after-dose via `tad()`, so multi-dose
  simulation remains coherent in principle, although the source study
  is single-dose only.

* **Lag time applies only to the first-order arm.** The supplement
  declares ALAG1 only on COMP=DEPOT1; the transit arm has no lag
  time (the analytical input rate begins at t = 0). The packaged
  model preserves this via `alag(depot) <- alag1` (with no `alag()`
  call on `depot2`).
