SAF-189s (Liu 2024)

library(nlmixr2lib)
library(rxode2)
library(PKNCA)
library(dplyr)
library(ggplot2)

Liu Y, Tan Y, Hu L, Li J, Yang J, Diao L, Yang J. Population pharmacokinetics and exposure-response analyses of SAF-189s in Chinese patients with ALK+/ROS1+ non-small cell lung cancer. Front Pharmacol. 2024;15:1418549. doi:10.3389/fphar.2024.1418549

SAF-189s is an investigational second-generation ALK/ROS1 tyrosine kinase inhibitor developed from the ceritinib scaffold. This paper pools a food-effect study in healthy Chinese volunteers with a first-in-human phase I/II study in Chinese patients to build one population PK model, then uses the individual exposures from that model to fit six independent exposure-response regressions. All seven models are packaged here.

Model	Endpoint	Exposure driver
Liu_2024_saf189s	plasma concentration	(the PK model itself)
Liu_2024_saf189s_orr	overall response (IRC)	CTROUGH (Cmin,ss), log scale
Liu_2024_saf189s_pfs	progression-free survival	CTROUGH, linear
Liu_2024_saf189s_dor	duration of response	CTROUGH, linear
Liu_2024_saf189s_hyperglycemia	hyperglycemia, any grade	AUC_SAF189S, log scale
Liu_2024_saf189s_hyperglycemia_grade2	hyperglycemia, grade >= 2	AUC_SAF189S, log scale
Liu_2024_saf189s_proteinuria	proteinuria, any grade	AUC_SAF189S, log scale

pk    <- rxode2::rxode(readModelDb("Liu_2024_saf189s"))
orr   <- rxode2::rxode(readModelDb("Liu_2024_saf189s_orr"))
pfs   <- rxode2::rxode(readModelDb("Liu_2024_saf189s_pfs"))
dor   <- rxode2::rxode(readModelDb("Liu_2024_saf189s_dor"))
hyper <- rxode2::rxode(readModelDb("Liu_2024_saf189s_hyperglycemia"))
hyp2  <- rxode2::rxode(readModelDb("Liu_2024_saf189s_hyperglycemia_grade2"))
prot  <- rxode2::rxode(readModelDb("Liu_2024_saf189s_proteinuria"))

Population

pop <- pk$population
flatten_field <- function(x) {
  nms <- names(x)
  if (is.null(nms)) {
    paste(as.character(x), collapse = "; ")
  } else {
    paste(paste0(nms, ": ", as.character(x)), collapse = "; ")
  }
}
tibble::tibble(
  Field = names(pop),
  Value = vapply(pop, flatten_field, character(1))
) |>
  knitr::kable(caption = "Population metadata of the SAF-189s population PK model.")

Population metadata of the SAF-189s population PK model.

Field	Value
species	human
n_subjects	317
n_studies	2
n_observations	3,173 measurable plasma concentrations retained of 3,538 acquired (89.68%); 329 pre-dose BQL, 6 pre-dose non-BQL, 30 post-dose BQL samples and 6 subjects with no samples were excluded
age_range	18.0-84.0 years (pooled; mean 51.6, median 53). Healthy volunteers median 27 (18.0-45.0); SAF001 phase I median 51 (28.0-68.0); SAF001 phase II median 54.1 (20.0-84.0)
age_median	53 years
weight_range	37.3-92.5 kg
weight_median	63.2 kg
sex_female_pct	48.2
race_ethnicity	Asian: 100
disease_state	ALK-positive or ROS1-positive advanced non-small cell lung cancer (190 ALK+ only, 103 ROS1+ only), plus 24 healthy volunteers; 136 of 317 (43%) had brain metastases at enrolment and 221 (69.8%) were disease stage IV
dose_range	SAF-189s 20, 40, 80, 120, 160 or 210 mg orally once daily in 21-day cycles (SAF001, with a single-dose 3-day PK run-in); a single 160 mg oral dose in each of two crossover periods (STL31147, fed and fasted)
regions	China
renal_function	207 normal (65.2%), 94 mild (29.6%), 1 moderate (0.4%), 15 with eGFR >= 130 (4.8%)
hepatic_function	276 normal (87%), 41 mild dysfunction (13%); no moderate or severe
notes	Baseline characteristics from Liu 2024 Supplementary document 1, Supplementary Table 1 (a separate publisher file, Table1.docx). Sex is reported inconsistently: the count-backed Supplementary Table 1 gives 153 of 317 female (48.2%) while the Results narrative says 56.5%; the count is used here. Food had no meaningful effect (geometric mean ratio Cmax 109.1%, AUC0-t 105.1%), so the fed and fasted STL31147 periods are pooled without a food covariate.

The PK analysis set pools 24 healthy Chinese volunteers (study STL31147, a two-period fed/fasted crossover of a single 160 mg dose) with 293 patients from the phase I/II study SAF001 (NCT04237805, 20-210 mg once daily in 21-day cycles): 317 subjects and 3,173 measurable concentrations. Baseline characteristics are in Liu 2024 Supplementary document 1, Supplementary Table 1. Median age is 53 years (range 18.0-84.0) and median weight 63.2 kg (37.3-92.5); the cohort is entirely Chinese.

The three exposure-response analysis sets are all different from the PK set and from each other: 244 phase II patients with IRC assessment (ORR, PFS), 192 responders among those (DOR), and 296 patients with safety follow-up (the three adverse-event models).

Source trace

