Skip to contents

Oxycodone + acetaminophen IR/ER combination tablets (Franke 2015)

Source: vignettes/articles/Franke_2015_oxycodone_acetaminophen.Rmd
Franke_2015_oxycodone_acetaminophen.Rmd

Model and source

Franke 2015 is a pooled post hoc population PK analysis of the two analytes delivered by one biphasic immediate-release / extended-release (IR/ER) oxycodone / acetaminophen 7.5/325 mg tablet. The authors fitted the two analytes separately and with different structural models, and they ran two separate analyses on two different cohorts:

  • a pooled fasted analysis of 151 completers from four phase 1 studies, and
  • a fed-vs-fasted food-effect analysis of the 31 completers of study 4.

That is four independent final models. Three of them are packaged here; the fourth is not, for the reason given under Errata and known deviations.

Packaged model Analyte Cohort Structure
Franke_2015_oxycodone oxycodone pooled fasted, n = 151 1-compartment, first-order absorption
Franke_2015_acetaminophen acetaminophen pooled fasted, n = 151 2-compartment, first-order absorption
Franke_2015_oxycodone_food oxycodone food effect, n = 31 1-compartment, first-order absorption
  • Article: https://doi.org/10.2147/DDDT.S79499
  • Citation: Franke RM, Morton T, Devarakonda K. Pooled post hoc analysis of population pharmacokinetics of oxycodone and acetaminophen following a single oral dose of biphasic immediate-release/extended-release oxycodone/acetaminophen tablets. Drug Des Devel Ther. 2015;9:4587-4597. doi:10.2147/DDDT.S79499 
mod_oc      <- readModelDb("Franke_2015_oxycodone")
mod_apap    <- readModelDb("Franke_2015_acetaminophen")
mod_oc_food <- readModelDb("Franke_2015_oxycodone_food")

Population

The pooled fasted analysis (Table 1) comprises the 151 of 251 enrolled participants who completed all treatment periods of four US phase 1 studies (protocols COV15000170, COV15000172, COV15000255, COV15000244). Participants were healthy adults aged 18-55 years with BMI 19 to <33 kg/m^2; study 3 enrolled healthy nondependent recreational users of prescription opioids. Mean age was 29.1 years (SD 8.6); 105 (69.5%) were men; 112 (74.2%) were White, 38 (25.2%) Black and 1 (0.6%) Asian. Median weight was 73.5 kg (range 50.0-120.9), median height 172.7 cm (range 153.6-194.0) and median BMI 25.1 kg/m^2 (range 18.6-32.9). Single oral doses were one, two or four intact tablets, i.e. 7.5/325, 15/650 or 30/1,300 mg of oxycodone/acetaminophen.

The food-effect analysis (Table 2) uses the 31 of 48 enrolled completers of study 4, a three-period six-sequence crossover in which each participant received two tablets (15/650 mg) fasted, after a low-fat meal (~25%-30% of kilocalories from fat, 800 +/- 80 kcal) and after a high-fat meal (~50% of kilocalories from fat, 1,000 +/- 100 kcal). Mean age was 30.8 years (SD 10.1); 21 (67.7%) were men; 23 (74.2%) White and 8 (25.8%) Black. Median weight was 72.3 kg (range 59.1-101.8) and median height 168.5 cm (range 155.1-185.5).

The same information is available programmatically from each model’s population metadata, e.g. readModelDb("Franke_2015_oxycodone")()$population.

Source trace

Every ini() entry carries an in-file comment naming its source location. They are collected here for review.

Franke 2015 prints its variability in a consistent and self-pinning style: each interindividual and residual entry is written as value (parenthetical) where the parenthetical is 100 * sqrt(value). For oxycodone, sqrt(0.0673) = 0.259, sqrt(0.0694) = 0.263 and sqrt(0.631) = 0.794 reproduce the printed 25.9, 26.3 and 79.4 exactly; for acetaminophen, sqrt(4.35) = 2.09 reproduces the printed 209.0. The first number of each pair is therefore a log-scale variance, and is used directly as the ~ value in ini().

Franke_2015_oxycodone (pooled fasted, n = 151)

Parameter Value Source location
lcl 92.4 L/h Table 3, final model, CL/F (L/hour); printed CL/F equation
lvc 772 L Table 3, final model, V/F (L); printed V/F equation
lka 1.15 1/h Table 3, final model, Ka (hour-1)
e_wt_cl 0.75 (fixed) printed CL/F equation, ((WT/73.45)**0.75)
WT exponent on V/F 1 (structural) printed V/F equation, bare (WT/73.45) ratio
e_race_black_cl 0.831 Table 3, CL/F ~ race; printed CL/F equation
e_race_black_vc 0.827 Table 3, V/F ~ race; printed V/F equation
race coding White 0 / Black 1 / Asian 2 n/a Figure 4 note
etalcl 0.0673 Table 3, interindividual CL/F = 0.0673 (25.9)
etalvc 0.0694 Table 3, interindividual V/F = 0.0694 (26.3)
etalka 0.631 Table 3, interindividual Ka = 0.631 (79.4)
addSd 2.63 ng/mL Table 3, intraindividual additive (SD) = 6.94 (2.63)
propSd 0.125 Table 3, intraindividual proportional (CV %) = 0.0156 (12.5)
d/dt(depot), d/dt(central) n/a Figure 1A, one-compartment schematic

Franke_2015_acetaminophen (pooled fasted, n = 151)

Parameter Value Source location
lcl 20.5 L/h Table 4, final model, CL/F (L/hour); printed CL/F equation
lvc 58.9 L Table 4, final model, V2/F (L); printed V2/F equation
lq 29.9 L/h Table 4, final model, Q (L/hour)
lvp 99.7 L Table 4, final model, V3/F (L)
lka 5.4 1/h (fixed) Table 4, Ka (hour-1) (fixed); Methods: “the Ka was fixed at 5.4 hour-1”
e_wt_cl 0.75 (fixed) printed CL/F equation, ((WT/73.45)**0.75)
WT exponent on V2/F 1 (structural) printed V2/F equation, bare (WT/73.45) ratio
etalcl 0.0554 printed CL/F equation exp(0.0554); Table 4 rounds to 0.06 (23.5)
etalvc 0.453 printed V2/F equation exp(0.453); Table 4 rounds to 0.45 (67.3)
etalq 4.35 Table 4, interindividual Q = 4.35 (209.0)
etalvp 0.75 Table 4, interindividual V3/F = 0.75 (86.3)
addSd 39.9 ng/mL Table 4, intraindividual additive (SD) = 1.59e+3 (39.9)
propSd 0.171 Table 4, intraindividual proportional (CV %) = 0.03 (17.1)
2-compartment ODE system n/a Figure 1B schematic; Methods names V1 = absorption, V2 = central, V3 = peripheral

