---
title: "Cetagliptin (Zhou 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cetagliptin (Zhou 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

* Citation: `r rxode2::rxode(readModelDb("Zhou_2024_cetagliptin"))$reference`
* Description: `r rxode2::rxode(readModelDb("Zhou_2024_cetagliptin"))$description`
* Article: <https://doi.org/10.3389/fendo.2024.1359407>
* Supplement (Supplementary Table 1 -- the only source of the final parameter
  estimates): <https://www.frontiersin.org/articles/10.3389/fendo.2024.1359407/full#supplementary-material>

Cetagliptin (CAS 2243737-33-7) is a dipeptidyl peptidase-4 (DPP-4) inhibitor
under development for type 2 diabetes mellitus (T2DM). Zhou 2024 reports a
sequential two-step population PK/PD analysis performed in Phoenix NLME: a
population PK model was fit first, its parameters were fixed, and a
direct-effect sigmoid Emax model was then fit linking cetagliptin plasma
concentration to plasma DPP-4 inhibition.

The structural PK model is two-compartment with first-order oral absorption and
**saturable Michaelis-Menten elimination**. The main text describes only "the
two-compartment model"; the Vmax/Km parameterisation is visible solely in
Supplementary Table 1, which reports `tvKm` and `tvVmax` and no clearance term.
Saturation is material at therapeutic exposures: Km is 171.5 ng/mL while the
observed steady-state concentrations span roughly 50-300 ng/mL.

# Population

Thirty-two Chinese adults with T2DM were enrolled at a single centre (the First
Affiliated Hospital with Nanjing Medical University; CTR20190599) into two dose
groups of 16. Within each group subjects were randomised 10:2:4 to cetagliptin
(50 or 100 mg), matching placebo, or open-label sitagliptin 100 mg, so 20
subjects received cetagliptin (10 per dose level) and contribute the cetagliptin
concentrations that the population PK model was fit to. Dosing was oral, once
daily under fasting conditions, for 14 consecutive days.

Baseline characteristics (Table 1) were comparable across arms: mean age 47.8
(SD 4.3) and 45.2 (SD 9.8) years, mean weight 69.6 (SD 7.6) and 72.6 (SD 7.6) kg,
and mean BMI 25.5 (SD 1.9) and 25.8 (SD 2.3) kg/m^2 in the 50 mg and 100 mg
cetagliptin arms respectively. Protocol inclusion required age 18-65 years, BMI
19.00-30.00 kg/m^2, HbA1c from 6.5% to below 9%, and fasting blood glucose below
13.4 mmol/L. Baseline HbA1c was 8.21 (SD 0.66)% and 7.79 (SD 0.53)%, and baseline
fasting plasma glucose 7.96 (SD 1.29) and 6.87 (SD 1.40) mmol/L. Four of the 20
cetagliptin recipients were female (90% and 70% male in the 50 mg and 100 mg
arms). No subject had clinically significant abnormal liver function.

The same information is available programmatically via the model's `population`
metadata (`readModelDb("Zhou_2024_cetagliptin")()$population`).

# Source trace