tibble::tribble(
  ~Quantity, ~Value, ~`Source location`,
  "V/F",                       "5,210 L",        "Table 2, theta1",
  "KA",                        "0.501 1/h",      "Table 2, theta2",
  "CL1/F (time-independent)",  "64 L/h",         "Table 2, theta3",
  "CL2/F (time-dependent)",    "40.3 L/h",       "Table 2, theta4",
  "Kout",                      "1.35 1/day",     "Table 2, theta5",
  "ALAG1",                     "0.483 h",        "Table 2, theta6",
  "Age exponent on CL/F",      "-0.314",         "Table 2, theta7; reference 53 y from Figure 4 caption",
  "ALKPOT 2,3,4 on V/F",       "0.734",          "Table 2, theta8",
  "ALKPOT 5 on V/F",           "0.917",          "Table 2, theta9",
  "ALKPOT 6 on V/F",           "0.784",          "Table 2, theta10",
  "IIV V/F, KA, CL/F, ALAG1",  "0.125, 0.279, 0.138, 0.120 (variances)", "Table 2, eta1-eta4; CV column confirms variance scale",
  "Proportional residual",     "0.0468 (variance) -> SD 0.2163", "Table 2, epsilon1",
  "Additive residual",         "0.0938 (variance) -> SD 0.3063 ng/mL", "Table 2, epsilon2",
  "CL(t) equation",            "CL = CL1 + CL2*exp(-Kout*(DAY-1))", "Equation 2; arm roles confirmed by Figure 1",
  "One-compartment, lagged first-order absorption", "structure", "Figure 1; Results, PopPK analysis",
  "ORR slope",                 "OR 0.556 per e-fold Cmin,ss",  "Results, Exposure-response analysis",
  "ORR intercept",             "3.871 (digitised)",            "Figure 6 fitted curve -- NOT printed",
  "PFS hazard coefficient",    "0.0071 per ng/mL",             "Table 4, PFS row",
  "DOR hazard coefficient",    "0.0059 per ng/mL",             "Table 4, DOR row",
  "Hyperglycemia slope",       "OR 3.521 per e-fold AUCss",    "Results, Exposure-safety analysis",
  "Hyperglycemia intercept",   "-9.430 (digitised)",           "Figure 10A fitted curve -- NOT printed",
  "Hyperglycemia gr>=2 slope", "OR 7.662 per e-fold AUCss",    "Results, Exposure-safety analysis",
  "Hyperglycemia gr>=2 intercept", "-16.812 (digitised)",      "Figure 10B fitted curve -- NOT printed",
  "Proteinuria slope",         "OR 2.031 per e-fold AUCss",    "Results, Exposure-safety analysis",
  "Proteinuria intercept",     "-6.297 (digitised)",           "Figure 10C fitted curve -- NOT printed"
) |>
  knitr::kable(caption = "Source location for every model equation and ini() parameter.")

Source location for every model equation and ini() parameter.

Quantity	Value	Source location
V/F	5,210 L	Table 2, theta1
KA	0.501 1/h	Table 2, theta2
CL1/F (time-independent)	64 L/h	Table 2, theta3
CL2/F (time-dependent)	40.3 L/h	Table 2, theta4
Kout	1.35 1/day	Table 2, theta5
ALAG1	0.483 h	Table 2, theta6
Age exponent on CL/F	-0.314	Table 2, theta7; reference 53 y from Figure 4 caption
ALKPOT 2,3,4 on V/F	0.734	Table 2, theta8
ALKPOT 5 on V/F	0.917	Table 2, theta9
ALKPOT 6 on V/F	0.784	Table 2, theta10
IIV V/F, KA, CL/F, ALAG1	0.125, 0.279, 0.138, 0.120 (variances)	Table 2, eta1-eta4; CV column confirms variance scale
Proportional residual	0.0468 (variance) -> SD 0.2163	Table 2, epsilon1
Additive residual	0.0938 (variance) -> SD 0.3063 ng/mL	Table 2, epsilon2
CL(t) equation	CL = CL1 + CL2exp(-Kout(DAY-1))	Equation 2; arm roles confirmed by Figure 1
One-compartment, lagged first-order absorption	structure	Figure 1; Results, PopPK analysis
ORR slope	OR 0.556 per e-fold Cmin,ss	Results, Exposure-response analysis
ORR intercept	3.871 (digitised)	Figure 6 fitted curve – NOT printed
PFS hazard coefficient	0.0071 per ng/mL	Table 4, PFS row
DOR hazard coefficient	0.0059 per ng/mL	Table 4, DOR row
Hyperglycemia slope	OR 3.521 per e-fold AUCss	Results, Exposure-safety analysis
Hyperglycemia intercept	-9.430 (digitised)	Figure 10A fitted curve – NOT printed
Hyperglycemia gr>=2 slope	OR 7.662 per e-fold AUCss	Results, Exposure-safety analysis
Hyperglycemia gr>=2 intercept	-16.812 (digitised)	Figure 10B fitted curve – NOT printed
Proteinuria slope	OR 2.031 per e-fold AUCss	Results, Exposure-safety analysis
Proteinuria intercept	-6.297 (digitised)	Figure 10C fitted curve – NOT printed

Structural checks against closed forms

The time-dependent clearance is the one piece of structure worth checking before anything stochastic. Liu 2024 Equation 2 is

$$\mathrm{CL}_i = \mathrm{CL_1/F} + \mathrm{CL_2/F}\cdot e^{-K_{out}(\mathrm{DAY}-1)}$$

with DAY the integer post-administration day, so CL/F starts at 64 + 40.3 = 104.3 L/h on day 1 and decays stepwise to 64 L/h. These are deterministic identities, so they are asserted tightly.

ref_cov <- data.frame(AGE = 53, PRIOR_ALKI = 0, TUM_ALK_MUT = 1, DIS_HEALTHY = 0)
pk_typ  <- rxode2::zeroRe(pk)

ev_sd <- rxode2::et(amt = 160, cmt = "depot") |>
  rxode2::et(seq(0, 336, by = 0.25))
sd <- rxode2::rxSolve(pk_typ, ev_sd, ref_cov, returnType = "data.frame")
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalka', 'etalcl', 'etaltlag'

# Half-life from a log-linear slope well after Tmax (pattern 11: a
# time-to-50% read includes the absorption transient and reads long).
half_life <- function(d) {
  fit <- stats::lm(log(Cc) ~ time, data = d)
  log(2) / (-stats::coef(fit)[["time"]])
}
t_half_ss  <- half_life(subset(sd, time >= 200 & time <= 336))

