---
title: "Febuxostat (Chen 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Febuxostat (Chen 2024)}
  %\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 model}
mod <- readModelDb("Chen_2024_febuxostat")
ui <- rxode2::rxode2(mod)
```

* Citation: `r ui$reference`
* Description: `r ui$description`
* Article: <https://doi.org/10.1186/s40360-024-00783-1>
* Supplementary Material 1 (study design, model-building table, Caucasian
  simulation table): <https://doi.org/10.1186/s40360-024-00783-1>

Febuxostat is a BCS class II xanthine-oxidase inhibitor used for hyperuricaemia.
Its plasma profiles frequently show a second absorption peak, which Chen and
co-authors set out to characterise explicitly. After screening Weibull,
sequential zero-to-first-order, simultaneous zero-and-first-order, variable
exponential, variable binomial and double-gamma absorption models
(Supplementary Information section 2), they selected a **dual-depot parallel
first-order absorption model**: fraction `F1` of the dose is absorbed from a
first depot at rate `ka1` after a lag `Tlag1`, and the remaining `1 - F1` from a
second depot at rate `ka2` after a lag `Tlag2`. Both depots feed a
two-compartment disposition model with first-order elimination. The final model
improved the objective function by 187 points over the base model
(OFV 30082.76 to 29895.08, Supplementary Information section 2).

# Population

The model was built on 2455 plasma concentrations from 128 healthy Chinese adult
volunteers (81 male, 47 female) pooled from two open-label, single-dose,
randomised, two-period crossover bioequivalence studies conducted at the Second
Affiliated Hospital of Zhejiang University School of Medicine (Hangzhou). Only
the reference-product periods entered the popPK analysis. Subjects received a
single 20 mg or 80 mg febuxostat tablet with 240 mL of water, either fasted (at
least 10 h pre-dose fast, first meal permitted 4 h post-dose) or with a
standardised high-fat, high-calorie breakfast (800-1000 kcal) consumed within
30 min of dosing.

Per Table 1 of the source, the four arms were 20 mg fasting (n = 36), 20 mg fed
(n = 31), 80 mg fasting (n = 31) and 80 mg fed (n = 30); median age ranged from
26.0 to 29.0 years (overall range 18-44), median body weight from 59.5 to
63.4 kg (overall range 45.4-79.0), and median BMI from 21.3 to 22.8 kg/m^2.

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

```{r population}
str(mod()$population, max.level = 1)
```

# Source trace

Every `ini()` entry in
`inst/modeldb/specificDrugs/Chen_2024_febuxostat.R` carries an in-file comment
pointing at its origin. They are collected here for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` (ka1) | 4.02 1/h | Table 3, "ka1 (h-1)", RSE 11.4% |
| `lka2` (ka2) | 4.34 1/h | Table 3, "ka2 (h-1)", RSE 15.6% |
| `lfdepot` (F1) | 0.51 | Table 3, "F1", RSE 3.7% |
| `ltlag` (Tlag1) | 0.22 h | Table 3, "Tlag1 (h)", RSE 8.93% |
| `ltlag2` (Tlag2) | 0.92 h | Table 3, "Tlag2 (h)", RSE 7.63% |
| `lcl` (CL/F) | 6.14 L/h | Table 3, "CL/F (L.h-1)", RSE 2.26% |
| `lvc` (Vc/F) | 10.38 L | Table 3, "Vc/F (L)", RSE 2.07% |
| `lq` (Q/F) | 1.93 L/h | Table 3, "Q/F (L.h-1)", RSE 4.25% |
| `lvp` (Vp/F) | 10.22 L | Table 3, "Vp/F (L)", RSE 3.21% |
| `e_fed_highfat_ka` | -2.04 | Table 3, "Food_ka1", RSE 8% |
| `e_fed_highfat_ka2` | -1.08 | Table 3, "Food_ka2", RSE 19.9% |
| `e_fed_highfat_tlag` | 0.73 | Table 3, "Food_Tlag1", RSE 18.2% |
| `e_wt_vc` | 0.58 | Table 3, "WT_Vc/F", RSE 27.9% |
| `e_wt_vp` | 1.15 | Table 3, "WT_Vp/F", RSE 22.5% |
| `addSd` | 2.14 ng/mL | Table 3, "Add. (ng/mL)", RSE 5.82% |
| `propSd` | 0.12 | Table 3, "Prop.", RSE 2.27% |
| All nine `eta*` variances | CV% column | Table 3, "IIV (CV%) (RSE%)"; converted as omega^2 = log(CV^2 + 1) |
| Categorical covariate form | n/a | Eq. (1): `log(P) = log(P_TV) + theta_COV * Food` |
| Continuous covariate form | n/a | Eq. (2): `log(P) = log(P_TV) + theta_COV * log(COV_i / COV_m)` |
| Dual-depot structure, dose split F1 / (1 - F1) | n/a | Results, "PopPK analysis" paragraph 1, and Fig. 2 |
| Two-compartment disposition, first-order elimination | n/a | Results, "PopPK analysis"; Fig. 2 |