Franke_2015_oxycodone_food (food effect, n = 31)

Parameter Value Source location
lcl 77.6 L/h Table 5, OC CL/F (L/hour), final model; printed CL/F equation
lvc 640 L Table 5, OC V/F (L), final model; printed V/F equation
lka 0.555 1/h Table 5, OC Ka (hour-1), final model; printed Ka equation
HT ratio on CL/F, reference 168.5 cm n/a printed CL/F equation, (HT/168.5)
WT ratio on V/F, reference 72.3 kg n/a printed V/F equation, (WT/72.3)
e_fed_lowfat_vc 0.608 printed V/F equation, (1 or 0.608 or 0.648)
e_fed_highfat_vc 0.648 printed V/F equation, (1 or 0.608 or 0.648)
e_fed_lowfat_ka 0.518 Table 5, OC Ka low fat; printed Ka equation
e_fed_highfat_ka 0.497 Table 5, OC Ka high fat; printed Ka equation
etalcl 0.0517 printed CL/F equation, exp(0.0517)
etalvc 0.0208 printed V/F equation, exp(0.0208)
etalka 0.104 printed Ka equation, exp(0.104)
addSd, propSd fixed(0) not reported; Table 5 footnote a states the structure only

Structural parameters reproduce the published tables

This first gate compares each packaged model’s parameters, evaluated at the paper’s own reference covariates, against the values transcribed by hand from Tables 3, 4 and 5. It is the transcription gate: a mistyped clearance, volume or rate constant fails here.

ref_grid <- function(mod, covs, amt) {
  ev <- rxode2::et(amt = amt, cmt = "depot") |>
    rxode2::et(c(0, 1), cmt = "central")
  d <- as.data.frame(ev)
  for (nm in names(covs)) d[[nm]] <- covs[[nm]]
  rxode2::rxSolve(mod, d, omega = NA) |> as.data.frame()
}

ref_oc <- ref_grid(mod_oc,
  list(WT = 73.45, RACE_BLACK = 0, RACE_ASIAN = 0), amt = 15)
#> ℹ parameter labels from comments will be replaced by 'label()'
ref_apap <- ref_grid(mod_apap, list(WT = 73.45), amt = 650)
#> ℹ parameter labels from comments will be replaced by 'label()'
ref_food <- ref_grid(mod_oc_food,
  list(HT = 168.5, WT = 72.3, FED_LOWFAT = 0, FED_HIGHFAT = 0), amt = 15)
#> ℹ parameter labels from comments will be replaced by 'label()'

structural <- tibble::tribble(
  ~Model,                        ~Parameter, ~Published, ~Model_value,
  "Franke_2015_oxycodone",       "CL/F (L/h)",   92.4,  ref_oc$cl[1],
  "Franke_2015_oxycodone",       "V/F (L)",     772,    ref_oc$vc[1],
  "Franke_2015_oxycodone",       "Ka (1/h)",      1.15, ref_oc$ka[1],
  "Franke_2015_acetaminophen",   "CL/F (L/h)",   20.5,  ref_apap$cl[1],
  "Franke_2015_acetaminophen",   "V2/F (L)",     58.9,  ref_apap$vc[1],
  "Franke_2015_acetaminophen",   "Q (L/h)",      29.9,  ref_apap$q[1],
  "Franke_2015_acetaminophen",   "V3/F (L)",     99.7,  ref_apap$vp[1],
  "Franke_2015_acetaminophen",   "Ka (1/h)",      5.4,  ref_apap$ka[1],
  "Franke_2015_oxycodone_food",  "CL/F (L/h)",   77.6,  ref_food$cl[1],
  "Franke_2015_oxycodone_food",  "V/F (L)",     640,    ref_food$vc[1],
  "Franke_2015_oxycodone_food",  "Ka (1/h)",      0.555, ref_food$ka[1]
) |>
  mutate(pct_diff = 100 * (Model_value - Published) / Published)

# Deterministic: these are the typical values with IIV suppressed, so they must
# match the printed table to within display rounding, not to a noise tolerance.
stopifnot(nrow(structural) == 11L)
stopifnot(max(abs(structural$pct_diff)) < 1e-8)

structural |>
  rename("Model file" = Model, "Published" = Published,
         "Model value" = Model_value, "% diff" = pct_diff) |>
  knitr::kable(digits = 4,
    caption = "Packaged typical values vs Franke 2015 Tables 3, 4 and 5.")
Packaged typical values vs Franke 2015 Tables 3, 4 and 5.
Model file Parameter Published Model value % diff
Franke_2015_oxycodone CL/F (L/h) 92.400 92.400 0
Franke_2015_oxycodone V/F (L) 772.000 772.000 0
Franke_2015_oxycodone Ka (1/h) 1.150 1.150 0
Franke_2015_acetaminophen CL/F (L/h) 20.500 20.500 0
Franke_2015_acetaminophen V2/F (L) 58.900 58.900 0
Franke_2015_acetaminophen Q (L/h) 29.900 29.900 0
Franke_2015_acetaminophen V3/F (L) 99.700 99.700 0
Franke_2015_acetaminophen Ka (1/h) 5.400 5.400 0
Franke_2015_oxycodone_food CL/F (L/h) 77.600 77.600 0
Franke_2015_oxycodone_food V/F (L) 640.000 640.000 0
Franke_2015_oxycodone_food Ka (1/h) 0.555 0.555 0

Covariate effects reproduce the published equations

The second gate re-implements each printed covariate equation independently, in plain R, and compares it against what the compiled model produces. Because the expected side is written from the paper rather than read back from the model, this gate can go red on a mis-transcribed exponent, reference value, or multiplier.

cov_check <- function(mod, covs, get, expected) {
  got <- ref_grid(mod, covs, amt = 15)[[get]][1]
  tibble(Model_value = got, Expected = expected,
         pct_diff = 100 * (got - expected) / expected)
}

