---
title: "Brepocitinib (Maleki 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Brepocitinib (Maleki 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("Maleki_2024_brepocitinib"))$reference`
* Description: `r rxode2::rxode(readModelDb("Maleki_2024_brepocitinib"))$description`
* Article: <https://doi.org/10.1002/psp4.13100>
* Supplement (Tables S1-S3, Figures S1-S6): <https://doi.org/10.1002/psp4.13100> (Supporting Information)

Brepocitinib (PF-06700841) is an oral dual TYK2/JAK1 inhibitor. Maleki 2024
updates an earlier brepocitinib population PK model by adding the ulcerative
colitis, psoriatic arthritis, hidradenitis suppurativa, and vitiligo cohorts
that the prior model did not cover.

# Population

The analysis pooled nine clinical studies -- three phase I studies in healthy
participants (NCT02310750 single- and multiple-ascending dose plus a plaque
psoriasis arm, NCT03236493 in healthy Japanese participants, NCT03656952 the
thorough-QT crossover) and six phase II studies in plaque psoriasis
(NCT02969018), alopecia areata (NCT02974868), psoriatic arthritis
(NCT03963401), ulcerative colitis (NCT02958865), vitiligo (NCT03715829), and
hidradenitis suppurativa (NCT04092452).

Of 11,078 available plasma concentrations from 827 individuals, the analysis
retained 9532 records from 775 individuals (8552 of them plasma concentrations,
per the Abstract). Baseline demographics are Maleki 2024 Table 1: median age 43
years (range 18-75), median baseline body weight 81 kg (42-204), median baseline
albumin 4.5 g/dL (3.2-5.6), median baseline creatinine clearance 119 mL/min
(41.4-357), and median baseline aspartate aminotransferase 19 U/L (8-133). The
population was 41% female and predominantly White (85%, versus 7% Asian and 8%
Other). Patients contributed 88% of the individuals, split across psoriasis
(26%), psoriatic arthritis (24%), ulcerative colitis (19%), alopecia areata
(8%), vitiligo (7%), and hidradenitis suppurativa (5%); healthy participants
were the remaining 12%. Daily doses spanned under 10 mg to over 100 mg, with
60 mg (44%) and 30 mg (29%) the most common, and 92% of individuals received the
tablet formulation. The lower limit of quantification was 0.2 ng/mL; fewer than
10% of samples fell below it, so no BLQ imputation was applied.

Estimation used the stochastic approximation expectation maximization algorithm
in MonolixSuite 2021R1 -- not NONMEM -- which matters for reading the parameter
table: the Box-Cox random-effect transformation and the combined residual error
below are written in Monolix conventions.

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

```{r population}
str(readModelDb("Maleki_2024_brepocitinib")()$population)
```

# Source trace

Every `ini()` entry in
`inst/modeldb/specificDrugs/Maleki_2024_brepocitinib.R` carries an in-file
comment naming its source location. They are collected here for review.

Maleki 2024 Table 3 is internally inconsistent: its "Estimate" column disagrees
with the covariate equations printed in the same table's own footnotes c, d, e,
and f. The **Source used** column below records which one was taken and why; the
full reconciliation is in *Assumptions and deviations*.