Every final parameter estimate comes from **Supplementary Table 1**
("Parameter estimates and bootstrap results of the final population
pharmacokinetic model"), which is the only place in the publication where the
model parameters appear. The main article text reports no parameter values.
The supplement expresses amounts in ug and concentrations in ug/L (= ng/mL),
so the values below are transcribed verbatim; the mg-to-ug dose conversion is
applied once in `model()` via `f(depot) <- 1000`.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` | 0.06521 1/h | Supplementary Table 1, `tvKa` (RSE 1.623%) |
| `lvc` | 8.668 L | Supplementary Table 1, `tvV` (RSE 2.234%) |
| `lvp` | 558.0 L | Supplementary Table 1, `tvV2` (RSE 2.023%) |
| `lq` | 9.671 L/h | Supplementary Table 1, `tvCl2` (RSE 2.050%) |
| `lkm` | 171.5 ng/mL | Supplementary Table 1, `tvKm` (RSE 6.694%) |
| `lvmax` | 9373 ug/h | Supplementary Table 1, `tvVmax` (RSE 5.044%) |
| `e_tbili_vp` | 0.3723 | Supplementary Table 1, `dV2dTBIL` (RSE 3.066%) |
| `lemax` | 91.78 % | Supplementary Table 1, `tvEmax` (RSE 0.6484%) |
| `lec50` | 5.120 ng/mL | Supplementary Table 1, `tvEC50` (RSE 3.850%) |
| `lhill` | 1.008 | Supplementary Table 1, `tvGam` (RSE 4.296%) |
| `etalvc` | 0.8112 | Supplementary Table 1, `omega2 V` (shrinkage 18.24%) |
| `etalvmax` | 0.0260 | Supplementary Table 1, `omega2 Vmax` (shrinkage 8.642%) |
| `etalka` | 0.01462 | Supplementary Table 1, `omega2 Ka` (shrinkage 11.31%) |
| `etalvp` | 0.02874 | Supplementary Table 1, `omega2 V2` (shrinkage 12.68%) |
| `etalkm` | 0.01278 | Supplementary Table 1, `omega2 Km` (shrinkage 43.52%) |
| `etalec50` | 0.01259 | Supplementary Table 1, `omega2 EC50` (shrinkage 18.13%) |
| `etalemax` | 7.262e-05 | Supplementary Table 1, `omega2 Emax` (shrinkage 17.58%) |
| `propSd` | 0.2241 | Supplementary Table 1, "Multiplicative residual variability PK(sigma)", `Stdev0` |
| `addSd_dpp4Inhibition` | 11.62 % | Supplementary Table 1, "MixRatio residual variability PD(sigma)", `Stdev0` |
| Exponential IIV, `Pij = Pj * exp(eta_ij)` | n/a | Methods section 2.7.1, Equation 1 |
| Two-compartment structural PK model | n/a | Results section 3.6.1 |
| Michaelis-Menten elimination `Vmax * Cc / (Km + Cc)` | n/a | Supplementary Table 1 parameterisation (`tvVmax`, `tvKm`; no clearance term) |
| TBIL effect on V2, `V2 increases with TBIL` | n/a | Discussion page 10; Results section 3.6.1 |
| Direct-effect sigmoid Emax on DPP-4 inhibition | n/a | Methods section 2.7.2; Results section 3.6.2 |
| PK sampling schedule used for the NCA | see below | Methods section 2.3.1 |
| Reference NCA values (Table 2) | see below | Table 2, "Pharmacokinetic parameters after single and multiple oral doses" |
| Reference PD values (Table 3) | see below | Table 3, "Pharmacodynamic parameters of DPP-4 inhibition" |

# Virtual cohort

Original observed data are not publicly available. The simulations below use
virtual populations matching the published trial design: two arms at the
protocol dose levels, 50 mg and 100 mg once daily for 14 days.

Sampling deliberately reproduces **the study's own PK/PD schedule** (Methods
section 2.3.1) rather than a convenience grid, because every reference value in
Tables 2 and 3 is a non-compartmental estimate computed from those specific
times. Cmax read off a dense grid is biased upward relative to Cmax read off a
sparse clinical schedule, so matching the schedule is what makes the comparison
an apples-to-apples test of the model rather than of the sampling density. A
denser grid is added on top for the figures only, and is excluded from the NCA.

Total bilirubin is the model's only covariate, and the paper reports neither
its distribution nor its population median (TBIL does not appear in the Table 1
demographics). Because no subject had clinically significant abnormal liver
function, the virtual cohort draws TBILI from a log-normal distribution centred
on the model's assumed reference of 10 umol/L and truncated to the adult normal
range. See "Assumptions and deviations".

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

n_per_arm <- 100L
tbili_ref <- 10 # umol/L; assumed centring value, see Assumptions
tau <- 24       # dosing interval (h)
t_d14 <- 13 * tau # time of the day-14 (final) dose

# The study's PK / PD sampling schedule (Methods section 2.3.1): dense over
# day 1, troughs on days 7 and 10, and dense over day 14 with a 120 h washout
# tail. These are the ONLY times used for the NCA comparisons.
t_nca <- sort(unique(c(
  c(0, 0.5, 1, 2, 3, 4, 5, 6, 8, 12, 24),
  6 * tau, 9 * tau,
  t_d14 + c(0, 0.5, 1, 2, 3, 4, 5, 6, 8, 12, 24, 48, 72, 96, 120)
)))

# Extra times for smooth figures only.
t_fig <- sort(unique(c(
  seq(0, 24, by = 0.25),
  seq(24, t_d14, by = 6),
  seq(t_d14, t_d14 + 24, by = 0.25),
  seq(t_d14 + 24, t_d14 + 120, by = 4)
)))

t_obs <- sort(unique(c(t_nca, t_fig)))

make_arm <- function(n, dose_mg, label, id_offset = 0L) {
  subj <- tibble(
    id = id_offset + seq_len(n),
    TBILI = pmin(pmax(rlnorm(n, meanlog = log(tbili_ref), sdlog = 0.30), 3), 21),
    treatment = label,
    dose_mg = dose_mg
  )
  dosing <- subj |>
    tidyr::crossing(time = seq(0, t_d14, by = tau)) |>
    mutate(amt = dose_mg, evid = 1L, cmt = "depot", dvid = NA_integer_)
  obs <- subj |>
    tidyr::crossing(time = t_obs, dvid = c(1L, 2L)) |>
    mutate(amt = NA_real_, evid = 0L, cmt = NA_character_)
  bind_rows(dosing, obs) |>
    arrange(id, time, desc(evid))
}

events <- bind_rows(
  make_arm(n_per_arm,  50, "Cetagliptin 50 mg",  id_offset = 0L),
  make_arm(n_per_arm, 100, "Cetagliptin 100 mg", id_offset = n_per_arm)
)
```

This model declares two endpoints (`Cc` and `dpp4Inhibition`), so observation
rows are routed by DV id: `dvid = 1` selects `Cc` and `dvid = 2` selects
`dpp4Inhibition`. Because both endpoints are algebraic expressions of the ODE
states, **both** appear as columns on **every** returned row; the `dvid` column
is therefore carried through the solve and used to select the rows belonging to
each endpoint. Filtering on `!is.na(Cc)` instead would silently retain each
time point twice.

# Simulation

```{r simulate}
mod <- readModelDb("Zhou_2024_cetagliptin")

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("treatment", "TBILI", "dose_mg", "dvid")
) |>
  as.data.frame()

stopifnot(dplyr::n_distinct(sim$id) == 2L * n_per_arm)
stopifnot(all(c("Cc", "dpp4Inhibition", "ipredSim", "sim") %in% names(sim)))

# `Cc` / `dpp4Inhibition` are the individual predictions (IPRED); `sim` is the
# same quantity with the model's residual error added, i.e. the scale on which
# the paper's observed NCA values were measured. Split the two endpoints and
# carry both scales.
lloq_pk <- 0.5 # ng/mL, lower limit of the validated LC-MS/MS range (Methods 2.3.1)

pk <- sim |>
  filter(dvid == 1L) |>
  transmute(
    id, time, treatment, dose_mg, TBILI,
    ipred = Cc,
    observed = pmax(sim, 0) * (pmax(sim, 0) >= lloq_pk)
  )

pd <- sim |>
  filter(dvid == 2L) |>
  transmute(id, time, treatment, dose_mg, ipred = dpp4Inhibition, observed = sim)

stopifnot(nrow(pk) == nrow(pd), nrow(pk) == 2L * n_per_arm * length(t_obs))
```

The published figures are mean profiles, so the typical-value (random-effects
zeroed) prediction is also carried through for direct overlay. `omega = NA` is
used rather than `zeroRe()`, which mutates shared model state.

```{r simulate-typical}
events_typical <- events |>
  filter(id %in% c(1L, n_per_arm + 1L)) |>
  mutate(TBILI = tbili_ref)

sim_typical <- rxode2::rxSolve(
  mod, events = events_typical,
  keep = c("treatment", "TBILI", "dose_mg", "dvid"),
  omega = NA, sigma = NA
) |>
  as.data.frame() |>
  filter(dvid == 1L)
```

## Deterministic checks

Two properties follow from the model structure alone and hold regardless of the
simulated cohort, so they are asserted rather than merely plotted.

**Steady-state mass balance.** Under once-daily dosing the amount cleared over
one steady-state interval must equal the administered dose. With Michaelis-Menten
elimination that is a genuine constraint linking `tvVmax`, `tvKm` and the
observed exposure, and it is the check that confirms `tvVmax` is in ug/h (not
mg/h) and that `f(depot) <- 1000` converts the mg dose correctly.

```{r gate-mass-balance}
mass_balance <- sim_typical |>
  filter(time >= t_d14, time <= t_d14 + tau) |>
  arrange(dose_mg, time) |>
  group_by(dose_mg) |>
  summarise(
    # ug eliminated over the interval = integral of vmax * Cc / (km + Cc)
    eliminated_ug = {
      rate <- vmax * Cc / (km + Cc)
      sum(diff(time) * (head(rate, -1) + tail(rate, -1)) / 2)
    },
    dose_ug = unique(dose_mg) * 1000,
    .groups = "drop"
  ) |>
  mutate(pct_of_dose = 100 * eliminated_ug / dose_ug)

knitr::kable(
  mass_balance |>
    dplyr::rename(
      "Dose (mg)" = dose_mg, "Eliminated over tau (ug)" = eliminated_ug,
      "Dose (ug)" = dose_ug, "Percent of dose" = pct_of_dose
    ),
  digits = 1,
  caption = "Steady-state mass balance over the day-14 dosing interval."
)

# By day 14 (>= 7 elimination half-lives) accumulation is essentially complete,
# so the interval must close to within a couple of percent.
stopifnot(all(abs(mass_balance$pct_of_dose - 100) < 3))
```

**Saturable elimination is active.** Km (171.5 ng/mL) sits inside the observed
concentration range, so doubling the dose must raise steady-state AUC by more
than two-fold. A linear-clearance encoding would give exactly two-fold, so this
assertion distinguishes the published Vmax/Km parameterisation from the
"two-compartment model" the main text describes.

```{r gate-nonlinearity}
auc_ss <- sim_typical |>
  filter(time >= t_d14, time <= t_d14 + tau) |>
  arrange(dose_mg, time) |>
  group_by(dose_mg) |>
  summarise(
    auc = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups = "drop"
  )

dose_ratio <- max(auc_ss$dose_mg) / min(auc_ss$dose_mg)
auc_ratio <- max(auc_ss$auc) / min(auc_ss$auc)

cat(sprintf(
  "Dose ratio %.1f-fold gives a steady-state AUC ratio of %.2f-fold.\n",
  dose_ratio, auc_ratio
))
stopifnot(auc_ratio > dose_ratio * 1.05)
```

# Replicate published figures

```{r figure-1}
# Replicates Figure 1 of Zhou 2024: mean cetagliptin plasma concentration-time
# profiles after single (day 1) and multiple (day 14) oral doses.
day_panels <- function(df) {
  df |>
    mutate(
      panel = case_when(
        time <= tau ~ "Day 1 (single dose)",
        time >= t_d14 ~ "Day 14 (steady state)",
        TRUE ~ NA_character_
      ),
      tad = if_else(time >= t_d14, time - t_d14, time)
    ) |>
    filter(!is.na(panel), tad <= tau)
}

day_panels(pk) |>
  group_by(panel, treatment, tad) |>
  summarise(
    Q05 = quantile(ipred, 0.05), Mean = mean(ipred),
    Q95 = quantile(ipred, 0.95), .groups = "drop"
  ) |>
  ggplot(aes(tad, Mean, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~panel) +
  labs(
    x = "Time after dose (h)", y = "Cetagliptin concentration (ng/mL)",
    colour = NULL, fill = NULL,
    title = "Figure 1 - cetagliptin plasma concentration-time profiles",
    caption = "Replicates Figure 1 of Zhou 2024. Line = mean, band = 5th-95th percentile."
  ) +
  theme_bw() + theme(legend.position = "bottom")
```

```{r figure-2}
# Replicates Figure 2 of Zhou 2024: mean plasma DPP-4 inhibition-time profiles.
day_panels(pd) |>
  group_by(panel, treatment, tad) |>
  summarise(
    Q05 = quantile(ipred, 0.05), Mean = mean(ipred),
    Q95 = quantile(ipred, 0.95), .groups = "drop"
  ) |>
  ggplot(aes(tad, Mean, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = 80, linetype = "dashed", colour = "grey40") +
  facet_wrap(~panel) +
  labs(
    x = "Time after dose (h)", y = "DPP-4 inhibition (%)",
    colour = NULL, fill = NULL,
    title = "Figure 2 - plasma DPP-4 inhibition-time profiles",
    caption = paste(
      "Replicates Figure 2 of Zhou 2024. Dashed line = the paper's 80%",
      "inhibition threshold used for DUR80%."
    )
  ) +
  theme_bw() + theme(legend.position = "bottom")
```

The paper's `DUR80%` -- the time DPP-4 inhibition stays above 80% -- is the PD
descriptor that most directly separates the two dose levels, so it is
recomputed from the simulation and compared numerically.

```{r dur80}
dur80 <- pd |>
  filter(time >= t_d14, time <= t_d14 + tau) |>
  arrange(id, time) |>
  group_by(id, treatment) |>
  summarise(
    # Linear interpolation between adjacent samples, as WinNonlin would do.
    dur = {
      above <- ipred >= 80
      seg <- diff(time)
      w <- (head(above, -1) + tail(above, -1)) / 2
      sum(seg * w)
    },
    .groups = "drop"
  ) |>
  group_by(treatment) |>
  summarise(`Simulated DUR80% (h)` = mean(dur), .groups = "drop") |>
  mutate(`Published DUR80% (h)` = c(32.3, 21.9)[match(
    treatment, c("Cetagliptin 100 mg", "Cetagliptin 50 mg")
  )])

knitr::kable(
  dur80 |> dplyr::rename("Treatment" = treatment),
  digits = 1,
  caption = paste(
    "Duration of DPP-4 inhibition above 80% over the day-14 dosing interval,",
    "against Zhou 2024 Table 3 (DUR80%, day 14). The published values exceed",
    "the 24 h interval because they were computed over the full day-14",
    "profile including the 120 h washout; the simulated values are censored",
    "at tau = 24 h and so are bounded above by 24."
  )
)
```

```{r figure-3}
# Replicates Figure 3 of Zhou 2024: the DPP-4 inhibition vs. cetagliptin
# concentration Emax relationship.
emax_curve <- tibble(
  Cc = 10^seq(-1, 3, length.out = 200)
) |>
  mutate(dpp4Inhibition = 91.78 * Cc^1.008 / (5.120^1.008 + Cc^1.008))

tibble(Cc = pk$ipred, dpp4Inhibition = pd$ipred) |>
  filter(Cc > 0) |>
  ggplot(aes(Cc, dpp4Inhibition)) +
  geom_point(alpha = 0.05, size = 0.5, colour = "steelblue") +
  geom_line(data = emax_curve, colour = "black", linewidth = 0.9) +
  scale_x_log10() +
  labs(
    x = "Cetagliptin plasma concentration (ng/mL)", y = "DPP-4 inhibition (%)",
    title = "Figure 3 - concentration vs. DPP-4 inhibition",
    caption = paste(
      "Replicates Figure 3 of Zhou 2024. Points = virtual subjects,",
      "line = typical-value sigmoid Emax curve."
    )
  ) +
  theme_bw()
```

The paper's own concentration-effect fit of the observed data (Results section
3.4.1, a non-compartmental Phoenix WinNonlin Emax fit that is distinct from the
population model) gave Emax 92.47% and EC50 5.37 ng/mL, closely matching the
population PK/PD estimates of 91.78% and 5.120 ng/mL encoded here.