cov_tab <- bind_rows(
  # --- pooled fasted oxycodone -------------------------------------------
  cov_check(mod_oc, list(WT = 73.45 * 1.1, RACE_BLACK = 0, RACE_ASIAN = 0),
            "cl", 92.4 * 1.1^0.75) |>
    mutate(Effect = "OC CL/F, WT +10% (power 0.75)"),
  cov_check(mod_oc, list(WT = 73.45 * 1.1, RACE_BLACK = 0, RACE_ASIAN = 0),
            "vc", 772 * 1.1) |>
    mutate(Effect = "OC V/F, WT +10% (power 1)"),
  cov_check(mod_oc, list(WT = 73.45, RACE_BLACK = 1, RACE_ASIAN = 0),
            "cl", 92.4 * 0.831) |>
    mutate(Effect = "OC CL/F, Black (RACE code 1)"),
  cov_check(mod_oc, list(WT = 73.45, RACE_BLACK = 1, RACE_ASIAN = 0),
            "vc", 772 * 0.827) |>
    mutate(Effect = "OC V/F, Black (RACE code 1)"),
  cov_check(mod_oc, list(WT = 73.45, RACE_BLACK = 0, RACE_ASIAN = 1),
            "cl", 92.4 * 0.831^2) |>
    mutate(Effect = "OC CL/F, Asian (RACE code 2, theta^RACE squares)"),
  cov_check(mod_oc, list(WT = 73.45, RACE_BLACK = 0, RACE_ASIAN = 1),
            "vc", 772 * 0.827^2) |>
    mutate(Effect = "OC V/F, Asian (RACE code 2, theta^RACE squares)"),
  # --- pooled fasted acetaminophen ---------------------------------------
  cov_check(mod_apap, list(WT = 73.45 * 1.1), "cl", 20.5 * 1.1^0.75) |>
    mutate(Effect = "APAP CL/F, WT +10% (power 0.75)"),
  cov_check(mod_apap, list(WT = 73.45 * 1.1), "vc", 58.9 * 1.1) |>
    mutate(Effect = "APAP V2/F, WT +10% (power 1)"),
  cov_check(mod_apap, list(WT = 73.45 * 1.1), "q", 29.9) |>
    mutate(Effect = "APAP Q, WT +10% (no weight term)"),
  cov_check(mod_apap, list(WT = 73.45 * 1.1), "vp", 99.7) |>
    mutate(Effect = "APAP V3/F, WT +10% (no weight term)"),
  # --- food-effect oxycodone ---------------------------------------------
  cov_check(mod_oc_food,
            list(HT = 168.5 * 1.1, WT = 72.3, FED_LOWFAT = 0, FED_HIGHFAT = 0),
            "cl", 77.6 * 1.1) |>
    mutate(Effect = "OC-food CL/F, HT +10% (linear ratio)"),
  cov_check(mod_oc_food,
            list(HT = 168.5, WT = 72.3, FED_LOWFAT = 1, FED_HIGHFAT = 0),
            "ka", 0.555 * 0.518) |>
    mutate(Effect = "OC-food Ka, low-fat meal"),
  cov_check(mod_oc_food,
            list(HT = 168.5, WT = 72.3, FED_LOWFAT = 0, FED_HIGHFAT = 1),
            "ka", 0.555 * 0.497) |>
    mutate(Effect = "OC-food Ka, high-fat meal"),
  cov_check(mod_oc_food,
            list(HT = 168.5, WT = 72.3, FED_LOWFAT = 1, FED_HIGHFAT = 0),
            "vc", 640 * 0.608) |>
    mutate(Effect = "OC-food V/F, low-fat meal"),
  cov_check(mod_oc_food,
            list(HT = 168.5, WT = 72.3, FED_LOWFAT = 0, FED_HIGHFAT = 1),
            "vc", 640 * 0.648) |>
    mutate(Effect = "OC-food V/F, high-fat meal")
) |>
  select(Effect, Expected, Model_value, pct_diff)

stopifnot(nrow(cov_tab) == 15L)
stopifnot(max(abs(cov_tab$pct_diff)) < 1e-8)

cov_tab |>
  rename("Covariate effect" = Effect, "Expected (paper equation)" = Expected,
         "Model" = Model_value, "% diff" = pct_diff) |>
  knitr::kable(digits = 4,
    caption = "Model covariate effects vs the printed Franke 2015 equations, recomputed independently.")
Model covariate effects vs the printed Franke 2015 equations, recomputed independently.
Covariate effect Expected (paper equation) Model % diff
OC CL/F, WT +10% (power 0.75) 99.2468 99.2468 0
OC V/F, WT +10% (power 1) 849.2000 849.2000 0
OC CL/F, Black (RACE code 1) 76.7844 76.7844 0
OC V/F, Black (RACE code 1) 638.4440 638.4440 0
OC CL/F, Asian (RACE code 2, theta^RACE squares) 63.8078 63.8078 0
OC V/F, Asian (RACE code 2, theta^RACE squares) 527.9932 527.9932 0
APAP CL/F, WT +10% (power 0.75) 22.0190 22.0190 0
APAP V2/F, WT +10% (power 1) 64.7900 64.7900 0
APAP Q, WT +10% (no weight term) 29.9000 29.9000 0
APAP V3/F, WT +10% (no weight term) 99.7000 99.7000 0
OC-food CL/F, HT +10% (linear ratio) 85.3600 85.3600 0
OC-food Ka, low-fat meal 0.2875 0.2875 0
OC-food Ka, high-fat meal 0.2758 0.2758 0
OC-food V/F, low-fat meal 389.1200 389.1200 0
OC-food V/F, high-fat meal 414.7200 414.7200 0

Virtual cohort

Original observed data are not publicly available. The cohorts below sample covariates from the published baseline distributions (Table 1 for the pooled fasted analysis, Table 2 for the food-effect analysis), truncated to the reported ranges.

# set.seed() seeds R's RNG. It does NOT seed rxode2's simulation RNG, whose
# streams are partitioned PER SOLVER THREAD -- so this cohort is reproducible
# on this machine and different on one with a different thread count. Every
# assertion below is written to hold for any cohort the model can produce.
set.seed(20150801)

n_arm <- 100L

rtrunc_norm <- function(n, mean, sd, lo, hi) {
  x <- rnorm(n, mean, sd)
  pmin(pmax(x, lo), hi)
}

# Sampling grid: dense through absorption (acetaminophen Ka is 5.4 1/h, so the
# peak is early and a coarse grid understates AUC), then out to 48 h, matching
# the studies' 48 h sampling window.
tgrid <- unique(c(seq(0, 2, by = 0.1), seq(2.25, 12, by = 0.25),
                  seq(13, 48, by = 1)))