cl_at <- function(tt) sd$cl[which.min(abs(sd$time - tt))]
cl_expected <- function(day_index) 64 + 40.3 * exp(-1.35 * day_index)

structural <- tibble::tibble(
  Quantity = c(
    "CL/F on day 1 (L/h)", "CL/F on day 2 (L/h)", "CL/F on day 3 (L/h)",
    "CL/F on day 5 (L/h)", "Asymptotic half-life (h)"
  ),
  Model = c(cl_at(0), cl_at(24), cl_at(48), cl_at(96), t_half_ss),
  `Closed form` = c(
    cl_expected(0), cl_expected(1), cl_expected(2), cl_expected(4),
    log(2) * 5210 / 64
  ),
  Paper = c(104.3, NA, NA, NA, 56.4)
)
knitr::kable(structural, digits = 3,
             caption = "Time-dependent clearance and the steady-state half-life against their closed forms.")

Time-dependent clearance and the steady-state half-life against their closed forms.

Quantity	Model	Closed form	Paper
CL/F on day 1 (L/h)	104.300	104.300	104.3
CL/F on day 2 (L/h)	74.447	74.447	NA
CL/F on day 3 (L/h)	66.708	66.708	NA
CL/F on day 5 (L/h)	64.182	64.182	NA
Asymptotic half-life (h)	56.426	56.427	56.4

# Deterministic: the solved CL must equal Equation 2 to solver precision.
stopifnot(
  max(abs(structural$Model[1:4] - structural$`Closed form`[1:4])) < 1e-6,
  abs(t_half_ss - log(2) * 5210 / 64) < 0.05,
  # The paper prints 56.4 h; ln(2)*5210/64 = 56.43 h.
  abs(t_half_ss - 56.4) < 0.5
)

The day-1 half-life the paper quotes (34.6 h) is log(2) * 5210 / 104.3 = 34.62 h. A regression over 12-23 h reads slightly longer because absorption has not fully died away inside the first dosing interval, so the clean check is the closed form above rather than a slope fit on day 1.

Steady-state exposure has an equally direct identity: at steady state CL/F has reached its asymptote, so AUC over one dosing interval is exactly dose / CL1.

ev_md <- rxode2::et(amt = 160, cmt = "depot", ii = 24, until = 24 * 20) |>
  rxode2::et(seq(0, 24 * 21, by = 0.1))
md <- rxode2::rxSolve(pk_typ, ev_md, ref_cov, returnType = "data.frame",
                      maxsteps = 200000L)
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalka', 'etalcl', 'etaltlag'

last_int <- subset(md, time >= 24 * 20 & time <= 24 * 21)
auc_trap <- sum(diff(last_int$time) *
                  (head(last_int$Cc, -1) + tail(last_int$Cc, -1)) / 2)
auc_closed <- 160 / 64 * 1000   # mg / (L/h) -> mg*h/L -> ng*h/mL

c(`AUCss simulated (ng*h/mL)` = auc_trap,
  `Dose / CL1 (ng*h/mL)` = auc_closed,
  `percent difference` = 100 * (auc_trap - auc_closed) / auc_closed) |>
  round(2)
#> AUCss simulated (ng*h/mL)      Dose / CL1 (ng*h/mL)        percent difference 
#>                   2494.13                   2500.00                     -0.23

# Same drawn parameters on both sides -> pure numerical error, so a tight
# bound is correct here (see CLAUDE.md on cohort-extreme assertions).
stopifnot(abs(100 * (auc_trap - auc_closed) / auc_closed) < 0.5)

Covariate effects: replicating Figure 4

Liu 2024 Figure 4 is a forest plot of the covariate effects on CL/F, V/F and AUCss relative to the reference subject (median age 53 years, ALK-positive patient with no prior ALK-inhibitor treatment). The caption prints the point estimates, and every one of them is a deterministic function of the ini() values.

cl_ratio <- function(age) (age / 53)^(-0.314)
vc_ratio <- function(alki = 0, nonalk = 0, healthy = 0) {
  0.734^alki * 0.917^nonalk * 0.784^healthy
}

fig4 <- tibble::tribble(
  ~Comparison,                        ~Parameter, ~Model,                         ~Paper,
  "Age 25 y vs 53 y",                 "CL/F",     cl_ratio(25),                    1.27,
  "Age 71 y vs 53 y",                 "CL/F",     cl_ratio(71),                    0.912,
  "Age 25 y vs 53 y",                 "AUCss",    1 / cl_ratio(25),                0.79,
  "Age 71 y vs 53 y",                 "AUCss",    1 / cl_ratio(71),                1.10,
  "ALKPOT 2,3,4 vs 1",                "V/F",      vc_ratio(alki = 1),              0.734,
  "ALKPOT 5 vs 1",                    "V/F",      vc_ratio(nonalk = 1),            0.917,
  "ALKPOT 6 vs 1",                    "V/F",      vc_ratio(healthy = 1),           0.784
) |>
  mutate(`% difference` = 100 * (Model - Paper) / Paper)

knitr::kable(fig4, digits = 3,
             caption = "Figure 4 covariate point estimates (ratio to the reference subject).")

Figure 4 covariate point estimates (ratio to the reference subject).

Comparison	Parameter	Model	Paper	% difference
Age 25 y vs 53 y	CL/F	1.266	1.270	-0.307
Age 71 y vs 53 y	CL/F	0.912	0.912	0.031
Age 25 y vs 53 y	AUCss	0.790	0.790	-0.022
Age 71 y vs 53 y	AUCss	1.096	1.100	-0.349
ALKPOT 2,3,4 vs 1	V/F	0.734	0.734	0.000
ALKPOT 5 vs 1	V/F	0.917	0.917	0.000
ALKPOT 6 vs 1	V/F	0.784	0.784	0.000

# All deterministic; the paper's own rounding is the only slack.
stopifnot(max(abs(fig4$`% difference`)) < 2)