| Equation / parameter | Value | Source used |
|---|---|---|
| Structural model: two compartments, first-order absorption, linear elimination, lag time on the tablet formulation | n/a | Results, "PK analysis"; Discussion |
| `lka` (Ka) | 2.66 1/h | Table 3 Estimate column (no footnote value adopted) |
| `lcl` (CL/F) | 17.49 L/h | Table 3 **footnote e** (Estimate column says 18; Abstract and Discussion say 17.5) |
| `lvc` (Vc/F) | 88.52 L | Table 3 **footnote f** (Estimate column says 91.1; Abstract says 88.5) |
| `lq` (Q/F) | 0.55 L/h | Table 3 Estimate column |
| `lvp` (Vp/F) | 232.5 L | Table 3 Estimate column; bootstrap 95% CI (12, 455) |
| `ltlag` (Tlag, tablet only) | 0.26 h | Table 3 Estimate column; footnote b: Tlag = 0 for suspension |
| `lfdepot` (Frel at <= 100 mg/day) | 1 (fixed) | Table 3 row "Frel" = 1 with no bootstrap entry |
| `e_wt_cl_q` | 0.75 (fixed) | Table 3 row "Effect of BWT on CL/F and Q/F", marked Fixed |
| `e_wt_vc_vp` | 1 (fixed) | Table 3 row "Effect of BWT on Vc/F and Vp/F", marked Fixed |
| `e_ast_cl` | -0.17 | Table 3 **footnote e**: `(BAST / 22)^-0.17` (Estimate column says -0.22) |
| `e_race_asian_cl` | -0.11 | Table 3 **footnote e**: `exp(-0.11 * RACE)` (Estimate column says -0.17) |
| `e_sexf_vc` | 0.12 | Table 3 **footnote f**: `exp(0.12 * SEX)`, SEX = 1 male (Estimate column says 0.10) |
| `e_dis_uc_ka` | -1.16 1/h | Table 3 **Estimate column** (1.16), sign from the bootstrap 95% CI (-1.58, -0.68); footnote d's -0.61 lies outside that CI |
| `e_dose_high_fdepot` | 0.3 | Table 3 **footnote c**: `F = 1 * (1 + 0.3 * D100)` (Estimate column says 0.26) |
| Box-Cox IIV: `P_i = theta_P * exp(eta_T)`, `eta_T = ((exp(eta))^theta_s - 1) / theta_s` | n/a | Methods, "Stochastic model development" (Petersson 2009 form) |
| `etalcl` SD (healthy participants) | 0.42 | Table 3 row "omega2 CL/F (SD)" |
| `etalvc` SD (healthy participants) | 0.19 | Table 3 row "omega2 Vc/F (SD)" |
| correlation(`etalcl`, `etalvc`) | 0.71 | Table 3 row "rho CL/F, Vc/F" |
| `boxcox_cl` | 0.26 | Table 3 row "Box-Cox transformation shape parameter on omega2 CL/F" |
| `boxcox_vc` | 5.86 | Table 3 row "Box-Cox transformation shape parameter on omega2 Vc/F" |
| `e_dis_healthy_etalcl` | 0.44 | Table 3 row "Effect of patients on omega2 CL/F" |
| `e_dis_healthy_etalvc` | 0.29 | Table 3 row "Effect of patients on omega2 Vc/F" |
| Residual error: `DV = IPRED + sqrt(sigma_add^2 + (sigma_pro * IPRED)^2) * eps` | n/a | Methods, "Stochastic model development" (= nlmixr2 `combined2`) |
| `addSdHv` / `propSdHv` | 0.89 ng/mL / 0.1 | Table 3 rows "Additive / Proportional residual error of HV (SD)" |
| `addSdPat` / `propSdPat` | 0.45 ng/mL / 0.05 | Table 3 rows "Additive / Proportional residual error of patients (SD)" |
| Covariate reference values 70 kg and 22 U/L | n/a | Discussion: "the typical values of 70 kg and 22 U/L were used for BWT and BAST" |
| Reference individual: 70 kg, White, non-UC, non-Asian, female, BAST 22 U/L | n/a | Methods, "Clinical trial simulations" |

```{r model-print}
mod <- readModelDb("Maleki_2024_brepocitinib")
rxode2::rxode(mod)
```

# Virtual cohort

Original observed data are not publicly available. The simulations below use
virtual cohorts whose covariates match either the paper's reference individual
or the specific covariate scenario under test. Population cohorts use 200
participants per arm.