# --- pooled fasted cohort (Table 1) --------------------------------------
# Oxycodone and acetaminophen are delivered by the SAME tablet to the SAME
# participants, so both analytes reuse one set of covariate draws.
make_fasted <- function(n, oc_mg, apap_mg, id_offset = 0L) {
  race <- sample(c("White", "Black", "Asian"), n, replace = TRUE,
                 prob = c(0.742, 0.252, 0.006))
  tibble(
    id         = id_offset + seq_len(n),
    WT         = rtrunc_norm(n, 75.6, 13.6, 50.0, 120.9),
    RACE_BLACK = as.integer(race == "Black"),
    RACE_ASIAN = as.integer(race == "Asian"),
    race       = race,
    treatment  = sprintf("%g/%g mg", oc_mg, apap_mg),
    oc_mg      = oc_mg,
    apap_mg    = apap_mg
  )
}

fasted_subj <- bind_rows(
  make_fasted(n_arm,  7.5,  325, id_offset =        0L),
  make_fasted(n_arm, 15.0,  650, id_offset =   n_arm   ),
  make_fasted(n_arm, 30.0, 1300, id_offset = 2L * n_arm)
)

expand_events <- function(subj, amt_col) {
  doses <- subj |>
    mutate(time = 0, amt = .data[[amt_col]], evid = 1L, cmt = "depot")
  obs <- subj |>
    tidyr::crossing(time = tgrid) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "central")
  bind_rows(doses, obs) |> arrange(id, time, desc(evid))
}

ev_oc   <- expand_events(fasted_subj, "oc_mg")
ev_apap <- expand_events(fasted_subj, "apap_mg")
stopifnot(!anyDuplicated(unique(ev_oc[, c("id", "time", "evid")])))

# --- food-effect cohort (Table 2) ----------------------------------------
# Study 4 was a three-period crossover: the SAME participant was dosed under
# all three prandial states. Each arm therefore reuses one set of covariate
# draws, with disjoint IDs so rxSolve keeps the arms separate.
food_base <- tibble(
  subject = seq_len(n_arm),
  WT      = rtrunc_norm(n_arm, 74.9, 12.0, 59.1, 101.8),
  HT      = rtrunc_norm(n_arm, 169.8, 8.5, 155.1, 185.5)
)

food_subj <- bind_rows(
  food_base |> mutate(treatment = "Fasted",   FED_LOWFAT = 0L, FED_HIGHFAT = 0L,
                      id = subject),
  food_base |> mutate(treatment = "Low fat",  FED_LOWFAT = 1L, FED_HIGHFAT = 0L,
                      id = subject + n_arm),
  food_base |> mutate(treatment = "High fat", FED_LOWFAT = 0L, FED_HIGHFAT = 1L,
                      id = subject + 2L * n_arm)
) |>
  mutate(oc_mg = 15, treatment = factor(treatment,
                                        levels = c("Fasted", "Low fat", "High fat")))

ev_food <- expand_events(food_subj, "oc_mg")
stopifnot(!anyDuplicated(unique(ev_food[, c("id", "time", "evid")])))

Simulation 

rxode2::rxSetSeed(20150801)
sim_oc <- rxode2::rxSolve(
  mod_oc, events = ev_oc,
  keep = c("treatment", "WT", "RACE_BLACK", "RACE_ASIAN", "race", "oc_mg")
) |> as.data.frame()

rxode2::rxSetSeed(20150802)
sim_apap <- rxode2::rxSolve(
  mod_apap, events = ev_apap,
  keep = c("treatment", "WT", "apap_mg")
) |> as.data.frame()

rxode2::rxSetSeed(20150803)
sim_food <- rxode2::rxSolve(
  mod_oc_food, events = ev_food,
  keep = c("treatment", "WT", "HT", "FED_LOWFAT", "FED_HIGHFAT", "oc_mg")
) |> as.data.frame()

# Concentrations must stay non-negative for PKNCA's terminal-slope log() fit.
stopifnot(all(sim_oc$Cc   >= 0, na.rm = TRUE))
stopifnot(all(sim_apap$Cc >= 0, na.rm = TRUE))
stopifnot(all(sim_food$Cc >= 0, na.rm = TRUE))

Replicate published figures

Figure 1 of Franke 2015 is a pair of structural schematics, and Figures 2, 3, 5 and 6 are goodness-of-fit plots against observed concentrations that are not publicly available. What can be replicated from the packaged models is the concentration-time behaviour implied by each final model, and the interindividual Ka variability by prandial state shown in Figure 7.

prof <- bind_rows(
  sim_oc   |> mutate(analyte = "Oxycodone (pooled fasted)"),
  sim_apap |> mutate(analyte = "Acetaminophen (pooled fasted)")
) |>
  filter(!is.na(Cc)) |>
  group_by(analyte, treatment, time) |>
  summarise(Q05 = quantile(Cc, 0.05), Q50 = quantile(Cc, 0.50),
            Q95 = quantile(Cc, 0.95), .groups = "drop")

ggplot(prof, aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.18, colour = NA) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~analyte, scales = "free_y") +
  labs(x = "Time (h)", y = "Concentration (ng/mL)",
       colour = "Dose (OC/APAP)", fill = "Dose (OC/APAP)",
       title = "Simulated single-dose profiles, pooled fasted models",
       caption = "Median with 5th-95th percentile band, 100 subjects per dose arm.")

# Replicates Figure 7 of Franke 2015: interindividual Ka variability for
# oxycodone under fasted vs low-fat vs high-fat conditions, final model.
sim_food |>
  distinct(id, treatment, ka) |>
  ggplot(aes(treatment, ka, fill = treatment)) +
  geom_boxplot(alpha = 0.6, outlier.alpha = 0.4) +
  labs(x = NULL, y = "Individual Ka (1/h)",
       title = "Figure 7 - interindividual Ka variability by prandial state",
       caption = "Replicates the final-model panel of Figure 7 of Franke 2015.") +
  theme(legend.position = "none")

The food contrast is shown at typical values, not as a median across the simulated arms. rxSolve draws random effects per id, and the three prandial arms necessarily carry disjoint IDs, so their etas are independent rather than paired as they were in study 4’s within-subject crossover. Clearance in this model does not depend on diet at all, yet independent eta draws still move the median AUC between arms by several percent – noise that would be read as a food effect. Zeroing the random effects removes it, leaving exactly the published covariate relationships.