# Covariate model checks

Two arithmetic identities in the source pin the functional form of the
covariate model without any simulation. The Discussion states that "food
reduced ka1 and ka2 by 87% and 66%, respectively". Under the exponential form
of Eq. (1) those reductions are `1 - exp(theta)`:

```{r covariate-form-gate}
theta_ka1  <- -2.04   # Table 3 Food_ka1
theta_ka2  <- -1.08   # Table 3 Food_ka2
theta_tlag <-  0.73   # Table 3 Food_Tlag1

reduction_ka1 <- 1 - exp(theta_ka1)
reduction_ka2 <- 1 - exp(theta_ka2)

c(ka1_reduction = reduction_ka1, ka2_reduction = reduction_ka2)

# The paper's own Discussion figures, to two significant figures.
stopifnot(
  abs(reduction_ka1 - 0.87) < 0.005,
  abs(reduction_ka2 - 0.66) < 0.005
)
```

Both reproduce the published percentages exactly, which confirms the
exponential (log-additive) form of Eq. (1) and the sign convention (`Food = 1`
in the fed state). The fed-state typical values implied by the model are
therefore:

```{r fed-typical-values}
tibble::tibble(
  Parameter = c("ka1 (1/h)", "ka2 (1/h)", "Tlag1 (h)", "Tlag2 (h)"),
  Fasted    = c(4.02, 4.34, 0.22, 0.92),
  Fed       = c(4.02 * exp(theta_ka1), 4.34 * exp(theta_ka2),
                0.22 * exp(theta_tlag), 0.92)
) |>
  knitr::kable(digits = 3, caption = "Typical absorption parameters by prandial state.")
```

# Virtual cohort

The original concentration data are not publicly available. The cohort below
reproduces the four study arms of Table 1 with 200 subjects each, using each
arm's published median and range for body weight and the **actual sampling
schedules** given in Supplementary Information section 1. Matching the real
sampling grid matters: the published NCA in Table 2 was computed on those
times, so a simulated NCA on the same grid carries the same trapezoidal and
Tmax-resolution bias and is directly comparable.