```{r cohort-helpers}
set.seed(20240415)

REF_COVS <- list(
  WT = 70, AST = 22, RACE_ASIAN = 0, SEXF = 1,
  DIS_UC = 0, DIS_HEALTHY = 0, FORM_TABLET = 1, DOSE_HIGH = 0
)

# Build one arm: q.d. dosing at `dose_times`, observations at `obs_times`.
# `covs` overrides individual entries of REF_COVS.
build_arm <- function(ids, dose, dose_times, obs_times, arm, covs = list()) {
  cv <- utils::modifyList(REF_COVS, covs)
  dosing <- tidyr::expand_grid(id = ids, time = dose_times) |>
    dplyr::mutate(amt = dose, evid = 1L, cmt = "depot")
  obs <- tidyr::expand_grid(id = ids, time = obs_times) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = "central")
  dplyr::bind_rows(dosing, obs) |>
    dplyr::mutate(
      WT = cv$WT, AST = cv$AST, RACE_ASIAN = cv$RACE_ASIAN, SEXF = cv$SEXF,
      DIS_UC = cv$DIS_UC, DIS_HEALTHY = cv$DIS_HEALTHY,
      FORM_TABLET = cv$FORM_TABLET, DOSE_HIGH = cv$DOSE_HIGH,
      arm = arm, dose_mg = dose
    ) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

mod_typ <- rxode2::zeroRe(mod)
```

# Replicate published figures

## Figure 1 -- concentration-time profiles after single and multiple doses

Maleki 2024 Figure 1 plots observed dose-normalized concentrations after single
or multiple doses across the pooled studies. The typical-value replication below
covers the phase I single-dose range (1 to 200 mg) so the dose-dependent
bioavailability step above 100 mg is visible on the dose-normalized scale.

```{r figure-1}
single_doses <- c(1, 10, 30, 100, 200)
sd_events <- dplyr::bind_rows(lapply(seq_along(single_doses), function(i) {
  d <- single_doses[i]
  build_arm(
    ids = i, dose = d, dose_times = 0,
    obs_times = seq(0, 96, by = 0.25),
    arm = paste0(d, " mg"),
    covs = list(DOSE_HIGH = as.integer(d > 100))
  )
}))
stopifnot(!anyDuplicated(unique(sd_events[, c("id", "time", "evid")])))

sd_sim <- rxode2::rxSolve(
  mod_typ, events = sd_events, omega = NA,
  keep = c("arm", "dose_mg"), returnType = "data.frame"
)
stopifnot(length(unique(sd_sim$id)) == length(single_doses))

sd_sim |>
  # Blank the pre-dose zero rather than filtering, so no row is dropped from
  # anything downstream; log10 simply omits the NA.
  dplyr::mutate(Cc = ifelse(Cc > 0, Cc, NA_real_),
                arm = factor(arm, levels = paste0(single_doses, " mg"))) |>
  ggplot(aes(time, Cc / dose_mg, colour = arm)) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)", y = "Dose-normalized Cc (ng/mL per mg)",
    colour = "Single dose",
    title = "Typical-value single-dose profiles",
    caption = paste(
      "Replicates the single-dose panel of Figure 1 of Maleki 2024.",
      "The 200 mg curve sits 30% above the others on the dose-normalized",
      "scale because DOSE_HIGH = 1 raises relative bioavailability."
    )
  ) +
  theme_bw()
```

## Figure 3 -- prediction-corrected visual predictive check

Maleki 2024 Figure 3 is a pcVPC of the index dataset. The replication below is a
straight VPC of the typical phase II regimen (30 mg q.d.) in patients, showing
the median and the 5th/95th percentiles of the simulated concentrations on day 1
and over the 14th dosing interval.

```{r figure-3}
n_pop <- 200
vpc_obs <- c(seq(0, 24, by = 1), seq(312, 336, by = 1))
vpc_events <- build_arm(
  ids = seq_len(n_pop), dose = 30, dose_times = seq(0, 13 * 24, by = 24),
  obs_times = vpc_obs, arm = "30 mg q.d. (patients)"
)

vpc_sim <- rxode2::rxSolve(
  mod, events = vpc_events, keep = c("arm"), returnType = "data.frame"
)
stopifnot(length(unique(vpc_sim$id)) == n_pop)

vpc_sim |>
  dplyr::mutate(day = ifelse(time <= 24, "Day 1 (first dose)", "Day 14 (last dose)"),
                tad = ifelse(time <= 24, time, time - 312),
                Cc = ifelse(Cc > 0, Cc, NA_real_)) |>
  dplyr::group_by(day, tad) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE), Q50 = median(Cc, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(tad, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(linewidth = 0.7) +
  facet_wrap(~day) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)", y = "Cc (ng/mL)",
    title = "VPC of 30 mg q.d. in patients (n = 200)",
    caption = paste(
      "Companion to Figure 3 of Maleki 2024 (pcVPC of the index dataset).",
      "Ribbon spans the 5th to 95th simulated percentiles."
    )
  ) +
  theme_bw()
```