# PKNCA validation

Both endpoints are run through PKNCA on the study's own sampling schedule. For
the PK endpoint the comparison is reported on **two scales**: the individual
prediction (`IPRED`, no residual error) and the simulated observation
(`Observed`, with the model's proportional residual error and the 0.5 ng/mL
assay lower limit applied). Table 2 of the paper reports NCA of *measured*
concentrations, so the `Observed` row is the like-for-like comparator; the
`IPRED` row shows how much of any discrepancy is attributable to residual error
rather than to the structural model. No parameter is tuned to either.

```{r pknca-setup}
# Restrict to the study's sampling times; drop the figure-only grid.
pk_nca_in <- pk |> filter(time %in% t_nca)
pd_nca_in <- pd |> filter(time %in% t_nca)

dose_df <- events |>
  filter(evid == 1) |>
  select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

intervals_pk <- data.frame(
  start     = c(0,     t_d14,        t_d14 + tau),
  end       = c(tau,   t_d14 + tau,  Inf),
  cmax      = c(TRUE,  TRUE,         FALSE),
  tmax      = c(TRUE,  TRUE,         FALSE),
  auclast   = c(TRUE,  TRUE,         FALSE),
  cav       = c(FALSE, TRUE,         FALSE),
  half.life = c(FALSE, FALSE,        TRUE)
)

run_nca <- function(df, value, intervals) {
  d <- df |>
    transmute(id, treatment, time, conc = .data[[value]]) |>
    filter(!is.na(conc)) |>
    arrange(id, treatment, time)
  res <- PKNCA::pk.nca(
    PKNCA::PKNCAdata(
      PKNCA::PKNCAconc(d, conc ~ time | treatment + id),
      dose_obj,
      intervals = intervals
    )
  )
  # PKNCA also returns dependency rows -- parameters it had to compute to
  # satisfy a request (e.g. tmax and lambda.z for half.life). Those carry the
  # interval they were computed on, so they must be filtered on
  # start/end AND parameter, not on parameter alone: otherwise the terminal
  # interval's tmax would be pooled with the requested dosing-interval tmax.
  requested <- intervals |>
    tidyr::pivot_longer(
      -c(start, end), names_to = "PPTESTCD", values_to = "asked"
    ) |>
    filter(asked) |>
    select(start, end, PPTESTCD)

  as.data.frame(res) |>
    filter(!is.na(PPORRES)) |>
    inner_join(requested, by = c("start", "end", "PPTESTCD")) |>
    mutate(
      group = paste0(
        treatment,
        if_else(start == 0, " - Day 1", " - Day 14")
      )
    ) |>
    select(group, id, PPTESTCD, PPORRES)
}
```