```{r cohort}
# Sampling schedules, Supplementary Information section 1.
times_20_fasting <- c(0, 0.17, 0.33, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25,
                      2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 24, 36, 48)
times_20_fed     <- c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25,
                      2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 24, 36, 48)
times_80_fasting <- c(0, 0.17, 0.33, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25,
                      2.5, 3, 3.5, 4, 6, 8, 10, 12, 24, 36)
times_80_fed     <- c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25,
                      2.5, 3, 3.5, 4, 6, 8, 10, 12, 24, 36)

n_per_arm <- 200L

# Deterministic weight cohort: evenly spaced normal quantiles centred on the
# arm's published median, with SD approximated as range / 4, then clamped to
# the published range. Deterministic quantiles (rather than a random draw)
# keep the arm-to-arm contrasts free of sampling noise.
weight_quantiles <- function(n, med, lo, hi) {
  w <- stats::qnorm(stats::ppoints(n), mean = med, sd = (hi - lo) / 4)
  pmin(pmax(w, lo), hi)
}

make_arm <- function(label, dose_mg, fed, times, med, lo, hi, id_offset) {
  subj <- tibble::tibble(
    id   = id_offset + seq_len(n_per_arm),
    WT   = weight_quantiles(n_per_arm, med, lo, hi),
    FED_HIGHFAT = fed,
    treatment   = label,
    dose_mg     = dose_mg
  )
  # The same dose record goes to both depots; f(depot) = F1 and
  # f(depot2) = 1 - F1 partition it inside the model.
  doses <- subj |>
    tidyr::crossing(cmt = c("depot", "depot2")) |>
    dplyr::mutate(time = 0, amt = dose_mg, evid = 1L)
  obs <- subj |>
    tidyr::crossing(time = times) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = "central")
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(
  make_arm("20 mg fasting", 20, 0, times_20_fasting, 61.7, 45.4, 73.6,   0L),
  make_arm("20 mg fed",     20, 1, times_20_fed,     63.4, 46.3, 79.0, 200L),
  make_arm("80 mg fasting", 80, 0, times_80_fasting, 59.5, 47.1, 73.7, 400L),
  make_arm("80 mg fed",     80, 1, times_80_fed,     60.4, 49.0, 78.4, 600L)
)

# Disjoint IDs across arms (a shared id silently merges arms into one subject).
stopifnot(
  dplyr::n_distinct(events$id) == 4L * n_per_arm,
  events |> dplyr::distinct(id, treatment) |> nrow() == 4L * n_per_arm
)
```

# Simulation

```{r simulate}
set.seed(20240824)
sim <- rxode2::rxSolve(
  mod, events,
  keep = c("treatment", "dose_mg", "WT", "FED_HIGHFAT")
) |>
  as.data.frame() |>
  dplyr::select(id, time, treatment, dose_mg, WT, FED_HIGHFAT, Cc, sim, fdepot)
```

`Cc` is the individual prediction; `sim` additionally carries the combined
residual error and is therefore the observed-concentration analogue used for
the NCA comparison below.

The lognormal `F1` distribution reported by the paper is unbounded above, so a
small fraction of simulated subjects draw `F1 > 1` and receive a negative dose
fraction into the second depot. The count is reported here rather than filtered
away, because filtering would misrepresent the published parameterisation:

```{r f1-tail}
f1_tail <- sim |>
  dplyr::distinct(id, fdepot) |>
  dplyr::summarise(n = dplyr::n(), n_gt1 = sum(fdepot > 1), pct = 100 * mean(fdepot > 1))
f1_tail
```

# Replicate published figures

```{r figure-1, fig.cap = "Replicates Figure 1 of Chen 2024."}
sim |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(mean_c = mean(sim), sd_c = stats::sd(sim), .groups = "drop") |>
  dplyr::mutate(
    prandial = ifelse(grepl("fasting", treatment), "(A) Fasted", "(B) Fed"),
    dose     = ifelse(grepl("^20", treatment), "20 mg", "80 mg")
  ) |>
  ggplot(aes(time, mean_c, colour = dose, fill = dose)) +
  geom_ribbon(aes(ymin = pmax(mean_c - sd_c, 1), ymax = mean_c + sd_c),
              alpha = 0.2, colour = NA) +
  geom_line() +
  facet_wrap(~prandial) +
  scale_y_log10() +
  coord_cartesian(xlim = c(0, 48)) +
  labs(x = "Time (h)", y = "Febuxostat concentration (ng/mL)",
       colour = "Dose", fill = "Dose",
       title = "Figure 1 - mean (SD) plasma profiles by dose and prandial state",
       caption = "Replicates Figure 1 of Chen 2024.")
```

The early-time detail is where the dual-depot structure does its work. Plotted
on a linear scale over the first 6 h, the typical-value profile shows the
shoulder that motivated the model:

```{r figure-1-early, fig.cap = "Early-time typical-value profiles showing the second absorption peak."}
# One typical subject per scenario; no IIV, so a single ID per curve suffices.
make_typical <- function(label, dose_mg, fed, wt, times, id) {
  subj <- tibble::tibble(id = id, WT = wt, FED_HIGHFAT = fed,
                         treatment = label, dose_mg = dose_mg)
  doses <- subj |>
    tidyr::crossing(cmt = c("depot", "depot2")) |>
    dplyr::mutate(time = 0, amt = dose_mg, evid = 1L)
  obs <- subj |>
    tidyr::crossing(time = times) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = "central")
  dplyr::bind_rows(doses, obs) |> dplyr::arrange(id, time, dplyr::desc(evid))
}

typ_events <- dplyr::bind_rows(
  make_typical("20 mg fasting", 20, 0, 60, seq(0, 12, by = 0.02), 1L),
  make_typical("20 mg fed",     20, 1, 60, seq(0, 12, by = 0.02), 2L)
)

typ <- rxode2::rxSolve(mod, typ_events, omega = NA, keep = "treatment") |>
  as.data.frame()

typ |>
  ggplot(aes(time, Cc, colour = treatment)) +
  geom_line(linewidth = 0.8) +
  coord_cartesian(xlim = c(0, 6)) +
  labs(x = "Time (h)", y = "Febuxostat concentration (ng/mL)", colour = NULL,
       title = "Typical-value 20 mg profiles - dual-depot shoulder",
       caption = "Structure of Figure 2 of Chen 2024; compare Figure 3(A).")
```

```{r figure-5, fig.cap = "Replicates Figure 5 of Chen 2024."}
# (A) dose ranging at 60 kg; (B) weight ranging at 60 mg.
fig5_times <- seq(0, 48, by = 0.05)

fig5a_events <- dplyr::bind_rows(lapply(seq_along(c(20, 40, 80, 120)), function(i) {
  d <- c(20, 40, 80, 120)[i]
  make_typical(paste0(d, " mg"), d, 0, 60, fig5_times, id = i)
}))

fig5a <- rxode2::rxSolve(mod, fig5a_events, omega = NA, keep = "treatment") |>
  as.data.frame() |>
  dplyr::mutate(panel = "(A) 20-120 mg, WT = 60 kg")

fig5b_events <- dplyr::bind_rows(lapply(seq_along(c(50, 60, 70, 80)), function(i) {
  w <- c(50, 60, 70, 80)[i]
  make_typical(paste0(w, " kg"), 60, 0, w, fig5_times, id = i)
}))

fig5b <- rxode2::rxSolve(mod, fig5b_events, omega = NA, keep = "treatment") |>
  as.data.frame() |>
  dplyr::mutate(panel = "(B) 60 mg, WT = 50-80 kg")

dplyr::bind_rows(fig5a, fig5b) |>
  ggplot(aes(time, Cc, colour = treatment)) +
  geom_line() +
  facet_wrap(~panel, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Febuxostat concentration (ng/mL)", colour = NULL,
       title = "Figure 5 - typical-value simulations in the Chinese population",
       caption = "Replicates Figure 5 of Chen 2024.")
```

Figure 5(B) reproduces the paper's observation that heavier subjects have
slightly lower exposure early on: weight enters only the two apparent volumes,
so peak concentration falls with weight while `AUC0-inf = Dose / (CL/F)` is
weight-invariant. Because total bioavailability is 1 by construction (the two
depots receive `F1` and `1 - F1` of the same dose), that identity is exact and
gives a free end-to-end gate on the whole ODE system.

```{r weight-invariance-gate}
auc_closed_form <- 1000 * 60 / 6.14   # 1000 * Dose(mg) / (CL/F, L/h) -> ng*h/mL

auc_by_weight <- fig5b |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    auc  = sum(diff(time) * (utils::head(Cc, -1) + utils::tail(Cc, -1)) / 2),
    cmax = max(Cc),
    .groups = "drop"
  ) |>
  dplyr::mutate(pct_of_closed_form = 100 * auc / auc_closed_form)
auc_by_weight

# Every weight's AUC0-48 reproduces Dose/(CL/F) to better than 0.5% (the small
# residual is the un-integrated terminal tail beyond 48 h, which grows with
# weight), while Cmax varies by nearly 18% across the same weight span.
stopifnot(
  all(abs(auc_by_weight$auc - auc_closed_form) / auc_closed_form < 0.005),
  diff(range(auc_by_weight$cmax)) / mean(auc_by_weight$cmax) > 0.05
)
```