The paper states these as percentage changes: CL/F “increased by 27% in the 5th percentile (25 years) and decreased by 8.8% in the 95th percentile (71 years)”, AUC “decreased by 21% and increased by 10%”, and V/F “26.6% lower when ALKPOT = 2, 3, 4; 8.3% lower when ALKPOT = 5; and 21.6% lower when ALKPOT = 6”. Every one reproduces.

Note that the age effect on AUCss is the reciprocal of its effect on CL/F, because AUCss = dose / CL/F and the age term multiplies total clearance. That is what pins the covariate to total CL rather than to the CL1 arm alone: if the exponent applied only to CL1, the day-1 CL/F ratio at age 25 would be 1.163 rather than 1.266, and the paper’s +27% could not be recovered.

Virtual cohort and simulation

The four phase II dose levels are simulated with the patient covariate mix from Supplementary Table 1. Healthy volunteers are simulated separately as the single-dose 160 mg design of STL31147.

n_per_arm <- 150L
doses <- c(80, 120, 160, 210)

# Age: median 53, range 18-84 (Supplementary Table 1). A truncated normal
# centred at the reported mean reproduces the median and the range.
draw_age <- function(n) {
  a <- stats::rnorm(n, mean = 51.6, sd = 11.5)
  pmin(pmax(a, 18), 84)
}

# ALKPOT mix among the 293 patients (Supplementary Table 1):
#   level 1     115  (ALK+, ALK-inhibitor naive)   -> reference
#   levels 2-4   64  (ALK+, prior ALK inhibitor)
#   level 5     114  (ROS1+ or unknown)
draw_alkpot <- function(n) {
  sample(c("naive", "pretreated", "ros1_or_unknown"), n, replace = TRUE,
         prob = c(115, 64, 114) / 293)
}

# One covariate template, replicated across the four arms. Subject j of every
# arm therefore has the same age and the same ALKPOT level, and (with the
# common random numbers set in the solve below) the same etas -- so the arms
# differ ONLY by dose. That is what turns the dose-proportionality check into
# a structural identity rather than a race between four noisy medians.
set.seed(20240716)
base_grp <- draw_alkpot(n_per_arm)
base_age <- draw_age(n_per_arm)

cohort <- lapply(seq_along(doses), function(k) {
  data.frame(
    id          = (k - 1L) * n_per_arm + seq_len(n_per_arm),
    subject     = seq_len(n_per_arm),
    dose        = doses[k],
    AGE         = base_age,
    PRIOR_ALKI  = as.integer(base_grp == "pretreated"),
    TUM_ALK_MUT = as.integer(base_grp %in% c("naive", "pretreated")),
    DIS_HEALTHY = 0L
  )
}) |>
  bind_rows() |>
  mutate(treatment = factor(paste0(dose, " mg QD"),
                            levels = paste0(doses, " mg QD")))

summary_cohort <- cohort |>
  group_by(treatment) |>
  summarise(
    n = n(),
    `median age` = round(stats::median(AGE), 1),
    `prior ALKi (%)` = round(100 * mean(PRIOR_ALKI), 1),
    `ALK+ (%)` = round(100 * mean(TUM_ALK_MUT), 1),
    .groups = "drop"
  )
knitr::kable(summary_cohort, caption = "Simulated cohort by dose arm. The arms are covariate-matched by construction, so the summaries are identical; only the dose differs.")

Simulated cohort by dose arm. The arms are covariate-matched by construction, so the summaries are identical; only the dose differs.

treatment	n	median age	prior ALKi (%)	ALK+ (%)
80 mg QD	150	55.6	25.3	56.7
120 mg QD	150	55.6	25.3	56.7
160 mg QD	150	55.6	25.3	56.7
210 mg QD	150	55.6	25.3	56.7

rxSolve() on an rxUi scales super-linearly in the number of subjects, so each arm is solved in its own call. rxSetSeed() is called inside the loop so every arm draws the same eta sequence: with common random numbers the dose arms differ only by dose, which is what makes the dose-proportionality check below meaningful.

tau <- 24
n_days <- 21L

obs_times <- sort(unique(c(
  seq(0, tau, by = 0.25),                                   # first interval
  seq(tau, (n_days - 1) * tau, by = tau),                   # daily troughs
  seq((n_days - 1) * tau, n_days * tau, by = 0.25)          # steady-state interval
)))

sim <- lapply(doses, function(d) {
  cov_d <- cohort |>
    filter(dose == d) |>
    select(id, AGE, PRIOR_ALKI, TUM_ALK_MUT, DIS_HEALTHY)

  ev <- rxode2::et(amt = d, cmt = "depot", ii = tau,
                   until = (n_days - 1) * tau,
                   id = cov_d$id) |>
    rxode2::et(obs_times, cmt = "central", id = cov_d$id)

  rxode2::rxSetSeed(4242)
  rxode2::rxSolve(pk, ev, cov_d, returnType = "data.frame",
                  maxsteps = 200000L, addDosing = FALSE) |>
    mutate(dose = d)
}) |>
  bind_rows() |>
  left_join(select(cohort, id, treatment, AGE), by = "id")

nrow(sim)
#> [1] 127200

band <- sim |>
  filter(time <= tau | time >= (n_days - 1) * tau) |>
  mutate(window = ifelse(time <= tau, "Day 1", "Day 21 (steady state)"),
         t_in_int = ifelse(time <= tau, time, time - (n_days - 1) * tau)) |>
  group_by(treatment, window, t_in_int) |>
  summarise(lo = stats::quantile(Cc, 0.05),
            md = stats::median(Cc),
            hi = stats::quantile(Cc, 0.95), .groups = "drop")