```{r pknca-pk}
# Aggregate to the group level using the SAME statistic the paper reports:
# Table 2 gives mean +/- SD for every parameter except Tmax, which is a median.
aggregate_like_paper <- function(x) {
  x |>
    group_by(group, PPTESTCD) |>
    summarise(
      PPORRES = if (unique(PPTESTCD) == "tmax") median(PPORRES) else mean(PPORRES),
      .groups = "drop"
    )
}

nca_pk <- bind_rows(
  run_nca(pk_nca_in, "ipred", intervals_pk) |>
    aggregate_like_paper() |> mutate(scale = "IPRED"),
  run_nca(pk_nca_in, "observed", intervals_pk) |>
    aggregate_like_paper() |> mutate(scale = "Observed")
)

published_pk <- tibble::tribble(
  ~group,                          ~cmax, ~tmax, ~auclast, ~cav,  ~half.life,
  "Cetagliptin 50 mg - Day 1",      80.5,  2.00,      717,   NA,          NA,
  "Cetagliptin 100 mg - Day 1",    219.0,  1.00,     1830,   NA,          NA,
  "Cetagliptin 50 mg - Day 14",    162.0,  1.00,     1530, 63.9,        41.9,
  "Cetagliptin 100 mg - Day 14",   300.0,  1.00,     3120,  130,        34.9
) |>
  tidyr::crossing(scale = c("IPRED", "Observed"))

cmp_pk <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_pk,
  reference = published_pk,
  by        = c("group", "scale"),
  units     = c(cmax = "ng/mL", auclast = "h*ng/mL", cav = "ng/mL",
                tmax = "h", half.life = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_pk,
  caption = paste(
    "Simulated vs. published cetagliptin NCA (Zhou 2024 Table 2), on both the",
    "IPRED and the residual-error (Observed) scale.",
    "* differs from reference by more than 20%."
  )
)
```