## Healthy participants versus patients -- the two variability layers

`DIS_HEALTHY` changes no typical structural parameter. It scales the Box-Cox
random effects on CL/F and Vc/F (Table 3 rows "Effect of patients on omega2")
and selects which pair of residual-error magnitudes applies. Simulating the same
regimen in each cohort isolates the first of those two effects.

```{r health-status-iiv}
hv_events <- build_arm(
  ids = 1000L + seq_len(n_pop), dose = 30,
  dose_times = seq(0, 13 * 24, by = 24), obs_times = vpc_obs,
  arm = "30 mg q.d. (healthy)", covs = list(DIS_HEALTHY = 1)
)
hv_sim <- rxode2::rxSolve(
  mod, events = hv_events, keep = c("arm"), returnType = "data.frame"
)
stopifnot(length(unique(hv_sim$id)) == n_pop)

iiv_summary <- dplyr::bind_rows(vpc_sim, hv_sim) |>
  dplyr::group_by(arm, id) |>
  dplyr::summarise(cl = dplyr::first(cl), vc = dplyr::first(vc), .groups = "drop") |>
  dplyr::group_by(arm) |>
  dplyr::summarise(
    `CL/F median (L/h)`  = round(median(cl), 2),
    `CL/F 5th-95th`      = paste0(round(quantile(cl, 0.05), 1), " - ",
                                  round(quantile(cl, 0.95), 1)),
    `Vc/F median (L)`    = round(median(vc), 1),
    `Vc/F 5th-95th`      = paste0(round(quantile(vc, 0.05), 1), " - ",
                                  round(quantile(vc, 0.95), 1)),
    .groups = "drop"
  ) |>
  dplyr::rename(Cohort = arm)

knitr::kable(
  iiv_summary,
  caption = paste(
    "Simulated individual CL/F and Vc/F by health status (n = 200 per arm).",
    "Patients carry a 1 + 0.44 uplift on the CL/F random effect and a",
    "1 + 0.29 uplift on the Vc/F random effect relative to healthy",
    "participants (Maleki 2024 Table 3)."
  )
)
```

# PKNCA validation

The paper reports no absolute noncompartmental Cmax / Tmax / AUC / half-life
values, so the validation target is Maleki 2024 Figure 4: geometric mean ratios
of dose-normalized steady-state AUCtau and Cmax for each covariate scenario
relative to the reference population (70 kg White non-UC non-Asian female with
BAST 22 U/L). Following the paper's Methods, each scenario is dosed q.d. for two
weeks and the exposure metrics are taken over the last dosing interval
(312 to 336 h). NCA is computed with PKNCA over that interval.

```{r nca-scenarios}
scenarios <- tibble::tribble(
  ~arm,                       ~dose, ~covs,
  "Reference",                   30, list(),
  "BWT 53 kg (5th pctile)",      30, list(WT = 53),
  "BWT 121 kg (95th pctile)",    30, list(WT = 121),
  "BAST 8 U/L (1st pctile)",     30, list(AST = 8),
  "BAST 65 U/L (99th pctile)",   30, list(AST = 65),
  "Asian",                       30, list(RACE_ASIAN = 1),
  "Male",                        30, list(SEXF = 0),
  "Ulcerative colitis",          30, list(DIS_UC = 1),
  "Daily dose > 100 mg",        200, list(DOSE_HIGH = 1)
)

ss_dose_times <- seq(0, 13 * 24, by = 24)
ss_obs_times <- seq(312, 336, by = 0.1)

ss_events <- dplyr::bind_rows(lapply(seq_len(nrow(scenarios)), function(i) {
  build_arm(
    ids = i, dose = scenarios$dose[i], dose_times = ss_dose_times,
    obs_times = ss_obs_times, arm = scenarios$arm[i],
    covs = scenarios$covs[[i]]
  )
}))
stopifnot(!anyDuplicated(unique(ss_events[, c("id", "time", "evid")])))

ss_sim <- rxode2::rxSolve(
  mod_typ, events = ss_events, omega = NA,
  keep = c("arm", "dose_mg"), returnType = "data.frame"
)
stopifnot(length(unique(ss_sim$id)) == nrow(scenarios))
```