ev_food_typ <- bind_rows(
  tibble(id = 1L, treatment = "Fasted",   FED_LOWFAT = 0L, FED_HIGHFAT = 0L),
  tibble(id = 2L, treatment = "Low fat",  FED_LOWFAT = 1L, FED_HIGHFAT = 0L),
  tibble(id = 3L, treatment = "High fat", FED_LOWFAT = 0L, FED_HIGHFAT = 1L)
) |>
  mutate(HT = 168.5, WT = 72.3, oc_mg = 15,
         treatment = factor(treatment, levels = c("Fasted", "Low fat", "High fat"))) |>
  expand_events("oc_mg")

sim_food_typ <- rxode2::rxSolve(
  rxode2::zeroRe(mod_oc_food), events = ev_food_typ,
  keep = c("treatment")
) |> as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalka'
#> Warning: multi-subject simulation without without 'omega'

# Clearance carries no diet term, so AUC0-inf = Dose/CL must be identical in
# all three arms. This is the paired contrast the arm medians cannot give.
auc_typ <- sim_food_typ |> filter(!is.na(Cc)) |> group_by(treatment) |>
  summarise(cl = first(cl), vc = first(vc),
            auc48 = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
            clast = tail(Cc, 1), .groups = "drop") |>
  mutate(aucinf = auc48 + clast / (cl / vc),
         auc_theory = 15 / cl * 1e3,
         pct_diff = 100 * (aucinf - auc_theory) / auc_theory)

stopifnot(nrow(auc_typ) == 3L)
stopifnot(length(unique(round(auc_typ$cl, 9))) == 1L)
# Deterministic (zeroRe), so a tight bound is correct: realised max |diff|
# 0.033% and an across-arm AUC0-inf spread of 0.012%. Note that AUC over the
# finite 0-48 h window is NOT equal across arms -- it spreads 0.37%, because
# slower absorption under a meal leaves more drug unabsorbed at 48 h. The
# published "total exposure is unchanged" claim is about AUC0-inf.
stopifnot(max(abs(auc_typ$pct_diff)) < 0.5)
stopifnot(max(auc_typ$aucinf) / min(auc_typ$aucinf) < 1.005)

sim_food_typ |>
  filter(!is.na(Cc)) |>
  ggplot(aes(time, Cc, colour = treatment)) +
  geom_line(linewidth = 0.7) +
  coord_cartesian(xlim = c(0, 24)) +
  labs(x = "Time (h)", y = "Oxycodone (ng/mL)", colour = NULL,
       title = "Food effect on the oxycodone profile, typical values",
       caption = paste("15 mg oxycodone at HT 168.5 cm and WT 72.3 kg. A meal slows",
                       "absorption and lowers the apparent volume, so the peak is",
                       "later; total exposure is unchanged."))

PKNCA validation 

run_nca <- function(sim, ev, amt_col) {
  sim_nca <- sim |>
    filter(!is.na(Cc)) |>
    select(id, time, Cc, treatment)

  # Guarantee a time-zero record per (id, treatment); pre-dose Cc = 0 is
  # correct for an extravascular single dose.
  sim_nca <- bind_rows(
    sim_nca,
    sim_nca |> distinct(id, treatment) |> mutate(time = 0, Cc = 0)
  ) |>
    distinct(id, treatment, time, .keep_all = TRUE) |>
    arrange(id, treatment, time)

  dose_df <- ev |>
    filter(evid == 1) |>
    mutate(amt = .data[[amt_col]]) |>
    select(id, time, amt, treatment)

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

  PKNCA::pk.nca(PKNCA::PKNCAdata(
    PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id),
    PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id),
    intervals = intervals
  ))
}

nca_oc   <- run_nca(sim_oc,   ev_oc,   "oc_mg")
nca_apap <- run_nca(sim_apap, ev_apap, "apap_mg")
nca_food <- run_nca(sim_food, ev_food, "oc_mg")

nca_wide <- function(res) {
  as.data.frame(res) |>
    filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "aucinf.obs", "half.life")) |>
    select(treatment, id, PPTESTCD, PPORRES) |>
    tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)
}

nca_oc_w   <- nca_wide(nca_oc)
nca_apap_w <- nca_wide(nca_apap)
nca_food_w <- nca_wide(nca_food)
stopifnot(nrow(nca_oc_w) == 3L * n_arm, nrow(nca_apap_w) == 3L * n_arm,
          nrow(nca_food_w) == 3L * n_arm)
nca_summary <- function(w, label) {
  w |>
    group_by(treatment) |>
    summarise(n = n(),
              cmax = median(cmax), tmax = median(tmax),
              aucinf.obs = median(aucinf.obs), half.life = median(half.life),
              .groups = "drop") |>
    mutate(Analysis = label, .before = 1)
}

bind_rows(
  nca_summary(nca_oc_w,   "Oxycodone, pooled fasted"),
  nca_summary(nca_apap_w, "Acetaminophen, pooled fasted"),
  nca_summary(nca_food_w, "Oxycodone, food effect (15 mg)")
) |>
  rename("Analysis" = Analysis, "Arm" = treatment, "N" = n,
         "Cmax (ng/mL)" = cmax, "Tmax (h)" = tmax,
         "AUC0-inf (ng*h/mL)" = aucinf.obs, "t1/2 (h)" = half.life) |>
  knitr::kable(digits = 2,
    caption = "Median NCA parameters by arm (PKNCA, 100 subjects per arm).")
Median NCA parameters by arm (PKNCA, 100 subjects per arm).
Analysis Arm N Cmax (ng/mL) Tmax (h) AUC0-inf (ng*h/mL) t1/2 (h)
Oxycodone, pooled fasted 15/650 mg 100 14.51 2.00 161.58 6.03
Oxycodone, pooled fasted 30/1300 mg 100 29.70 2.38 340.17 5.99
Oxycodone, pooled fasted 7.5/325 mg 100 7.93 2.12 91.01 6.12
Acetaminophen, pooled fasted 15/650 mg 100 6800.45 0.50 30497.86 9.21
Acetaminophen, pooled fasted 30/1300 mg 100 15109.37 0.40 60318.19 10.19
Acetaminophen, pooled fasted 7.5/325 mg 100 3061.06 0.45 15352.82 9.19
Oxycodone, food effect (15 mg) Fasted 100 14.91 3.50 180.40 5.41
Oxycodone, food effect (15 mg) Low fat 100 15.89 4.25 188.10 3.81
Oxycodone, food effect (15 mg) High fat 100 15.94 4.25 187.22 3.95