`r attr(cmp_pk, "footnote")`

## Plasma DPP-4 inhibition

The paper reports the DPP-4 inhibition endpoint with the same non-compartmental
descriptors (Rmax, TRmax, AUEC0-24h), so the same machinery applies with the
effect measure in place of concentration. Only the IPRED scale is compared here:
the published PD residual standard deviation of 11.62 percentage points is large
relative to the 80-90% inhibition plateau, so simulated observations routinely
exceed 100% inhibition, which is not a physically attainable measurement and
would make an "observed-scale" Rmax meaningless.

```{r pknca-pd}
intervals_pd <- data.frame(
  start   = c(0,   t_d14),
  end     = c(tau, t_d14 + tau),
  cmax    = TRUE,
  tmax    = TRUE,
  auclast = TRUE
)

nca_pd <- run_nca(pd_nca_in, "ipred", intervals_pd) |> aggregate_like_paper()

published_pd <- tibble::tribble(
  ~group,                         ~cmax, ~tmax, ~auclast,
  "Cetagliptin 50 mg - Day 1",    86.39,  2.00,     1820,
  "Cetagliptin 100 mg - Day 1",   88.78,  1.00,     2000,
  "Cetagliptin 50 mg - Day 14",   89.47,  2.00,     2010,
  "Cetagliptin 100 mg - Day 14",  89.99,  2.00,     2090
)

cmp_pd <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_pd,
  reference = published_pd,
  by        = "group",
  units     = c(cmax = "%", auclast = "h*%", tmax = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_pd,
  caption = paste(
    "Simulated vs. published DPP-4 inhibition NCA (Zhou 2024 Table 3).",
    "Rmax maps to Cmax, TRmax to Tmax, and AUEC0-24h to AUClast.",
    "* differs from reference by more than 20%."
  )
)
```