```{r pknca}
# Concentrations. Filter on !is.na(Cc) only -- a `time > 0` or `Cc > 0` filter
# would drop the interval-start record PKNCA anchors the AUC on.
ss_conc <- ss_sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, arm)

ss_dose <- ss_events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, arm)

conc_obj <- PKNCA::PKNCAconc(ss_conc, Cc ~ time | arm + id)
dose_obj <- PKNCA::PKNCAdose(ss_dose, amt ~ time | arm + id)

intervals <- data.frame(
  start   = 312, end = 336,
  cmax    = TRUE, tmax = TRUE, cmin = TRUE,
  auclast = TRUE, cav = TRUE
)

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

nca_wide <- as.data.frame(nca_res) |>
  dplyr::select(arm, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::left_join(dplyr::distinct(ss_events, arm, dose_mg), by = "arm") |>
  dplyr::mutate(
    cmax_dn    = cmax / dose_mg,
    auclast_dn = auclast / dose_mg,
    # PKNCA reports tmax as time after the interval's dose, not absolute time.
    tmax_pd    = tmax,
    arm        = factor(arm, levels = scenarios$arm)
  ) |>
  dplyr::arrange(arm)

knitr::kable(
  nca_wide |>
    dplyr::transmute(
      Scenario                          = as.character(arm),
      `Dose (mg q.d.)`                  = dose_mg,
      `Cmax,ss (ng/mL)`                 = round(cmax, 1),
      `Tmax,ss (h post-dose)`           = round(tmax_pd, 2),
      `Cmin,ss (ng/mL)`                 = round(cmin, 2),
      `AUCtau,ss (ng*h/mL)`             = round(auclast, 0),
      `Dose-norm. AUCtau,ss (ng*h/mL/mg)` = round(auclast_dn, 2)
    ),
  caption = paste(
    "PKNCA results over the 14th dosing interval (312-336 h) for the",
    "typical-value profile of each Maleki 2024 covariate scenario."
  )
)
```

## Comparison against published exposure ratios

```{r nca-ratios}
ref_row <- nca_wide |> dplyr::filter(arm == "Reference")

sim_ratios <- nca_wide |>
  dplyr::filter(arm != "Reference") |>
  dplyr::transmute(
    scenario = as.character(arm),
    cmax     = cmax_dn / ref_row$cmax_dn,
    auclast  = auclast_dn / ref_row$auclast_dn
  )

# The paper gives numeric geometric mean ratios only for the two body-weight
# scenarios (Discussion) and the dose-dependent bioavailability scenario
# (Discussion). The remaining scenarios are reported qualitatively or as CI
# bounds and are discussed in the narrative below.
published_ratios <- tibble::tribble(
  ~scenario,                    ~cmax, ~auclast,
  "BWT 53 kg (5th pctile)",       1.40,     1.30,
  "BWT 121 kg (95th pctile)",     0.65,     0.70,
  "Daily dose > 100 mg",          1.30,     1.30
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = sim_ratios,
  reference = published_ratios,
  by        = "scenario",
  units     = c(cmax = "ratio", auclast = "ratio"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste(
    "Simulated versus published geometric mean ratios of dose-normalized",
    "steady-state Cmax and AUCtau relative to the reference population",
    "(Maleki 2024 Discussion and Figure 4).",
    "* marks a difference greater than 20%."
  )
)
```

All three fully-quantified scenarios reproduce within 20%. The simulated
body-weight ratios (AUCtau 1.23 and 0.66; Cmax 1.30 and 0.59) sit just inside
the paper's one-decimal values of 1.3 / 0.7 and 1.4 / 0.65; the residual gap is
expected because the paper's numbers are geometric means over 1000 trials of 220
individuals carrying the full random-effect distribution, whereas these are
deterministic typical-value ratios. The dose-dependent bioavailability scenario
reproduces exactly (1.300 against 1.3), because at steady state
AUCtau is exactly proportional to relative bioavailability.