Franke 2015 reports no NCA table of its own – it is a population-PK analysis whose companion NCA appears in the separate primary PK publication (Devarakonda 2014, Drug Des Devel Ther 8:1125-1134), which is not on disk here. In place of a published-NCA comparison the sections below use closed-form identities that the model must satisfy exactly.

Mass balance: AUC0-inf x CL/F = dose

For a linear model with first-order elimination, AUC(0, Inf) * CL/F equals the administered dose exactly, per subject and independently of the absorption model. This validates the ODE system and the ng/mL unit scaling together; the residual gap is trapezoidal-integration and extrapolation error only.

mass_balance <- function(w, sim, dose_col) {
  per_id <- sim |> distinct(id, cl, .data[[dose_col]]) |>
    rename(dose_mg = all_of(dose_col))
  w |>
    left_join(per_id, by = "id") |>
    # aucinf.obs is ng*h/mL; cl is L/h; dose is mg. 1 ng/mL = 1e-3 mg/L, so
    # AUC(ng*h/mL) * CL(L/h) / 1e3 recovers the dose in mg.
    mutate(recovered_mg = aucinf.obs * cl / 1e3,
           pct_diff = 100 * (recovered_mg - dose_mg) / dose_mg)
}

mb <- bind_rows(
  mass_balance(nca_oc_w,   sim_oc,   "oc_mg")   |> mutate(Analysis = "Oxycodone, pooled fasted"),
  mass_balance(nca_apap_w, sim_apap, "apap_mg") |> mutate(Analysis = "Acetaminophen, pooled fasted"),
  mass_balance(nca_food_w, sim_food, "oc_mg")   |> mutate(Analysis = "Oxycodone, food effect")
)

stopifnot(nrow(mb) == 9L * n_arm, !anyNA(mb$pct_diff))

# The two one-compartment oxycodone models integrate essentially exactly: both
# sides use the same drawn parameters, so the gap is pure trapezoidal and
# extrapolation error and a tight bound is correct. Realised max |diff| 0.084%
# (pooled fasted) and 0.172% (food effect) at 16 threads. A 1% bound keeps
# headroom for a different cohort draw and still goes red on a unit-factor
# error, which is off by 1000x.
oc_rows <- mb$Analysis != "Acetaminophen, pooled fasted"
stopifnot(sum(oc_rows) == 6L * n_arm)
stopifnot(max(abs(mb$pct_diff[oc_rows])) < 1)

# Acetaminophen is different in kind. Its IIV on Q is CV 209% and on V3/F is
# CV 86% (Table 4), so a minority of subjects draw a deep peripheral
# compartment whose terminal phase the 48 h sampling window cannot resolve --
# realised terminal half-lives ran from 5 h to over 1,700 h -- and aucinf.obs
# then under-recovers the dose. That is a per-subject physical mechanism, not
# numerical noise, so the centre and a robust quantile are asserted rather than
# the extreme. Realised median -0.20%, 90th pct 0.52%, 99th pct 1.60%,
# max 2.99%.
ap_rows <- !oc_rows
stopifnot(sum(ap_rows) == 3L * n_arm)
stopifnot(abs(median(mb$pct_diff[ap_rows])) < 1.5)
stopifnot(quantile(abs(mb$pct_diff[ap_rows]), 0.9) < 3)
stopifnot(max(abs(mb$pct_diff[ap_rows])) < 25)

mb |>
  group_by(Analysis, treatment) |>
  summarise(`Median % diff` = median(pct_diff),
            `Max |% diff|` = max(abs(pct_diff)), .groups = "drop") |>
  rename("Arm" = treatment) |>
  knitr::kable(digits = 3,
    caption = "Dose recovered as AUC0-inf x CL/F, versus the administered dose.")
Dose recovered as AUC0-inf x CL/F, versus the administered dose.
Analysis Arm Median % diff Max |% diff|
Acetaminophen, pooled fasted 15/650 mg -0.204 16.101
Acetaminophen, pooled fasted 30/1300 mg -0.213 9.773
Acetaminophen, pooled fasted 7.5/325 mg -0.187 8.309
Oxycodone, food effect Fasted -0.016 0.067
Oxycodone, food effect High fat -0.022 0.059
Oxycodone, food effect Low fat -0.023 0.048
Oxycodone, pooled fasted 15/650 mg -0.019 0.185
Oxycodone, pooled fasted 30/1300 mg -0.017 0.091
Oxycodone, pooled fasted 7.5/325 mg -0.017 0.159

The acetaminophen arms recover the dose less tightly than the oxycodone arms, and that is a property of the published model rather than of the implementation. Table 4 gives the interindividual variability on Q as CV 209% and on V3/F as CV 86%, which are very large; a minority of simulated subjects therefore draw a deep, slowly-equilibrating peripheral compartment whose terminal phase a 48 h sampling window cannot characterise, and aucinf.obs under-recovers their dose by a few percent. Users simulating this model for exposure metrics that depend on the terminal phase should extend the observation window well past 48 h, or suppress the Q and V3/F random effects.

Closed-form Cmax and Tmax for the one-compartment models

Both oxycodone models are one-compartment with first-order absorption, so Tmax = ln(ka/kel)/(ka - kel) and Cmax = (Dose/V) * exp(-kel * Tmax) have exact analytical forms. Comparing the NCA of the solved ODE against those forms checks the ODE structure independently of the parameter transcription already gated above.

closed_form <- function(w, sim, dose_col, label) {
  per_id <- sim |> distinct(id, cl, vc, ka, .data[[dose_col]]) |>
    rename(dose_mg = all_of(dose_col))
  w |>
    left_join(per_id, by = "id") |>
    mutate(kel = cl / vc,
           tmax_cf = log(ka / kel) / (ka - kel),
           # dose_mg / vc is mg/L = ug/mL; x1000 gives ng/mL.
           cmax_cf = (dose_mg / vc) * 1000 *
             (exp(-kel * tmax_cf) - exp(-ka * tmax_cf)) * ka / (ka - kel),
           cmax_pct = 100 * (cmax - cmax_cf) / cmax_cf,
           Analysis = label)
}

cf <- bind_rows(
  closed_form(nca_oc_w,   sim_oc,   "oc_mg", "Oxycodone, pooled fasted"),
  closed_form(nca_food_w, sim_food, "oc_mg", "Oxycodone, food effect")
)