Linearity is a second free gate: the model has no saturable term, so
typical-value exposure must be exactly dose-proportional.

```{r dose-linearity-gate}
auc_by_dose <- fig5a |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    auc = sum(diff(time) * (utils::head(Cc, -1) + utils::tail(Cc, -1)) / 2),
    .groups = "drop"
  ) |>
  dplyr::mutate(dose = as.numeric(sub(" mg", "", treatment)),
                auc_per_mg = auc / dose)
auc_by_dose

stopifnot(diff(range(auc_by_dose$auc_per_mg)) / mean(auc_by_dose$auc_per_mg) < 1e-6)
```

# PKNCA validation

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

# Guarantee a time = 0 row per subject (extravascular pre-dose Cc = 0).
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment, dose_mg) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)

dose_df <- events |>
  dplyr::filter(evid == 1, cmt == "depot") |>
  dplyr::select(id, time, amt, treatment)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

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

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

## Comparison against published NCA

Table 2 of the source reports Cmax and AUC as arithmetic **means** and Tmax as a
**median**. `ncaComparisonTable()` collapses multiple rows per group with the
median, so the simulated values are pre-aggregated here with the matching
statistic per parameter.

```{r nca-aggregate}
nca_wide <- as.data.frame(nca_res) |>
  dplyr::filter(start == 0, end == Inf) |>
  dplyr::select(treatment, id, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

simulated <- nca_wide |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    cmax       = mean(cmax, na.rm = TRUE),
    auclast    = mean(auclast, na.rm = TRUE),
    aucinf.obs = mean(aucinf.obs, na.rm = TRUE),
    tmax       = stats::median(tmax, na.rm = TRUE),
    .groups = "drop"
  )

published <- tibble::tribble(
  ~treatment,      ~cmax, ~tmax, ~auclast,  ~aucinf.obs,
  "20 mg fasting", 1091,  1.25,  3499.48,   3552.00,
  "20 mg fed",      856.0, 1.50, 2926.99,   2984.70,
  "80 mg fasting", 4327,  1.00,  15437.21,  15549.11,
  "80 mg fed",     3432,  1.75,  13972.88,  14086.00
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = simulated,
  reference = published,
  by        = "treatment",
  units     = c(cmax = "ng/mL", auclast = "h*ng/mL",
                aucinf.obs = "h*ng/mL", tmax = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste(
    "Simulated vs. published NCA (Chen 2024 Table 2).",
    "Reference Cmax / AUC are arithmetic means, Tmax is a median.",
    "* marks a difference greater than 20%."
  )
)
```

All twelve Cmax and AUC comparisons agree to within 15% except the two 20 mg
fed AUC rows, and four rows exceed the 20% flag: `Tmax` in the 20 mg fed
(+50%) and 80 mg fasting (+25%) arms, and `AUC0-inf` / `AUClast` in the 20 mg
fed arm (+24%). Each traces to a property of the published model rather than to
the transcription, and no parameter was adjusted to improve any comparison.

**The 20 mg fed AUC stars are a deliberate feature of the published model.**
Stepwise covariate modelling retained food only on `ka1`, `ka2` and `Tlag1` --
never on `CL/F` and never on bioavailability -- so the model predicts
*identical* AUC fasted and fed by construction, whereas the observed 20 mg data
showed a significant food effect on AUC (p = 0.0021 in Table 2). That the
mismatch is structural rather than a mis-transcription can be checked
mechanically: the simulated fasted and fed AUCs must coincide at each dose.