ggplot(band, aes(t_in_int, md, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, colour = NA) +
  geom_line() +
  facet_wrap(~window) +
  labs(x = "Time within dosing interval (h)", y = "SAF-189s (ng/mL)",
       colour = NULL, fill = NULL) +
  theme_bw()

Simulated SAF-189s plasma concentrations over the first and the twenty-first dosing interval. Ribbons are the 5th-95th percentiles, lines the median.

Accumulation between day 1 and day 21 is the visible signature of the time-dependent clearance working in the opposite direction to accumulation: clearance rises from 64 L/h effective on day 1 (as 104.3 L/h) toward 64 L/h, which damps the accumulation that the 56.4 h terminal half-life would otherwise produce.

PKNCA validation

start_ss <- (n_days - 1) * tau
end_ss   <- n_days * tau

conc_df <- sim |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, treatment)

dose_df <- cohort |>
  transmute(id, time = start_ss, amt = dose, treatment)

conc_obj <- PKNCA::PKNCAconc(conc_df, Cc ~ time | treatment + id,
                             concu = "ng/mL", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg")

intervals <- data.frame(
  start   = start_ss,
  end     = end_ss,
  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) |>
  select(treatment, id, PPTESTCD, PPORRES)

nca_wide |>
  group_by(treatment, PPTESTCD) |>
  summarise(median = stats::median(PPORRES), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = median) |>
  knitr::kable(digits = 1,
               caption = "Median steady-state NCA parameters by dose arm (day 21 interval).")

Median steady-state NCA parameters by dose arm (day 21 interval).

treatment	auclast	cav	cmax	cmin	tmax
80 mg QD	1251.8	52.2	58.5	44.4	5.2
120 mg QD	1877.6	78.2	87.7	66.5	5.2
160 mg QD	2503.5	104.3	117.0	88.7	5.2
210 mg QD	3285.9	136.9	153.6	116.5	5.2

Comparison against the published values

Liu 2024 does not tabulate a full NCA table, but it prints four steady-state quantities that a simulation must reproduce: the geometric mean AUCss at the recommended 160 mg dose (2,374 ngh/mL), the median AUCss across the safety analysis set (2,233 ngh/mL), the median Cmin,ss of the phase II efficacy set (79 ng/mL, the Table 3 quartile boundary), and the median Tmax (6 h, from the food-effect study).

sim_nca <- nca_wide |>
  filter(PPTESTCD %in% c("auclast", "cmin", "tmax")) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

published <- tibble::tribble(
  ~treatment,   ~auclast, ~cmin, ~tmax,
  "160 mg QD",  2374,     NA,    6
) |>
  mutate(treatment = factor(treatment, levels = levels(cohort$treatment)))

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = filter(sim_nca, treatment == "160 mg QD"),
  reference     = filter(published, !is.na(auclast)),
  by            = "treatment",
  params        = c("auclast", "tmax"),
  units         = c(auclast = "ng*h/mL", tmax = "h"),
  tolerance_pct = 20
)
knitr::kable(cmp,
             caption = "Simulated vs published steady-state exposure at the recommended 160 mg once-daily dose. * marks rows differing by >20%.")

Simulated vs published steady-state exposure at the recommended 160 mg once-daily dose. * marks rows differing by >20%.

NCA parameter	treatment	Reference	Simulated	% diff
Tmax (h)	160 mg QD	6	5.25	-12.5%
AUClast (ng*h/mL)	160 mg QD	2370	2500	+5.5%

auc160  <- filter(sim_nca, treatment == "160 mg QD")$auclast
cmin_all <- sim_nca$cmin

# Dose proportionality, matched subject by subject. The arms are
# covariate-matched and share common random numbers, so subject j at 160 mg
# and subject j at 80 mg differ only by dose; a linear-elimination model must
# give a ratio of exactly 2 for every one of them.
prop_ratio <- sim_nca |>
  left_join(select(cohort, id, subject), by = "id") |>
  filter(treatment %in% c("80 mg QD", "160 mg QD")) |>
  select(subject, treatment, auclast) |>
  tidyr::pivot_wider(names_from = treatment, values_from = auclast) |>
  mutate(ratio = `160 mg QD` / `80 mg QD`)

gates <- tibble::tibble(
  Claim = c(
    "Geometric mean AUCss at 160 mg near the published 2,374 ng*h/mL",
    "Median Cmin,ss across 80-210 mg near the published 79 ng/mL",
    "AUCss is dose proportional (worst per-subject 160 mg / 80 mg ratio)"
  ),
  Achieved = c(
    exp(mean(log(auc160))),
    stats::median(cmin_all),
    prop_ratio$ratio[which.max(abs(prop_ratio$ratio - 2))]
  ),
  Target = c(2374, 79, 2)
)
gates$`% difference` <- 100 * (gates$Achieved - gates$Target) / gates$Target
knitr::kable(gates, digits = 2, caption = "Cohort-level exposure gates.")

Cohort-level exposure gates.

Claim	Achieved	Target	% difference
Geometric mean AUCss at 160 mg near the published 2,374 ng*h/mL	2451.33	2374	3.26
Median Cmin,ss across 80-210 mg near the published 79 ng/mL	73.55	79	-6.90
AUCss is dose proportional (worst per-subject 160 mg / 80 mg ratio)	2.00	2	0.00

stopifnot(
  # Centre, not extremes. Realised +2.2% and -9.3%. With CV ~38% on clearance
  # over 150 subjects the cohort statistic itself carries a couple of percent
  # of draw-to-draw noise (and rxode2's stream is partitioned per solver
  # thread, so CI draws a different cohort). These bounds keep several times
  # that headroom while still going red on a mis-transcribed dose, volume or
  # clearance, each of which moves the whole distribution by tens of percent.
  abs(gates$`% difference`[1]) < 20,
  abs(gates$`% difference`[2]) < 30,
  # Dose proportionality is a structural identity of a linear-elimination
  # model with covariate-matched, common-random-number arms: it holds for
  # EVERY subject, up to trapezoidal error on a shared time grid. Asserted on
  # the worst subject, not the median.
  abs(gates$`% difference`[3]) < 0.1
)