The remaining scenarios can only be checked against the paper's prose:

```{r nca-ratios-qualitative}
qual <- sim_ratios |>
  dplyr::filter(!scenario %in% published_ratios$scenario) |>
  dplyr::mutate(
    `Published claim (Maleki 2024)` = dplyr::case_when(
      scenario == "BAST 8 U/L (1st pctile)"   ~ "AUCtau CI reaches 0.78; Cmax not significantly different",
      scenario == "BAST 65 U/L (99th pctile)" ~ "AUCtau CI reaches 1.3; Cmax not significantly different",
      scenario == "Asian"                     ~ "~10% slower CL/F; AUCtau not significantly different",
      scenario == "Male"                      ~ "13% higher Vc/F; Cmax and AUCtau within 20%",
      scenario == "Ulcerative colitis"        ~ "AUCtau comparable; Cmax 'marginally lower' (Discussion) / 'average 20% lower' (Conclusions)"
    )
  ) |>
  dplyr::transmute(
    Scenario                = scenario,
    `Simulated Cmax ratio`  = round(cmax, 3),
    `Simulated AUCtau ratio` = round(auclast, 3),
    `Published claim (Maleki 2024)`
  )

knitr::kable(
  qual,
  caption = paste(
    "Scenarios for which Maleki 2024 reports a qualitative claim or a",
    "confidence-interval bound rather than a point ratio."
  )
)
```

Each simulated ratio is consistent with the corresponding claim, with one
exception discussed below: the ulcerative colitis Cmax ratio of 0.905 matches the
Discussion's "marginally lower" but not the Conclusions' "average 20% lower".
The two BAST point ratios (0.844 and 1.199) sit just inside the CI bounds of
0.78 and 1.3 that the paper reports for those scenarios, exactly as expected for
point estimates whose confidence intervals extend to those bounds.

## Figure 4 -- covariate forest plot

```{r figure-4, fig.height = 4}
sim_ratios |>
  tidyr::pivot_longer(c(cmax, auclast), names_to = "metric", values_to = "ratio") |>
  dplyr::mutate(
    metric = factor(metric, levels = c("cmax", "auclast"),
                    labels = c("Cmax,ss", "AUCtau,ss")),
    scenario = factor(scenario, levels = rev(sim_ratios$scenario))
  ) |>
  ggplot(aes(ratio, scenario, shape = metric, colour = metric)) +
  annotate("rect", xmin = 0.8, xmax = 1.25, ymin = -Inf, ymax = Inf,
           fill = "darkgreen", alpha = 0.12) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_point(size = 3, position = position_dodge(width = 0.4)) +
  scale_x_continuous(limits = c(0.5, 1.5)) +
  labs(
    x = "Ratio to reference population", y = NULL,
    shape = NULL, colour = NULL,
    title = "Dose-normalized steady-state exposure ratios",
    caption = paste(
      "Replicates Figure 4 of Maleki 2024 (typical-value ratios rather than",
      "geometric means over simulated trials). Shaded band is the 0.8 to 1.25",
      "clinical-indifference range."
    )
  ) +
  theme_bw()
```

# Assumptions and deviations

## Maleki 2024 Table 3 contradicts its own footnote equations

This is the load-bearing decision in the extraction. Table 3 reports an
"Estimate" column, a "Bootstrap median" column, and a bootstrap 95% CI. Four of
its footnotes (c, d, e, f) print the corresponding covariate equations with
*different* numbers. Neither set can be derived from the other, and the raw PDF
confirms the discrepancy is in the published table, not in text extraction.