stopifnot(nrow(cf) == 6L * n_arm, !anyNA(cf$cmax_pct))

# NCA Cmax is the maximum over a DISCRETE grid, so it sits at or below the
# analytical peak by an amount set by the local grid spacing -- a per-subject
# physical mechanism, so this is asserted on the centre and a robust quantile
# rather than on the extreme. Realised median -0.06% / 90th pct 0.9% at 16
# threads. A mis-specified absorption arm moves Cmax by tens of percent.
stopifnot(abs(median(cf$cmax_pct)) < 2)
stopifnot(quantile(abs(cf$cmax_pct), 0.9) < 6)

cf |>
  group_by(Analysis, treatment) |>
  summarise(`Median % diff in Cmax` = median(cmax_pct),
            `90th pct |% diff|` = quantile(abs(cmax_pct), 0.9),
            `Median Tmax, NCA (h)` = median(tmax),
            `Median Tmax, closed form (h)` = median(tmax_cf),
            .groups = "drop") |>
  rename("Arm" = treatment) |>
  knitr::kable(digits = 3,
    caption = "Simulated Cmax and Tmax versus the one-compartment closed form.")
Simulated Cmax and Tmax versus the one-compartment closed form.
Analysis Arm Median % diff in Cmax 90th pct |% diff| Median Tmax, NCA (h) Median Tmax, closed form (h)
Oxycodone, food effect Fasted -0.009 0.044 3.500 3.497
Oxycodone, food effect High fat -0.010 0.034 4.250 4.262
Oxycodone, food effect Low fat -0.007 0.035 4.250 4.225
Oxycodone, pooled fasted 15/650 mg -0.011 0.050 2.000 1.977
Oxycodone, pooled fasted 30/1300 mg -0.013 0.051 2.375 2.348
Oxycodone, pooled fasted 7.5/325 mg -0.010 0.044 2.125 2.104

Dose proportionality

All three models are linear, so exposure must scale exactly with dose.

dp <- nca_oc_w |>
  left_join(sim_oc |> distinct(id, oc_mg), by = "id") |>
  mutate(auc_per_mg = aucinf.obs / oc_mg)

dp_med <- dp |> group_by(treatment) |> summarise(med = median(auc_per_mg),
                                                 .groups = "drop")
stopifnot(nrow(dp_med) == 3L)
# Dose-normalised AUC is identical per subject across arms by linearity; the
# spread across arms here is only the between-cohort covariate draw.
stopifnot(max(dp_med$med) / min(dp_med$med) < 1.25)

dp_med |>
  rename("Arm" = treatment, "Median AUC0-inf / mg (ng*h/mL/mg)" = med) |>
  knitr::kable(digits = 1,
    caption = "Dose-normalised oxycodone exposure across the three dose arms.")
Dose-normalised oxycodone exposure across the three dose arms.
Arm Median AUC0-inf / mg (ng*h/mL/mg)
15/650 mg 10.8
30/1300 mg 11.3
7.5/325 mg 12.1

Published claims

Franke 2015 states several quantitative covariate effects in prose. The table below evaluates each against the packaged models. Rows marked as deviations are places where the paper’s prose disagrees with its own tables and equations; in every such case the model follows the table and the equation, which agree with one another, as documented under Errata below.

pct_change <- function(mod, base, alt, what) {
  100 * (ref_grid(mod, alt, 15)[[what]][1] / ref_grid(mod, base, 15)[[what]][1] - 1)
}

oc_base   <- list(WT = 73.45, RACE_BLACK = 0, RACE_ASIAN = 0)
food_base_c <- list(HT = 168.5, WT = 72.3, FED_LOWFAT = 0, FED_HIGHFAT = 0)

claims <- tibble::tribble(
  ~Claim, ~Source, ~Stated, ~Achieved, ~Deviation,
  "OC: +10% WT changes CL/F by ~7.5%", "Abstract / Results", 7.5,
    pct_change(mod_oc, oc_base, modifyList(oc_base, list(WT = 73.45 * 1.1)), "cl"), FALSE,
  "OC: +10% WT changes V/F by ~7.5%", "Abstract / Results", 7.5,
    pct_change(mod_oc, oc_base, modifyList(oc_base, list(WT = 73.45 * 1.1)), "vc"), TRUE,
  "OC: Black CL/F is 17.3% lower", "Results", -17.3,
    pct_change(mod_oc, oc_base, modifyList(oc_base, list(RACE_BLACK = 1)), "cl"), TRUE,
  "OC: Black V/F is 16.9% lower", "Results", -16.9,
    pct_change(mod_oc, oc_base, modifyList(oc_base, list(RACE_BLACK = 1)), "vc"), TRUE,
  "APAP: +10% WT changes CL/F by 7.5%", "Results", 7.5,
    pct_change(mod_apap, list(WT = 73.45), list(WT = 73.45 * 1.1), "cl"), FALSE,
  "APAP: +10% WT changes V2/F by 7.5%", "Results", 7.5,
    pct_change(mod_apap, list(WT = 73.45), list(WT = 73.45 * 1.1), "vc"), TRUE,
  "OC: low-fat meal decreases Ka by 39%", "Discussion", -39,
    pct_change(mod_oc_food, food_base_c,
               modifyList(food_base_c, list(FED_LOWFAT = 1)), "ka"), TRUE,
  "OC: high-fat meal decreases Ka by 48%", "Discussion", -48,
    pct_change(mod_oc_food, food_base_c,
               modifyList(food_base_c, list(FED_HIGHFAT = 1)), "ka"), TRUE
) |>
  mutate(Gap = Achieved - Stated,
         Pass = Deviation | abs(Gap) < 0.5)

# A gate that cannot go red is worse than none: confirm there are non-deviation
# rows to test before asserting on them.
stopifnot(sum(!claims$Deviation) >= 2L)
stopifnot(all(claims$Pass[!claims$Deviation]))

claims |>
  rename("Stated (%)" = Stated, "Model (%)" = Achieved,
         "Gap (pp)" = Gap, "Known deviation" = Deviation) |>
  knitr::kable(digits = 2,
    caption = "Franke 2015 prose claims evaluated against the packaged models.")