The median Cmin,ss falls a little below the paper’s 79 ng/mL because the simulated cohort splits evenly across the four dose arms, whereas actual phase II enrolment was concentrated at 160 mg and above (56.8% of phase II subjects received 160 mg and a further 10% received 210 mg), which pulls the real cohort’s median up. The dose-proportionality check is the tight one: elimination is linear in this model and the arms are covariate-matched with common random numbers, so doubling the dose must double AUCss for every subject.

Exposure-response

Efficacy: Figure 6 and Table 4

The exposure-efficacy models are static regressions on CTROUGH, so they are driven directly by the simulated steady-state troughs.

# These models carry no ODE and no dose. The exposure column is supplied as a
# time-varying covariate on the event table, one observation row per grid
# point, which is the idiom the other static exposure-response vignettes use.
er_eval <- function(model, output, cov_name, grid) {
  ev <- data.frame(id = 1L, time = seq_along(grid), amt = 0, evid = 0L)
  ev[[cov_name]] <- grid
  s <- as.data.frame(rxode2::rxSolve(model, events = ev,
                                     returnType = "data.frame"))
  stopifnot(nrow(s) == length(grid))
  s[[output]]
}

ctrough_grid <- exp(seq(log(18), log(182), length.out = 200))
orr_curve <- data.frame(
  CTROUGH = ctrough_grid,
  prob_orr_central = er_eval(orr, "prob_orr_central", "CTROUGH", ctrough_grid)
)

quartiles <- c(18, 60.7, 79, 108, 182)
observed_q <- tibble::tibble(
  # Liu 2024 Table 3: CR + PR counts per Cmin,ss quartile, n = 61 each.
  mid = sqrt(head(quartiles, -1) * tail(quartiles, -1)),
  p   = c(54, 44, 47, 47) / 61
)

ggplot(orr_curve, aes(CTROUGH, prob_orr_central)) +
  geom_line(colour = "blue") +
  geom_point(data = observed_q, aes(mid, p), inherit.aes = FALSE,
             shape = 2, size = 3) +
  geom_vline(xintercept = quartiles, linetype = "dashed", alpha = 0.4) +
  scale_x_log10() +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "Cmin,ss (ng/mL, log scale)",
       y = "Probability of IRC-assessed response",
       caption = "Replicates Figure 6 of Liu 2024; triangles are the observed per-quartile response rates from Table 3.") +
  theme_bw()

ui_orr <- orr$iniDf
slope_orr <- ui_orr$est[ui_orr$name == "e_ctrough_logit"]
p_at_median <- orr_curve$prob_orr_central[which.min(abs(orr_curve$CTROUGH - 79))]

hazard_at <- function(mod, ctrough) {
  er_eval(mod, "hr", "CTROUGH", ctrough)
}

er_eff <- tibble::tibble(
  Claim = c(
    "ORR odds ratio per e-fold Cmin,ss",
    "P(response) at the median Cmin,ss (79 ng/mL) vs the observed ORR",
    "PFS hazard ratio per 1 ng/mL Cmin,ss",
    "DOR hazard ratio per 1 ng/mL Cmin,ss"
  ),
  Achieved = c(
    exp(slope_orr),
    p_at_median,
    hazard_at(pfs, 1) / hazard_at(pfs, 0),
    hazard_at(dor, 1) / hazard_at(dor, 0)
  ),
  Published = c(0.556, 0.7869, 1.007, 1.006)
)
er_eff$`% difference` <- 100 * (er_eff$Achieved - er_eff$Published) / er_eff$Published
knitr::kable(er_eff, digits = 4,
             caption = "Exposure-efficacy models against the published estimates.")

Exposure-efficacy models against the published estimates.

Claim	Achieved	Published	% difference
ORR odds ratio per e-fold Cmin,ss	0.5560	0.5560	0.0000
P(response) at the median Cmin,ss (79 ng/mL) vs the observed ORR	0.7871	0.7869	0.0285
PFS hazard ratio per 1 ng/mL Cmin,ss	1.0071	1.0070	0.0124
DOR hazard ratio per 1 ng/mL Cmin,ss	1.0059	1.0060	-0.0082

stopifnot(
  # Odds and hazard ratios are exp() of a stored coefficient: deterministic.
  abs(er_eff$`% difference`[c(1, 3, 4)]) < 0.1,
  # The intercept was digitised; agreement with the independently reported
  # observed ORR is the check that it was read correctly.
  abs(er_eff$Achieved[2] - 0.7869) < 0.01
)

The response curve is essentially flat across the observed exposure range, which is the paper’s finding: the odds ratio confidence interval (0.285-1.083) spans 1 and p = 0.084. The two Cox models are equally null (p = 0.059 and p = 0.251).

Safety: Figure 10

auc_grid <- exp(seq(log(174.0126), log(8338.102), length.out = 300))

safety_models <- list(
  "Hyperglycemia (any grade)"  = list(m = hyper, out = "prob_hyperglycemia"),
  "Hyperglycemia (grade >= 2)" = list(m = hyp2,  out = "prob_hyperglycemia_grade2"),
  "Proteinuria (any grade)"    = list(m = prot,  out = "prob_proteinuria")
)

safety_curve <- lapply(names(safety_models), function(nm) {
  s <- safety_models[[nm]]
  data.frame(AUC = auc_grid, Endpoint = nm,
             p = er_eval(s$m, s$out, "AUC_SAF189S", auc_grid))
}) |>
  bind_rows() |>
  mutate(Endpoint = factor(Endpoint, levels = names(safety_models)))

ggplot(safety_curve, aes(AUC, p)) +
  geom_line(colour = "blue") +
  facet_wrap(~Endpoint) +
  geom_vline(xintercept = 2233, linetype = "dashed", alpha = 0.5) +
  scale_x_log10() +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "AUCss (ng*h/mL, log scale)", y = "Probability of AE",
       caption = "Replicates Figure 10 of Liu 2024. Dashed line is the published median AUCss of 2,233 ng*h/mL.") +
  theme_bw()