```{r food-auc-identity-gate}
auc_food <- simulated |>
  dplyr::mutate(dose = sub(" (fasting|fed)$", "", treatment),
                state = ifelse(grepl("fasting$", treatment), "fasted", "fed")) |>
  dplyr::select(dose, state, aucinf.obs) |>
  tidyr::pivot_wider(names_from = state, values_from = aucinf.obs) |>
  dplyr::mutate(
    simulated_fed_fasted_ratio = fed / fasted,
    published_fed_fasted_ratio = c(2984.70 / 3552.00, 14086.00 / 15549.11)
  )
auc_food

# The model carries no food effect on exposure, so the simulated ratio is 1.
stopifnot(all(abs(auc_food$simulated_fed_fasted_ratio - 1) < 0.02))
```

The simulated fed/fasted AUC ratio is 1 at both doses, while the published
observed ratio is 0.84 at 20 mg and 0.91 at 80 mg. The 20 mg fed shortfall in
the comparison table is exactly that missing 16%; the 80 mg gap is smaller and
the paper itself reports the 80 mg AUC food effect as non-significant
(p = 0.14).

**The remaining AUC residuals are a dose-linearity compromise.** The observed
data are mildly supra-proportional -- observed AUC0-inf per mg is
178 h*ng/mL at 20 mg fasting versus 194 at 80 mg fasting -- whereas the model
is strictly linear. A single linear model fitted to both dose levels splits the
difference, so the 20 mg fasting arm sits above the published value (+3.6%) and
the 80 mg fasting arm below it (-14.6%).

**The Tmax stars are small absolute differences amplified by a ratio metric.**
They are 0.75 h (20 mg fed, 1.50 h to 2.25 h) and 0.25 h (80 mg fasting, 1.00 h
to 1.25 h) -- one to three steps on a sampling grid whose spacing is 0.25 h in
that window. With CV values above 100% on both absorption rate constants and on
both lag times, the median Tmax is a coarse summary. Note also that the 80 mg
fasting Tmax entry in Table 2 prints its range as "(4.00, 0.33)", with minimum
and maximum reversed relative to every other row of that table.

**Fed Cmax is under-predicted by about 12% at both doses.** This has the same
sign as the paper's own finding that the model "underpredicted" fasted-state
absorption in its Caucasian simulation, and reflects the very large `ka`
inter-individual variability interacting with a mean-of-maxima summary.

# Exploratory Caucasian scenario (Supplementary Table 1)

The paper's final section takes the Chinese model, adjusts it for reported
race differences, and compares against the only Caucasian food-effect study it
could find (Khosravan 2008, 40 mg and 120 mg fasted and fed, n = 24). The stated
adjustments are that "the bioavailability was decreased by 46%, and CL/F was
increased to 6.9 L/h"; all other parameters were unchanged.

Applied literally those two numbers do **not** reproduce the paper's own
Supplementary Table 1 predictions, so the relative bioavailability is treated
here as an unknown and back-solved from the published predicted AUC0-inf. That
leaves Cmax and Tmax as genuinely independent checks of the absorption model,
since neither was used in the back-solve.