`r attr(cmp_pd, "footnote")`

The pharmacodynamic endpoint -- the quantity the PK/PD model was actually built
to describe -- reproduces the published Rmax and AUEC values at both dose levels
and on both study days to within a few percent. The pharmacokinetic comparison
is good on day 1 (AUC within about 10% at both doses) and reasonable at steady
state, with two structural features of the published parameter set visible:

* Day-1 Cmax is underpredicted on the IPRED scale but much closer on the
  observed scale, which is the correct comparator for Table 2. Roughly half of
  the apparent Cmax gap is simply residual error, not structural misfit.
* At steady state the 50 mg arm is low and the 100 mg arm high, i.e. the model
  predicts more dose-dependent accumulation than the NCA found. This is inherent
  in the published estimates rather than in this encoding: Michaelis-Menten
  elimination necessarily accumulates *more* at the higher dose, whereas the
  observed accumulation ratios run the other way (R_AUC 2.13 at 50 mg versus
  1.75 at 100 mg, Table 2). A single Vmax/Km pair fit to both dose groups can
  only split the difference. The paper's own goodness-of-fit assessment
  (Results section 3.6.1) notes "individual data deviations", and six records
  with |CWRES| > 5 were excluded before the PK/PD step.

Every simulated Tmax / TRmax falls on an adjacent sampling time to the published
median (1 h versus 2 h, or 2 h versus 1 h) and inside the published range
(1.00-4.00 h on day 1 and 0.500-5.00 h at steady state for Cmax; 0.500-3.00 h and
0.500-6.00 h for the PD TRmax). Because the schedule is coarse near the peak, a
single-interval shift registers as a 50-100% relative difference and is flagged
by the 20% tolerance; the flag reflects the resolution of the sampling grid, not
a disagreement about when the peak occurs.