p_at <- function(df, endpoint, auc) {
  d <- filter(df, Endpoint == endpoint)
  d$p[which.min(abs(d$AUC - auc))]
}
slope_of <- function(mod) {
  d <- mod$iniDf
  d$est[d$name == "e_auc_saf189s_logit"]
}

er_saf <- tibble::tibble(
  Endpoint = c("Hyperglycemia (any grade)", "Hyperglycemia (grade >= 2)",
               "Proteinuria (any grade)"),
  `OR (model)` = c(exp(slope_of(hyper)), exp(slope_of(hyp2)), exp(slope_of(prot))),
  `OR (paper)` = c(3.521, 7.662, 2.031),
  `P at median AUCss` = c(
    p_at(safety_curve, "Hyperglycemia (any grade)", 2233),
    p_at(safety_curve, "Hyperglycemia (grade >= 2)", 2233),
    p_at(safety_curve, "Proteinuria (any grade)", 2233)
  ),
  `Observed incidence (Table 5)` = c(165 / 296, 82 / 296, 90 / 296)
)
knitr::kable(er_saf, digits = 4,
             caption = "Exposure-safety models against the published odds ratios and the Table 5 incidences.")

Exposure-safety models against the published odds ratios and the Table 5 incidences.

Endpoint	OR (model)	OR (paper)	P at median AUCss	Observed incidence (Table 5)
Hyperglycemia (any grade)	3.521	3.521	0.5679	0.5574
Hyperglycemia (grade >= 2)	7.662	7.662	0.2469	0.2770
Proteinuria (any grade)	2.031	2.031	0.3025	0.3041

stopifnot(
  # Odds ratios are exp() of a stored coefficient: deterministic.
  max(abs(100 * (er_saf$`OR (model)` - er_saf$`OR (paper)`) /
            er_saf$`OR (paper)`)) < 0.1,
  # The three intercepts were digitised from Figure 10. The check that they
  # were read correctly is that each model's fitted probability at the median
  # exposure lands on the independently counted Table 5 incidence. All three
  # agree within 3 percentage points; the bound admits that without becoming
  # un-failable (a mis-read intercept moves these by tens of points).
  max(abs(er_saf$`P at median AUCss` -
            er_saf$`Observed incidence (Table 5)`)) < 0.05
)

Finally, the three safety models can be driven by the simulated exposures rather than by a grid, which is how a dose-selection analysis would use them.

auc_by_subject <- nca_wide |>
  filter(PPTESTCD == "auclast") |>
  transmute(id, treatment, AUC_SAF189S = PPORRES)

pred_by_dose <- lapply(levels(auc_by_subject$treatment), function(tr) {
  d <- filter(auc_by_subject, treatment == tr)
  data.frame(
    treatment = tr,
    hyper  = mean(er_eval(hyper, "prob_hyperglycemia", "AUC_SAF189S", d$AUC_SAF189S)),
    hyper2 = mean(er_eval(hyp2, "prob_hyperglycemia_grade2", "AUC_SAF189S", d$AUC_SAF189S)),
    prot   = mean(er_eval(prot, "prob_proteinuria", "AUC_SAF189S", d$AUC_SAF189S))
  )
}) |>
  bind_rows()

pred_by_dose |>
  rename("Dose arm" = treatment,
         "Hyperglycemia, any grade" = hyper,
         "Hyperglycemia, grade >= 2" = hyper2,
         "Proteinuria, any grade" = prot) |>
  knitr::kable(digits = 3,
               caption = "Model-predicted mean adverse-event probability by dose arm, using the simulated steady-state AUCss of each subject.")

Model-predicted mean adverse-event probability by dose arm, using the simulated steady-state AUCss of each subject.

Dose arm	Hyperglycemia, any grade	Hyperglycemia, grade >= 2	Proteinuria, any grade
80 mg QD	0.388	0.106	0.224
120 mg QD	0.508	0.204	0.277
160 mg QD	0.593	0.305	0.320
210 mg QD	0.669	0.419	0.362

stopifnot(
  # Trend, not step-by-step monotonicity of a noisy cohort statistic.
  pred_by_dose$hyper[nrow(pred_by_dose)] > pred_by_dose$hyper[1],
  pred_by_dose$hyper2[nrow(pred_by_dose)] > pred_by_dose$hyper2[1]
)

This reproduces the paper’s dose-selection argument quantitatively: predicted grade >= 2 hyperglycemia rises steeply from the 80 mg arm to the 210 mg arm, while the 160 mg arm sits close to the observed cohort-wide rate. Liu 2024 selected 160 mg once daily as the recommended phase III dose on exactly this basis, noting that the geometric mean AUCss at 160 mg (2,374 ngh/mL) is close to the median exposure (2,233 ngh/mL) below which adverse-event incidence was significantly lower.

Assumptions and deviations

Figure-derived intercepts (four models). Liu 2024 reports the slope of every exposure-response regression but no intercept, and prints no coefficient table for any of them. The intercepts of Liu_2024_saf189s_orr, Liu_2024_saf189s_hyperglycemia, Liu_2024_saf189s_hyperglycemia_grade2 and Liu_2024_saf189s_proteinuria were therefore recovered by digitising the fitted curves of Figures 6 and 10. In each case the slope was digitised too and reproduced the printed odds ratio to within 0.15% (0.5553 vs 0.556; 3.524 vs 3.521; 7.672 vs 7.662; 2.032 vs 2.031), with the logit linear in log exposure to R-squared >= 0.99995 – which is what licenses reading the intercept off the same curve. Each intercept was then cross-checked against a source that is not the figure: the fitted probability at the median exposure against the independently counted event proportion (0.787 vs the reported ORR 78.69%; 0.569 vs 55.74%; 0.248 vs 27.70%; 0.303 vs 30.41%). The axis calibration itself was validated by the median-quartile dashed line landing at 2,231 ngh/mL against the 2,233 ngh/mL printed in the Discussion. The slopes in the model files are the printed values; only the intercepts are figure-derived.