```{r caucasian-setup}
mod_cauc <- rxode2::ini(rxode2::rxode2(mod), lcl = log(6.9))

khosravan_times <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 8,
                     10, 12, 16, 24, 36, 48)

make_cauc <- function(label, dose_mg, fed, id_offset, frel = 1) {
  subj <- tibble::tibble(
    id = id_offset + seq_len(24L), WT = 70, FED_HIGHFAT = fed,
    treatment = label, dose_mg = dose_mg
  )
  doses <- subj |>
    tidyr::crossing(cmt = c("depot", "depot2")) |>
    dplyr::mutate(time = 0, amt = frel * dose_mg, evid = 1L)
  obs <- subj |>
    tidyr::crossing(time = khosravan_times) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = "central")
  dplyr::bind_rows(doses, obs) |> dplyr::arrange(id, time, dplyr::desc(evid))
}

cauc_events <- function(frel) {
  dplyr::bind_rows(
    make_cauc("40 mg fasting",  40, 0,    0L, frel),
    make_cauc("40 mg fed",      40, 1,   50L, frel),
    make_cauc("120 mg fasting",120, 0,  100L, frel),
    make_cauc("120 mg fed",    120, 1,  150L, frel)
  )
}

cauc_nca <- function(frel, seed = 424242) {
  set.seed(seed)
  s <- rxode2::rxSolve(mod_cauc, cauc_events(frel),
                       keep = c("treatment", "dose_mg")) |>
    as.data.frame()
  d <- s |>
    dplyr::filter(!is.na(sim)) |>
    dplyr::select(id, time, treatment, Cc = sim)
  d <- dplyr::bind_rows(
    d, d |> dplyr::distinct(id, treatment) |> dplyr::mutate(time = 0, Cc = 0)
  ) |>
    dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
    dplyr::arrange(id, time)
  dose_d <- cauc_events(frel) |>
    dplyr::filter(evid == 1, cmt == "depot") |>
    dplyr::select(id, time, amt, treatment)
  r <- PKNCA::pk.nca(PKNCA::PKNCAdata(
    PKNCA::PKNCAconc(d, Cc ~ time | treatment + id),
    PKNCA::PKNCAdose(dose_d, amt ~ time | treatment + id),
    intervals = data.frame(start = 0, end = Inf,
                           cmax = TRUE, tmax = TRUE, aucinf.obs = TRUE)
  ))
  as.data.frame(r) |>
    dplyr::filter(start == 0, end == Inf) |>
    dplyr::select(treatment, id, PPTESTCD, PPORRES) |>
    tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
    dplyr::group_by(treatment) |>
    dplyr::summarise(cmax = mean(cmax, na.rm = TRUE),
                     aucinf.obs = mean(aucinf.obs, na.rm = TRUE),
                     tmax = mean(tmax, na.rm = TRUE), .groups = "drop")
}
```

```{r caucasian-backsolve}
# Supplementary Table 1, "Predicted" columns.
supp_predicted <- tibble::tribble(
  ~treatment,       ~cmax, ~aucinf.obs, ~tmax,
  "40 mg fasting",  1376,  4364,        1.5,
  "40 mg fed",      1005,  4397,        2.3,
  "120 mg fasting", 4667,  13801,       1.6,
  "120 mg fed",     3207,  14020,       2.4
)

unit_frel <- cauc_nca(frel = 1)

frel_by_arm <- supp_predicted |>
  dplyr::select(treatment, published = aucinf.obs) |>
  dplyr::inner_join(dplyr::select(unit_frel, treatment, unit_F = aucinf.obs),
                    by = "treatment") |>
  dplyr::mutate(implied_relative_F = published / unit_F)
frel_by_arm

frel_hat <- mean(frel_by_arm$implied_relative_F)
stopifnot(is.finite(frel_hat), frel_hat > 0, frel_hat < 2)

c(back_solved_relative_F = frel_hat,
  literal_reading_1_minus_0.46 = 0.54)
```

The AUC-implied relative bioavailability is about `r round(frel_hat, 2)`, not
the 0.54 that "decreased by 46%" implies. With that value, Cmax and Tmax
(neither used in the back-solve) can be checked against the paper's own
predictions:

```{r caucasian-compare}
cauc_cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = cauc_nca(frel = frel_hat),
  reference = supp_predicted,
  by        = "treatment",
  units     = c(cmax = "ng/mL", aucinf.obs = "h*ng/mL", tmax = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cauc_cmp,
  caption = paste(
    "Caucasian scenario vs. Chen 2024 Supplementary Table 1 'Predicted' columns.",
    "Relative bioavailability back-solved from AUC0-inf, so only the Cmax and",
    "Tmax rows are independent checks. * marks a difference greater than 20%."
  )
)
```

Cmax agrees to within 9% in all four arms and Tmax to within 16%, with the
fed-versus-fasted Tmax ordering reproduced in both dose groups. Because the
back-solve consumed only the AUC information, this is a genuine independent
check on the absorption sub-model against a target -- the supplement's
Caucasian simulation -- that shares no numbers with Table 2 or Table 3.

```{r caucasian-gate}
cauc_cmax_diff <- cauc_cmp |>
  dplyr::filter(grepl("^Cmax", .data[["NCA parameter"]])) |>
  dplyr::pull(`% diff`) |>
  sub(pattern = "[*]$", replacement = "") |>
  sub(pattern = "%$", replacement = "") |>
  as.numeric()
cauc_cmax_diff

stopifnot(all(abs(cauc_cmax_diff) < 15))
```

