Skip to contents

Sumatriptan (Lee 2015)

Source: vignettes/articles/Lee_2015_sumatriptan.Rmd
Lee_2015_sumatriptan.Rmd

Model and source

  • Citation: Lee J, Lim M, Seong SJ, Park S-M, Gwon M-R, Han S, Lee SM, Kim W, Yoon Y-R, Yoo H-D. Population pharmacokinetic analysis of the multiple peaks phenomenon in sumatriptan. Transl Clin Pharmacol. 2015;23(2):66-74. doi:10.12793/tcp.2015.23.2.66.
  • Description: One-compartment population PK model for oral sumatriptan in healthy Korean male volunteers (Lee 2015): two parallel absorption routes (first-order absorption with lag time, and a transit-compartment chain with the Savic 2007 analytical input form) into a single central compartment with linear elimination. Captures the multiple-peaks absorption phenomenon reported in oral sumatriptan.
  • 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:

str(rxode2::rxode(readModelDb("Lee_2015_sumatriptan"))$population)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> List of 13
#>  $ species       : chr "human"
#>  $ n_subjects    : int 26
#>  $ n_studies     : int 1
#>  $ age_range     : chr "22-28 years"
#>  $ age_median    : chr "23.9 years (mean)"
#>  $ weight_range  : chr "51-84 kg"
#>  $ weight_median : chr "66.7 kg (mean)"
#>  $ sex_female_pct: num 0
#>  $ race_ethnicity: chr "Korean (single-center Korean cohort)"
#>  $ disease_state : chr "Healthy adult male volunteers"
#>  $ dose_range    : chr "Single 50 mg oral dose of sumatriptan succinate (50 mg expressed as sumatriptan free base)"
#>  $ regions       : chr "Republic of Korea (Kyungpook National University Hospital Clinical Trial Center, Daegu)"
#>  $ notes         : chr "Retrospective re-analysis of the reference-formulation arm of a single-center, randomized, open-label, two-peri"| __truncated__

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

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 

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.

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.

# 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))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalka2', 'etalmtt', 'etalfr'

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.

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

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

nca_summary <- summary(nca_res)
knitr::kable(nca_summary,
             caption = "Simulated NCA parameters, 50 mg oral sumatriptan (n = 26 virtual subjects).")
Simulated NCA parameters, 50 mg oral sumatriptan (n = 26 virtual subjects).
Interval Start Interval End cohort N AUClast (h*ng/mL) Cmax (ng/mL) Tmax (h) Half-life (h) AUCinf,obs (h*ng/mL)
0 Inf 50 mg oral, single dose 26 NC 25.9 [24.2] 0.750 [0.500, 4.00] 2.44 [0.558] NC

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.
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))
#> Analytical AUCinf at CL/F = 418 L/h: 119.6 ng*h/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).