Franke 2015 prose claims evaluated against the packaged models.
Claim Source Stated (%) Model (%) Known deviation Gap (pp) Pass
OC: +10% WT changes CL/F by ~7.5% Abstract / Results 7.5 7.41 FALSE -0.09 TRUE
OC: +10% WT changes V/F by ~7.5% Abstract / Results 7.5 10.00 TRUE 2.50 TRUE
OC: Black CL/F is 17.3% lower Results -17.3 -16.90 TRUE 0.40 TRUE
OC: Black V/F is 16.9% lower Results -16.9 -17.30 TRUE -0.40 TRUE
APAP: +10% WT changes CL/F by 7.5% Results 7.5 7.41 FALSE -0.09 TRUE
APAP: +10% WT changes V2/F by 7.5% Results 7.5 10.00 TRUE 2.50 TRUE
OC: low-fat meal decreases Ka by 39% Discussion -39.0 -48.20 TRUE -9.20 TRUE
OC: high-fat meal decreases Ka by 48% Discussion -48.0 -50.30 TRUE -2.30 TRUE

Errata and known deviations

The acetaminophen food-effect model is not packaged. Franke 2015 states (Table 5 footnote a) that a two-compartment additive-and-proportional error model was used for acetaminophen in the fed-vs-fasted analysis, and Table 5 reports its CL/F (22.9 L/h), V2/F (140 L), Ka (3.17 1/h), sex factors on CL/F (0.278) and V2/F (0.295), and meal factors on Ka (0.802 low fat, 0.825 high fat). It does not report Q or V3/F – neither in the base model column nor in the final model column, and neither appears in the printed equations. A two-compartment model cannot be simulated without them, and the values cannot be borrowed from the pooled fasted acetaminophen model because that fit’s central volume (58.9 L) differs from this one’s (140 L) by a factor of 2.4, so its disposition parameters are not transferable. The model is therefore left out rather than completed with invented or transplanted numbers.

The paper’s prose percentages disagree with its own tables and equations in four places. In each case the table and the printed equation agree with each other, and the packaged models follow them.

  • Race effects are transposed. The Results text says Black participants had “a 17.3% lower CL/F and a 16.9% lower V/F”. Table 3 gives CL/F ~ race = 0.831 (16.9% lower) and V/F ~ race = 0.827 (17.3% lower), and the printed CL/F and V/F equations carry the same two multipliers in the same assignment. Two printed sources agree against one prose sentence.
  • The weight effect on volume is 10%, not 7.5%. The Abstract and Results say a 10% change in weight gives “~7.5% change in total body clearance (CL/F) and volume of distribution (V/F)”. The printed equations give CL/F a power of 0.75 (so 1.1^0.75 = +7.4%) but give V/F a bare (WT/73.45) ratio, i.e. a power of 1 (so +10.0%). The 7.5% figure is the clearance result quoted for both parameters. This applies to the acetaminophen model equally.
  • The oxycodone food effect on Ka does not match the stated percentages. The Discussion says the oxycodone Ka “decreased 39%” on a low-fat meal and “48%” on a high-fat meal. Table 5 and the printed Ka equation both give multipliers of 0.518 (low fat, i.e. 48.2% lower) and 0.497 (high fat, i.e. 50.3% lower). The corresponding acetaminophen multipliers, 0.802 and 0.825, reproduce that model’s stated 20% and 18% reductions to within 0.2 and 0.5 percentage points – which is what establishes that these numbers are dimensionless multipliers rather than the absolute rate constants that Table 5’s (hour-1) column label would suggest. Read as absolute rates they would imply oxycodone food effects of only 6.7% and 10.5%, contradicting the Discussion far more severely. The models use the multiplier reading.

Residual error is not reported for the food-effect analysis. Table 5 names the error structure but gives no magnitude, in either the base or the final model column. Franke_2015_oxycodone_food therefore declares addSd <- fixed(0) and propSd <- fixed(0) rather than carrying over the unrelated 151-subject pooled-fasted estimates. Simulations from that model are IPRED-only; the interindividual variability, which is reported, is intact.

The acetaminophen absorption rate constant is fixed, not estimated. Methods state it “was fixed at 5.4 hour-1 for the single-dose population PK of APAP” to stabilise the model, so lka is wrapped in fixed(). It was not fixed in the food-effect analysis.

Race is coded 0/1/2 and enters as a power. The oxycodone fasted model’s race term is written by the authors as theta^RACE, and the Figure 4 note gives the dataset coding as White = 0, Black = 1, Asian = 2. An Asian participant therefore receives the square of the Black multiplier. The packaged model reproduces this exactly, rebuilding the integer code from the canonical binary indicators as RACE_BLACK + 2 * RACE_ASIAN. This is an artifact of applying a power form to an ordinal code, not a fitted Asian effect: only one of 151 participants was Asian, and the paper’s Discussion states that PK variability in Asian participants cannot be ruled out from these data.

Reference weight is quoted two ways. The printed equations use 73.45 kg and Table 1 gives the cohort median as 73.5 kg; the acetaminophen Results text quotes 73.5 while its equation uses 73.45. The models use 73.45, the value in the equations.

Table 1’s weight range contradicts the stated inclusion criterion. The Participants section requires a body weight of at least 59 kg, but Table 1 reports a minimum of 50.0 kg among completers. The virtual cohort follows Table 1.

Assumptions and deviations

  • Observed concentrations are not publicly available, so all figures use virtual cohorts of 100 subjects per arm drawn from the Table 1 and Table 2 baseline distributions, truncated to the reported ranges. Weight and height are sampled independently from truncated normal distributions matching the published means and SDs; the paper reports no correlation structure between them.
  • Race is sampled at the Table 1 proportions (74.2% White, 25.2% Black, 0.6% Asian). With 100 subjects per arm, a given arm may contain no Asian participant at all – as in the source cohort, which contained one.
  • The food-effect cohort reuses one set of covariate draws across the three prandial arms, echoing study 4’s three-period crossover design in which the same participant received all three treatments. Arm IDs are offset so rxSolve keeps them separate – which means the random effects are drawn independently per arm and the arms are not a paired crossover. Any between-arm comparison of arm medians therefore carries eta noise: the model has no diet term on clearance, yet the arm-median AUC still moves by several percent. The food contrast is consequently shown at typical values (rxode2::zeroRe()), and the diet effects are gated deterministically in the covariate-effects section rather than read off the simulated arms.
  • Sex, age and BMI were screened by the authors but retained in none of the three packaged models; they are recorded in each model’s covariatesDataExcluded metadata rather than in covariateData, and are not sampled in the virtual cohorts.
  • All parameter values come from the paper’s own tables and printed equations. No value was digitised from a figure, obtained by correspondence, or carried from another publication.