| Quantity | Table 3 Estimate | Bootstrap median | Footnote equation | Paper's prose |
|---|---|---|---|---|
| CL/F | 18 L/h | 17.5 | **17.49** (fn e) | 17.5 (Abstract, Discussion) |
| Vc/F | 91.1 L | 91.3 | **88.52** (fn f) | 88.5 (Abstract), 88.2 (Discussion) |
| BAST exponent on CL/F | -0.22 | -0.22 | **-0.17** (fn e) | not stated numerically |
| Asian effect on CL/F | -0.17 | -0.19 | **-0.11** (fn e) | "~10% lower CL/F"; exp(-0.11) = 0.896 |
| Male effect on Vc/F | 0.10 | 0.09 | **0.12** (fn f) | "13% higher Vc/F"; exp(0.12) = 1.128 |
| Frel effect above 100 mg | 0.26 | 0.29 | **0.30** (fn c) | "30% higher bioavailability", "1.3-fold" |
| UC effect on Ka | **1.16** | 1.14, CI (-1.58, -0.68) | -0.61 (fn d) | "46% slower absorption rate" |

The footnote equations were adopted for CL/F, Vc/F, BAST, race, sex, and Frel,
because on every quantity the paper states numerically the footnote value -- not
the Estimate column -- is the one the Abstract and Discussion reproduce. That is
five independent corroborations (17.5 L/h, 88.5 L, ~10%, 13%, 30%/1.3-fold).
Simulating the alternative shows the Estimate column would contradict the prose:
it gives a 15.6% Asian clearance reduction against "~10%", a 10.5% male Vc/F
increase against "13%", and a 1.26-fold bioavailability ratio against "1.3-fold".

**Footnote d is the exception and was rejected.** Its coefficient of -0.61 falls
outside the bootstrap 95% CI of (-1.58, -0.68) printed on the same table row,
which is internally impossible, and it implies only 23% slower absorption against
the paper's stated 46%. The Table 3 Estimate of 1.16 was used instead, with the
negative sign taken from the CI (the Estimate column has lost the sign: every
value in that CI is negative). Ka = 2.66 - 1.16 = 1.50 1/h gives 44% slower
absorption, matching the stated 46% to within rounding.

Parameters with no footnote equation (Ka, Q/F, Vp/F, Tlag, the Box-Cox shapes,
the omegas, rho, and the residual-error SDs) take the Table 3 Estimate column,
which is their only source. Ka's typical value therefore comes from the Estimate
column (2.66) rather than footnote d's 2.62, keeping it consistent with the UC
coefficient taken from the same column.

## Residual error is encoded in linear concentration space

The Methods section is headed "log-transform combined error model", but the
equation it prints is

```
DV_ij = IPRED_ij + sqrt(sigma_add^2 + (sigma_pro * IPRED_ij)^2) * eps_ij
```