It follows that the absorption sub-model is transcribed correctly and that the
only irreconcilable quantity is the single scalar bioavailability adjustment,
which the paper under-specifies.

# Assumptions and deviations

* **Weight-centering value (`WT_ref = 60 kg`) is assumed.** Eq. (2) centres
  continuous covariates on `COV_m`, "the median value of the covariates", but
  the pooled median weight is never printed. Table 1 gives per-arm medians of
  61.7, 63.4, 59.5 and 60.4 kg (n-weighted mean 61.3 kg), and the authors'
  own typical-subject simulations in Figure 5 use 60 kg, so 60 kg is used here.
  The choice is nearly inconsequential: at `WT = 60` kg versus a 61.3 kg
  reference the implied `Vc/F` differs by 1.2% and `AUC` not at all.
* **IIV column read as a coefficient of variation.** Table 3's "IIV (CV%)"
  values are converted with `omega^2 = log(CV^2 + 1)`, Monolix's lognormal CV
  convention. The alternative reading (the column being `100 * omega`) was
  tested against the observed Cmax variability in Table 2 (27.6-33.4% CV across
  the four arms): the CV conversion gives 27.9% (fasted) / 35.5% (fed), while
  the `100 * omega` reading gives 31.9% / 41.7%. The CV conversion is the better
  match in every arm and is also what Monolix reports natively for
  lognormally-distributed parameters.
* **`F1` is encoded with the lognormal IIV the paper states, and is therefore
  unbounded above.** `F1 = 0.51` with a 32.87% CV puts 1.8% of subjects above
  `F1 = 1`, where `f(depot2) = 1 - F1` turns negative; about 0.4% of subjects
  show negative concentrations as a result (see the `f1-tail` chunk, which
  reports the realised count for this cohort). A logit-normal distribution
  would be the usual Monolix choice for a bounded fraction, but the paper
  neither states one nor reports a logit-scale omega, so no substitute value
  was invented. Typical-value simulations (`omega = NA`) are unaffected.
* **Combined residual-error variance form.** The paper reports an additive term
  (2.14 ng/mL) and a proportional term (0.12) but does not say whether Monolix's
  `combined1` (`a + b*f`) or `combined2` (`sqrt(a^2 + (b*f)^2)`) form was used.
  The rxode2 default (`combined2`) is used. The additive term is three orders of
  magnitude below the observed concentration range, so the choice has no
  practical effect.
* **`FED_HIGHFAT` rather than the generic `FED`.** The paper's covariate is
  called "Food", and its fed arm is specifically the FDA-style 800-1000 kcal
  high-fat breakfast described in Supplementary Information section 1, which
  matches the `FED_HIGHFAT` register entry exactly. Orientation and reference
  category are unchanged (1 = fed, 0 = fasted).
* **Virtual-cohort weight distributions** are deterministic normal quantiles
  centred on each arm's published median with SD approximated as range / 4,
  clamped to the published range. The paper reports only medians and ranges.
* **Caucasian scenario body weight assumed to be 70 kg.** Chen 2024 states that
  the simulated demographics matched Khosravan 2008 but does not reprint them.
  Weight affects only the apparent volumes, so this assumption moves Cmax but
  leaves AUC0-inf untouched, and the back-solved relative bioavailability is
  therefore unaffected by it.
* **Caucasian relative bioavailability back-solved rather than taken as 0.54.**
  See the section above. The model file itself contains no Caucasian
  adjustment; it encodes only the final Chinese model of Table 3.
* **Sampling schedule for the Caucasian scenario** is a generic dense grid;
  Chen 2024 does not reprint the Khosravan 2008 sampling times.
* **Not carried into the model file**: sex, age and BMI were screened by
  stepwise covariate modelling and not retained in the final model (recorded in
  the model's `covariatesDataExcluded` metadata); the base-model and
  intermediate absorption structures of Supplementary Information section 2 are
  model-development steps rather than final results.