No baseline hazard for PFS and DOR. A Cox regression is semiparametric, so h0(t) is not an unreported parameter but a quantity the method never produces. Liu_2024_saf189s_pfs and Liu_2024_saf189s_dor therefore return the relative hazard only and offer no survivor function. Liu 2024 Figure 8 gives Kaplan-Meier curves stratified at the median exposure as an empirical description of the absolute time course; both log-rank tests were non-significant.

Tumour-growth model not extracted. Liu 2024 Equation 1 gives a Stein-Fojo biexponential model for the sum of longest lesion diameters, but reports no estimate for any of its parameters (SOD0, KG, KS) anywhere in the paper or either supplementary document, and prints no figure of them. The Results say only that “there was no trend observed between the exposure and tumor growth model parameters”. With no point estimates there is nothing to transcribe, so no model file was created for it.

Endpoints without an exposure relationship. Liu 2024 screened hypercholesterolemia, nausea, vomiting and diarrhea against exposure and found no significant trend, reporting incidences (Table 5) but no coefficients. No model files exist for those endpoints for the same reason.

Residual error entered as variances. Table 2 prints the residual terms as 0.0468 and 0.0938 without stating a scale. They are entered here as variances, so propSd = sqrt(0.0468) = 0.2163 and addSd = sqrt(0.0938) = 0.3063 ng/mL. Read as standard deviations instead, the proportional error would be 4.68% – tighter than the assay on pooled sparse phase I/II oncology data – and the additive error would be 0.094 ng/mL, a fifth of the 0.5 ng/mL lower limit of quantitation. Both readings of the variance interpretation are ordinary; both readings of the SD interpretation are not.

Interindividual variability entered as variances. Unambiguous: Table 2 prints a percent CV beside each omega-squared, and sqrt(exp(omega^2) - 1) reproduces all four printed CVs (36.5%, 56.7%, 38.5%, 35.7%) to the printed decimal.

A single eta on clearance. The Results text says an exponential IIV model was implemented “including CL1/F, CL2/F, V/F, KA, and tlag”, but Table 2 lists only four eta rows with a single clearance term. Equation 2 settles it: exp(eta_i) multiplies the whole bracket, so one eta is shared by both clearance arms. The model file applies etalcl to both, which is algebraically identical.

High shrinkage on KA and the lag time. Table 2 reports 53.3% and 66.2% shrinkage on those two etas, so their individual estimates are weakly informed. Simulated absorption-phase variability should be read with that in mind; the clearance eta (2.2% shrinkage) and the volume eta (41.1%) are better determined.

Covariates screened but not retained are recorded in the model file’s covariatesDataExcluded list rather than covariateData, because the paper reports no point estimate for any of them: body weight, sex, creatinine clearance, albumin, ALT, ALP, total bilirubin, serum creatinine and concomitant medications.

Deterministic probability outputs. The four logistic models emit a probability with a placeholder additive residual of 0.001 fixed, because rxode2 requires an observation declaration; the source likelihood is Bernoulli and estimates no residual. Downstream callers should sample outcomes with rbinom(n, 1, prob).

Cohort composition. The simulated cohort splits evenly across the four phase II dose arms, whereas actual enrolment concentrated on 160 mg (54.8% of PK-evaluable subjects). Cohort-wide exposure summaries here are therefore not directly comparable with the paper’s pooled medians; the per-arm comparison at 160 mg is.

Errata

The CL1/CL2 prose is reversed. Liu 2024 Results describes CL1/F as “the initial apparent clearance after single dosing” and CL2/F as “the subsequent time-varying apparent clearance at steady state”. Equation 2, Figure 1 and the printed half-lives all say the opposite: exp(-Kout*(DAY-1)) multiplies TVCL2, Figure 1 labels the two elimination arms “Time independent CL1” and “Time dependent CL2”, and only the assignment used here reproduces both printed half-lives (ln(2)*5210/(64+40.3) = 34.6 h after the first dose; ln(2)*5210/64 = 56.4 h at steady state). The reversed assignment reproduces neither. The parameter values are fine; the sentence is the error.

Table 2 unit labels. Table 2 prints “KA (L h-1)” and “Kout (L d-1)”. Both are first-order rate constants; the litre is spurious. The Results narrative gives them correctly as 1/h and 1/day.

The odds-ratio unit statement. The Results say the exposure-safety odds ratios correspond to “an increase in SAF-189s exposure of 1 nghmL-1d-1”. An odds ratio of 3.5 per 1 ngh/mL across a range spanning thousands of ng*h/mL is arithmetically impossible. The same paragraph states that “a linear effect of log AUCss described the data well”, and digitising Figure 10 confirms the unit is one natural-log unit (an e-fold). The models use the e-fold reading.

Sex is reported two ways. Supplementary Table 1 counts 153 of 317 subjects female (48.2%); the Results narrative says “approximately half of them were women (56.5%)”. The count-backed 48.2% is used in the population metadata. Sex was screened and not retained, so nothing in the model depends on it.

Supplementary Table 1 has misaligned laboratory rows. That table reports a median serum creatinine of 163 umol/L (range 139-190) and a median total bilirubin of 68 umol/L (range 35.0-110). Neither is credible for a cohort the same table describes as 65% renally normal and 87% hepatically normal, and the creatinine value contradicts the creatinine clearance of 96.2 mL/min in the row above it. Median height of 134 cm and median hemoglobin of 146 g/L with a range starting at 20 g/L are similarly implausible. These look like rows shifted against their labels. None of the affected covariates was retained in the final model, so the extraction does not depend on them; they are recorded in covariatesDataExcluded with the discrepancy flagged.

Kout RSE. Kout is the least precisely estimated structural parameter (RSE 31.2%, bootstrap 95% CI 0.927-2.24 per day). Its value controls only how fast CL/F reaches its asymptote (half-life ln(2)/1.35 = 0.51 days), not the asymptote itself, so steady-state exposure is insensitive to it; day 1 to day 3 predictions are not.