which is a combined error in the *linear* concentration space -- exactly
nlmixr2's default `combined2` form for `add(addSd) + prop(propSd)`. The printed
equation was followed. Reading the additive SDs on the log scale instead would
make them implausible (an additive SD of 0.89 in log units is over 100% CV,
which is incompatible with the paper's diagnostic plots); read in ng/mL they are
entirely sensible next to the 0.2 ng/mL LLOQ, giving roughly 10% total residual
CV for healthy participants at a typical 30 mg Cmax of about 275 ng/mL. Because
the model's amounts are in mg and volumes in L, `model()` multiplies
`central / vc` by 1000 to put `Cc` in ng/mL, the analysis-dataset unit in which
those SDs are reported.

The direction of the residual-error split is reported as-is and contradicts the
paper's stated rationale. The Results text says "Historical patient data elicits
a higher degree of random variability compared to phase I HV studies", but
Table 3 gives healthy participants the *larger* residual error on both
components (0.89 / 0.10 versus 0.45 / 0.05 for patients). The tabulated values
were used. The direction is defensible on sampling grounds -- the phase I studies
sampled out to 96 h, deep into the near-LLOQ range where the additive term
dominates -- and the higher patient variability the text describes does appear,
in the interindividual rather than the residual layer.

## Box-Cox random effects

The Box-Cox transformation of the random effect is applied explicitly inside
`model()` because rxode2 has no native Box-Cox eta distribution. The
implementation reproduces the paper's printed equations exactly, including the
health-status uplift on the eta magnitude (scaling a mean-zero normal deviate
multiplies its standard deviation by the same factor, so
`etalcl * (1 + 0.44)` is precisely the patient omega). This means the etas are
not mu-referenced; the model is intended for simulation, and re-estimation with
it would need FOCEi rather than SAEM.

A consequence worth flagging to users: the shape parameter on the Vc/F random
effect is 5.86, which under the Petersson form produces a very heavy right tail.
The 5th-to-95th percentile spread of individual Vc/F is a reasonable 77 to
about 217 L, but roughly 1% of simulated participants exceed 700 L and rare draws
run into the thousands. Summarise population simulations from this model with
medians and percentiles, not arithmetic means. The IIV magnitudes and the shape
parameters are reproduced exactly as published; no smoothing was applied.

## Interindividual variability that the paper did not carry forward

The paper estimated interindividual variability on Ka and reports that it
significantly reduced AIC, but "the random effect on Ka was not carried to the
final model due to a high shrinkage". The final model in Table 3 therefore has
random effects only on CL/F and Vc/F, and that is what is encoded. No random
effects were invented for Q/F, Vp/F, Tlag, or Frel.

## Ulcerative colitis Cmax

The paper makes two different quantitative claims about the same result. The
Discussion says the UC Cmax was "marginally lower" after 30 mg daily for two
weeks; the Conclusions say the lower absorption rate "resulted in an average 20%
lower Cmax". The packaged model gives a 9.5% reduction on a typical-value basis,
consistent with "marginally lower" and with Figure 4 placing UC inside the
0.8-1.25 band. Sweeping the coefficient shows that no value inside the reported
bootstrap CI reaches 20%: the CI's lower bound of -1.58 gives 16%, and reaching
20% would require about -1.8, outside the CI. The published coefficient was kept.

## Covariate reference values and scenario definitions

* The allometric and AST reference values are 70 kg and 22 U/L, which the
  Discussion states were chosen deliberately in preference to the observed
  medians of 81 kg and 19 U/L.
* The paper's Methods say the dose-dependent bioavailability scenario was
  simulated at 100 mg, but footnote c defines the indicator as dose strictly
  greater than 100 mg, so a 100 mg subject would have `DOSE_HIGH = 0`. The
  scenario above is therefore simulated at 200 mg with `DOSE_HIGH = 1`, matching
  the Discussion's "individuals who received more than 100 mg q.d.". Because the
  ratio is taken on dose-normalized metrics, the choice of high dose does not
  affect the result.
* The body-weight scenarios use the 5th and 95th percentiles that the Discussion
  names explicitly (53 and 121 kg). Maleki 2024 Figure 4 uses the 5th and 95th
  percentiles for BAST too, but does not print them; the 1st and 99th percentiles
  it does print (8 and 65 U/L) are used instead, which is why the BAST rows are
  compared against CI bounds rather than point ratios.
* Sex is stored under the canonical `SEXF` (1 = female) while the paper's source
  column is a male indicator, so the effect is applied as
  `exp(0.12 * (1 - SEXF))`. This preserves both the paper's positive coefficient
  and its female reference category, under which the typical Vc/F of 88.52 L is a
  female value. The paper notes the reference sex "was selected alphabetically",
  not by group size -- males were 59% of the cohort.

## Validation scope

* The paper reports no absolute noncompartmental Cmax, Tmax, AUC, or half-life
  values for brepocitinib, so no absolute NCA comparison is possible. The PKNCA
  section therefore validates the published *exposure ratios* from Figure 4 and
  the Discussion, which are the only quantitative exposure predictions the paper
  makes.
* Figure 4's published values are geometric means over 1000 simulated trials of
  220 individuals each, carrying the full random-effect distribution. The
  replication uses deterministic typical-value profiles, which is why agreement
  is close but not exact for the scenarios whose effect is nonlinear in the
  exposure metric.
* Supplement Tables S2 and S3 record the base- and final-model development steps
  (OFV changes for each forward-addition and backward-elimination step). They
  contain no final parameter estimates and were not used. Figures S1 to S6
  (external validation of the prior model, IIV distributions, and stratified
  VPCs) are likewise diagnostic rather than parameter-bearing.
* The virtual cohorts fix every covariate at the reference or scenario value
  rather than sampling the joint covariate distribution of the analysis
  population, which Maleki 2024 does not publish. This is appropriate here
  because the validation targets are ratios defined against a fixed reference
  individual.