No parameter has been tuned to improve any of these comparisons.

# Covariate effect

Total bilirubin is the only covariate retained in the final model, acting on the
peripheral volume of distribution.

```{r covariate-tbili}
tibble(TBILI = c(5, 10, 15, 20)) |>
  mutate(vp = 558.0 * (TBILI / tbili_ref)^0.3723) |>
  dplyr::rename(
    "Total bilirubin (umol/L)" = TBILI,
    "Peripheral volume V2 (L)" = vp
  ) |>
  knitr::kable(
    digits = 1,
    caption = paste(
      "Effect of total bilirubin on the peripheral volume of distribution,",
      "V2 = 558.0 * (TBILI / 10)^0.3723. Consistent with the paper's",
      "statement that V2 increases with increasing TBIL."
    )
  )
```

# Assumptions and deviations

- **Total bilirubin centring value (10 umol/L) is an assumption.** Zhou 2024
  reports the covariate coefficient (`dV2dTBIL` = 0.3723, Supplementary
  Table 1) and the direction of the effect (Discussion page 10: "V2 increases
  with the increase of TBIL") but publishes **neither the covariate equation nor
  the population median total bilirubin** -- TBIL does not appear in the Table 1
  demographics. Two things had to be settled. (i) *Functional form*: the
  parameter name follows Phoenix NLME's `dPARAMdCOVARIATE` convention, and
  Phoenix's default continuous-covariate transformation is the median-centred
  log-ratio `V2 = tvV2 * exp(dV2dTBIL * log(TBIL / median(TBIL)))`, which is
  algebraically the power form `V2 = tvV2 * (TBIL / ref)^0.3723` used here. The
  alternative linear-exponential reading is untenable: over the normal
  bilirubin range it would inflate V2 by more than 40-fold. (ii) *Centring
  value*: 10 umol/L was chosen as a rounded mid-normal adult total bilirubin, so
  that the Supplementary Table 1 typical value `tvV2` = 558.0 L is the typical
  peripheral volume of the study population. Results section 3.2 states that no
  subject had clinically significant abnormal liver function, so observed TBIL
  values lay within the normal 5-21 umol/L range and the choice of centre only
  shifts V2 modestly (a 2-fold error in the reference changes V2 by 29%).
- **The PD residual error is encoded as additive only.** Supplementary Table 1
  labels the PD residual model "MixRatio residual variability PD(sigma)" and
  reports a single value, `Stdev0` = 11.62. Phoenix's mixed-ratio observation
  model carries both this additive standard deviation and a separate `MixRatio`
  coefficient scaling a proportional component; the `MixRatio` coefficient is
  not reported anywhere in the paper or the supplement. Only the additive term
  is therefore encoded, and no proportional coefficient has been invented. The
  additive reading is the only dimensionally sensible one: 11.62 is on the
  DPP-4 inhibition percentage-point scale, and as a proportional fraction it
  would imply a 1162% CV. Because 11.62 percentage points is large next to the
  80-90% inhibition plateau, the residual-error scale is not used for the PD
  NCA comparison (see that section).
- **The printed sigmoid Emax equation drops an exponent.** Methods section 2.7.2
  typesets the PD model as `Emax * C^g / (EC50 + C^g)`, which is dimensionally
  inconsistent and contradicts the paper's own definition of EC50 as the "plasma
  concentration of cetagliptin that achieves 50% of the maximum drug effect".
  The canonical Phoenix Sigmoid-Emax form `Emax * C^g / (EC50^g + C^g)` named in
  Results section 3.6.2 is used instead. With the estimated Hill coefficient of
  1.008 the two forms differ by under 0.1%, so nothing material turns on the
  choice.
- **`DUR80%` is compared on a censored basis.** The published day-14 `DUR80%`
  values (21.9 h at 50 mg, 32.3 h at 100 mg) exceed the 24 h dosing interval,
  so they must have been computed over the whole day-14 profile including the
  120 h washout rather than over one interval. The simulated value is reported
  over the dosing interval only and is therefore bounded above by 24 h; the two
  are not directly comparable at 100 mg and the table says so.
- **Eight screened covariates carry no encodable effect.** Discussion page 10
  lists the full stepwise candidate set: sex, body weight, ALT
  ("glutamic-pyruvic transaminase"), total bilirubin, triglycerides, LDL
  cholesterol, glucose, urea and creatinine. Only total bilirubin was retained,
  and the paper reports no point estimate, standard error or confidence interval
  for any of the other eight. They are therefore recorded in the model's
  `covariatesDataExcluded` metadata -- which documents the screen without
  declaring an unused covariate -- rather than in `covariateData`. Units for the
  screened-only entries are given per the SI conventions used by Chinese
  clinical-chemistry laboratories and are flagged as assumed in each entry's
  notes, because the paper never states them (no effect was retained, so it had
  no reason to).
- **Subject-count discrepancy in the source.** Methods section 2.7.1 states that
  the 560 cetagliptin concentrations used for the population PK analysis came
  from "32 patients with T2DM". Only 20 of the 32 enrolled subjects received
  cetagliptin (Table 1 and Table 2; the remainder received sitagliptin or
  placebo), so the 32 cannot be the number of subjects contributing cetagliptin
  concentrations. The model's `population$n_subjects` records the 20 cetagliptin
  recipients and the paper's statement is preserved in `population$notes`.
- **Elimination is Michaelis-Menten, which the main text does not state.** The
  article describes only a "two-compartment model"; that elimination is
  saturable is inferable solely from Supplementary Table 1 reporting `tvKm` and
  `tvVmax` with no clearance parameter. This is encoded as published, and the
  steady-state mass-balance check above confirms the parameterisation is
  internally consistent.
- **Dose units.** Supplementary Table 1 parameterises the model in ug (Vmax in
  ug/h) and ug/L (Km, EC50), while doses are naturally entered in mg. The model
  applies `f(depot) <- 1000` purely as a mg-to-ug unit conversion; it is not a
  bioavailability estimate. All disposition parameters remain apparent
  (oral, /F) because the paper never estimated bioavailability.
- **Virtual total bilirubin distribution.** The paper reports no TBIL summary
  statistics, so the virtual cohort draws TBILI from a log-normal distribution
  with median 10 umol/L and log-scale SD 0.30, truncated to 3-21 umol/L. This
  affects only the width of the simulated prediction intervals, not the typical
  profile.
- **Sitagliptin is not extracted.** Sitagliptin 100 mg was the active comparator
  and the paper reports its NCA parameters (Table 2) and a non-compartmental
  concentration-effect Emax fit (Emax 91.68%, EC50 6.73 ng/mL, Results section
  3.4.1). No population PK model was developed for sitagliptin, so no
  simulatable model exists to extract.
- **Other endpoints are descriptive only.** Active GLP-1 (Table 4), the OGTT
  glucose / insulin / C-peptide / glucagon responses (Table 5), FPG, 2 h PPG,
  HbA1c and glycated albumin (Table 6) were analysed non-compartmentally and by
  ANOVA. The paper builds no structural model for any of them, so DPP-4
  inhibition is the only PD endpoint with a model to extract.
- **Half-life comparison.** The published t1/2 (34.9-41.9 h) is an NCA estimate
  from the day-14 washout. Under saturable elimination the apparent terminal
  half-life is concentration-dependent, so the simulated value is expected to
  differ somewhat from the observed NCA value even when the model is encoded
  correctly.
