Skip to contents

MIT-001 (Bae 2026)

Source: vignettes/articles/Bae_2026_MIT001.Rmd
Bae_2026_MIT001.Rmd

Model and source 

#> ℹ parameter labels from comments will be replaced by 'label()'

  • Citation: Bae S, Kim E, Rhee SJ, Kim S, Yu KS, Lee S. Population Pharmacokinetic Analysis of MIT-001, a Novel Ferroptosis Inhibitor, for Dose Optimization. J Clin Pharmacol. 2026;66(4):e70189. doi:10.1002/jcph.70189

  • Description: Three-compartment population PK model for MIT-001, a mitochondria-targeting ferroptosis inhibitor, with first-order subcutaneous absorption in healthy Korean adults (Bae 2026)

  • Article: https://doi.org/10.1002/jcph.70189

  • Supplement (Supplementary Table 1, Supplementary Figure 1): retrieved from the Europe PMC open-access supplementary-files endpoint for PMC13084164.

MIT-001 is a mitochondria-targeting antioxidant and ferroptosis inhibitor. Bae 2026 pooled three phase 1 trials in healthy Korean adults to build a population PK model, then used it to pick a subcutaneous (SC) dose range for a planned trial in patients with oral mucositis (OM). The model is a three-compartment linear model with first-order SC absorption and no absorption lag.

Population

119 healthy Korean adults contributed 2859 plasma concentrations across three randomized, double-blind, placebo-controlled phase 1 trials (Bae 2026 Table 1; per-study designs and sampling schedules in Supplementary Table 1):

  • Study 1 (NCT01737424) – single ascending dose, IV. Cohort 1 received 3-200 mg as 30-min infusions; cohort 2 received 0.3-200 mg as 5-min infusions. 1185 samples.
  • Study 2 (NCT03196804) – multiple ascending dose, IV. 3-30 mg as 30-min infusions once daily for 7 days. 1038 samples.
  • Study 3 (NCT05389696) – SC dose escalation. Cohort 1 received single SC doses of 10-40 mg, and the 40 mg group additionally received a single 40 mg IV infusion after a two-week washout, which is what identifies bioavailability. Cohort 2 received 20 or 40 mg SC once daily for 7 days. 636 samples.

95 participants contributed IV data and 30 contributed SC data. The cohort was almost entirely male (3 of 119 female), aged 20-44 years, 55-87 kg, with transaminases in the healthy range (ALT 6-53 IU/L).

str(ui$population)
#> List of 14
#>  $ species       : chr "human"
#>  $ n_subjects    : num 119
#>  $ n_studies     : num 3
#>  $ n_observations: num 2859
#>  $ age_range     : chr "20-44 years"
#>  $ age_median    : chr "26, 28 and 28.5 years in studies 1, 2 and 3"
#>  $ weight_range  : chr "55-87 kg"
#>  $ weight_median : chr "67.6, 67.7 and 70.6 kg in studies 1, 2 and 3"
#>  $ sex_female_pct: num 2.5
#>  $ race_ethnicity: Named num 100
#>   ..- attr(*, "names")= chr "Asian"
#>  $ disease_state : chr "healthy adults (no disease); the model was developed to support dose selection for patients with oral mucositis"
#>  $ dose_range    : chr "IV 0.3-200 mg single dose and 3-30 mg once daily for 7 days; SC 10-40 mg single dose and 20-40 mg once daily for 7 days"
#>  $ regions       : chr "Republic of Korea"
#>  $ notes         : chr "Pooled from three phase 1 randomized, double-blind, placebo-controlled trials in healthy Korean adults: study 1"| __truncated__

Source trace

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

Equation / parameter Value Source location
Three-compartment structure, first-order SC absorption, no lag n/a Figure 1; Results “Pharmacokinetic Model”
lcl (CL, typical) 2.17 L/h Table 2, theta1
e_alt_cl (ALT exponent on CL) -0.176 Table 2, theta2
CL = 2.17 * (ALT / 16)^-0.176 n/a Table 2, footnote a
lvc (V1) 6.89 L Table 2, V1
lvp (V2, typical) 23.6 L Table 2, theta3
e_wt_vp (WT exponent on V2) 1 (structural, not estimated) Table 2, footnote b – prints no exponent
V2 = 23.6 * (WT / 68.4) n/a Table 2, footnote b
lvp2 (V3, typical) 68.3 L Table 2, theta4
e_wt_vp2 (WT exponent on V3) 0.742 Table 2, theta5
V3 = 68.3 * (WT / 68.4)^0.742 n/a Table 2, footnote c
lq (Q12) 108 L/h Table 2, Q12
lq2 (Q13) 32.9 L/h Table 2, Q13
lka (Ka) 0.757 1/h Table 2, Ka
lfdepot (F, SC) 0.82 Table 2, F
IIV on CL, V1, V2, V3, Q12, Q13, F 29.3, 49.4, 32.3, 16.2, 32.8, 26.6, 9.0 %CV Table 2, “Inter-individual variability”
propSd 0.108 Table 2, “Residual error”
Efficacy target Cavg > 87 ng/mL n/a Methods “Simulation”; mouse OM model
Safety targets Cmax < 10777.2 ng/mL, AUCtau < 16167.4 ng*h/mL n/a Methods “Simulation”
PTA table reproduced below n/a Table 3
Sampling schedules used for the NCA replication n/a Supplementary Table 1

Structural verification (typical values)

These checks use rxode2::zeroRe() so both sides of each comparison use the same parameter values. The differences are pure numerical error, so the bounds are tight on purpose.

ui_typ <- rxode2::zeroRe(ui)

# Reference covariate values reproduce the typical parameter values exactly.
WT_REF <- 68.4
ALT_REF <- 16

# Per-study medians (Bae 2026 Table 1), used for the NCA replication so each
# published observed value is compared against its own study's covariates.
study_cov <- tibble::tribble(
  ~study,    ~WT,  ~ALT,
  "study 1", 67.6, 14,
  "study 2", 67.7, 20,
  "study 3", 70.6, 18
)

solve_typ <- function(events) {
  out <- rxode2::rxSolve(ui_typ, events = events,
                         keep = c("treatment", "WT", "ALT")) |>
    as.data.frame() |>
    dplyr::filter(!is.na(Cc))
  # rxSolve omits the `id` column entirely when the event table holds a single
  # subject; downstream PKNCA grouping needs it, so restore it.
  if (is.null(out$id)) {
    out$id <- if (is.null(events$id)) {
      1L
    } else {
      stopifnot(length(unique(events$id)) == 1L)
      events$id[[1]]
    }
  }
  out
}

Disposition eigenvalues and terminal half-life

The three disposition half-lives come from the eigenvalues of the three-compartment rate matrix built from the packaged parameter values.

th <- ui$theta
CL_typ <- exp(th[["lcl"]])
V1 <- exp(th[["lvc"]])
V2 <- exp(th[["lvp"]])
V3 <- exp(th[["lvp2"]])
Q12 <- exp(th[["lq"]])
Q13 <- exp(th[["lq2"]])

k10 <- CL_typ / V1
k12 <- Q12 / V1; k21 <- Q12 / V2
k13 <- Q13 / V1; k31 <- Q13 / V3
A <- matrix(c(-(k10 + k12 + k13), k21,  k31,
                            k12, -k21,    0,
                            k13,    0, -k31),
            nrow = 3, byrow = TRUE)
lambda <- sort(-Re(eigen(A)$values))
half_lives <- log(2) / lambda

Vss <- V1 + V2 + V3
tibble::tibble(
  phase = c("terminal", "intermediate", "rapid distribution"),
  `rate constant (1/h)` = lambda,
  `half-life (h)` = half_lives
) |>
  knitr::kable(digits = c(0, 5, 3),
               caption = "Disposition phases implied by the packaged parameters.")
Disposition phases implied by the packaged parameters.
phase rate constant (1/h) half-life (h)
terminal 0.02126 32.598
intermediate 1.33449 0.519
rapid distribution 24.46709 0.028 

cat(sprintf("Vss = V1 + V2 + V3 = %.1f L;  MRT = Vss/CL = %.1f h\n", Vss, Vss / CL_typ))
#> Vss = V1 + V2 + V3 = 98.8 L;  MRT = Vss/CL = 45.5 h

# Bae 2026 reports t1/2 of 27-45 h after single IV dosing and 27-40 h after SC
# dosing (Introduction, summarising studies 1, 2 and 3).
stopifnot(
  half_lives[1] > 27, half_lives[1] < 45,
  # The rapid distribution phase is sub-minute-scale, which is why the observed
  # Cmax after a short IV infusion is strongly sampling-limited (see below).
  half_lives[3] < 0.1
)

The terminal half-life of 32.6 h sits inside the 27-45 h range Bae 2026 reports across the three trials.

Bioavailability is applied to the SC depot only

f(depot) must scale SC exposure by exactly F and leave IV exposure untouched. The AUC ratio of a 40 mg SC dose to a 40 mg IV dose therefore has to equal F.

long_grid <- seq(0, 1000, by = 0.05)
ev_sc <- rxode2::et(amt = 40, cmt = "depot") |>
  rxode2::et(long_grid, cmt = "central") |>
  as.data.frame() |>
  dplyr::mutate(WT = WT_REF, ALT = ALT_REF, treatment = "40 mg SC")
ev_iv <- rxode2::et(amt = 40, dur = 0.5, cmt = "central") |>
  rxode2::et(long_grid, cmt = "central") |>
  as.data.frame() |>
  dplyr::mutate(WT = WT_REF, ALT = ALT_REF, treatment = "40 mg IV")

auc_trap <- function(d) sum(diff(d$time) * (head(d$Cc, -1) + tail(d$Cc, -1)) / 2)
f_ratio <- auc_trap(solve_typ(ev_sc)) / auc_trap(solve_typ(ev_iv))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalfdepot'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalfdepot'
F_model <- exp(th[["lfdepot"]])

cat(sprintf("AUC(SC) / AUC(IV) = %.5f;  packaged F = %.5f\n", f_ratio, F_model))
#> AUC(SC) / AUC(IV) = 0.81996;  packaged F = 0.82000
stopifnot(abs(f_ratio / F_model - 1) < 1e-3)

Steady-state Cavg closed form

For a linear model the average concentration over a complete dosing cycle is exactly (absorbed dose per cycle) / (CL x cycle length). Bae 2026’s simulated regimen is twice weekly (days 1 and 4), so a week is one complete cycle. This identity is what makes the Table 3 reproduction below closed-form, so it is worth checking against the solved ODE system.

DOSE_DAYS <- c(0, 72)             # days 1 and 4 of each week
WEEK <- 168
dose_times <- sort(as.vector(outer(DOSE_DAYS, WEEK * (0:12), "+")))

# Bae 2026 Methods "Simulation": virtual patients assumed 60 kg and ALT 40 IU/L.
WT_SIM <- 60
ALT_SIM <- 40

# Cavg depends only on CL and F: no weight covariate was retained on CL, so
# WT does not enter this closed form.
cavg_closed_form <- function(dose, ALT = ALT_SIM) {
  cl <- CL_typ * (ALT / 16)^th[["e_alt_cl"]]
  length(DOSE_DAYS) * dose * F_model * 1000 / (WEEK * cl)
}

ev_ss <- rxode2::et(amt = 30, cmt = "depot", time = dose_times) |>
  rxode2::et(seq(2 * WEEK, 3 * WEEK, by = 0.1), cmt = "central") |>
  as.data.frame() |>
  dplyr::mutate(WT = WT_SIM, ALT = ALT_SIM, treatment = "30 mg 2x/week")
sim_ss <- solve_typ(ev_ss) |> dplyr::filter(time >= 2 * WEEK)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalfdepot'

cavg_solved <- auc_trap(sim_ss) / WEEK
cat(sprintf("week 3 Cavg (30 mg twice weekly): solved %.3f, closed form %.3f ng/mL (%+.4f%%)\n",
            cavg_solved, cavg_closed_form(30),
            100 * (cavg_solved / cavg_closed_form(30) - 1)))
#> week 3 Cavg (30 mg twice weekly): solved 158.533, closed form 158.575 ng/mL (-0.0263%)
stopifnot(abs(cavg_solved / cavg_closed_form(30) - 1) < 2e-3)

Replicating the phase 1 designs and published NCA

Each arm below reproduces a published dose level on the sampling schedule that trial actually used (Supplementary Table 1) and with that study’s median covariates. Sampling fidelity matters here: the rapid distribution half-life is 1.7 min, so the observed Cmax after a short IV infusion depends heavily on when blood was drawn.

mn <- function(x) x / 60

grid_s1c1 <- c(0, mn(c(10, 20, 30, 32, 35, 40, 45)), 1, 1.5, 2, 3, 4, 6, 8, 12,
               24, 48, 72, 96, 120, 144)
grid_s1c2 <- c(0, mn(c(2, 5, 7, 10, 15, 20, 35)), 1, 1.5, 2, 3, 4, 6, 8, 12,
               24, 48, 72, 96, 120, 144)
grid_s3sc <- c(0, mn(30), 1, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 12, 24, 48, 96, 144)
grid_s3iv <- c(0, mn(c(10, 20, 30, 45)), 1, 1.5, 2, 3, 4, 6, 8, 12, 24, 48, 96, 144)
# Study 2: day 1 rich profile, days 2-6 pre-dose troughs, day 7 rich profile
# with a 144 h terminal tail. The day-2 pre-dose sample is the 24 h point of
# the first dosing interval.
grid_s2d <- c(0, mn(c(10, 20, 30, 32, 35, 40, 45)), 1, 1.5, 2, 3, 4, 6, 8, 12)
grid_s2 <- sort(unique(c(grid_s2d, 24, 48, 72, 96, 120,
                        144 + grid_s2d, 144 + c(24, 48, 72, 96, 120, 144))))

make_arm <- function(id, treatment, study, dose, cmt, obs, dur = NULL,
                     dose_time = 0) {
  ev <- if (is.null(dur)) {
    rxode2::et(amt = dose, cmt = cmt, time = dose_time)
  } else {
    rxode2::et(amt = dose, dur = dur, cmt = cmt, time = dose_time)
  }
  cv <- study_cov[study_cov$study == study, ]
  stopifnot(nrow(cv) == 1L)
  rxode2::et(ev, obs, cmt = "central") |>
    as.data.frame() |>
    dplyr::mutate(id = id, treatment = treatment,
                  WT = cv$WT, ALT = cv$ALT)
}

arms <- dplyr::bind_rows(
  make_arm(1L, "200 mg IV, 30-min (study 1 cohort 1)", "study 1",
           200, "central", grid_s1c1, dur = 0.5),
  make_arm(2L, "200 mg IV, 5-min (study 1 cohort 2)", "study 1",
           200, "central", grid_s1c2, dur = 5 / 60),
  make_arm(3L, "40 mg SC (study 3 cohort 1)", "study 3",
           40, "depot", grid_s3sc),
  make_arm(4L, "40 mg IV, 30-min (study 3 cohort 1)", "study 3",
           40, "central", grid_s3iv, dur = 0.5)
)
stopifnot(!anyDuplicated(arms[, c("id", "time", "evid")]))

arm_qd <- make_arm(5L, "30 mg IV QD x7 (study 2)", "study 2",
                   30, "central", grid_s2, dur = 0.5,
                   dose_time = 24 * (0:6))

sim_arms <- solve_typ(arms)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalfdepot'
#> Warning: multi-subject simulation without without 'omega'
sim_qd <- solve_typ(arm_qd)
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalvp2', 'etalq', 'etalq2', 'etalfdepot'
stopifnot(all(sim_arms$Cc >= 0), all(sim_qd$Cc >= 0), nrow(sim_arms) > 0)

PKNCA on the single-dose arms 

conc_single <- sim_arms |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

dose_single <- arms |>
  dplyr::filter(evid %in% c(1, 4)) |>
  dplyr::select(id, time, amt, treatment)

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

intervals_single <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, auclast = TRUE, aucinf.obs = TRUE, half.life = TRUE
)
res_single <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals_single)
)

PKNCA on the multiple-dose arm

AUC over the dosing interval on day 1 and on day 7 gives both the steady-state AUCtau that Bae 2026 used as a safety threshold and the accumulation ratio.

conc_qd <- sim_qd |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)
dose_qd <- arm_qd |>
  dplyr::filter(evid %in% c(1, 4)) |>
  dplyr::select(id, time, amt, treatment)

conc_obj_qd <- PKNCA::PKNCAconc(conc_qd, Cc ~ time | treatment + id,
                                concu = "ng/mL", timeu = "h")
dose_obj_qd <- PKNCA::PKNCAdose(dose_qd, amt ~ time | treatment + id,
                                doseu = "mg")

intervals_qd <- data.frame(
  start = c(0, 144), end = c(24, 168),
  auclast = TRUE, cmax = TRUE
)
res_qd <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_qd, dose_obj_qd, intervals = intervals_qd)
)

qd_tbl <- as.data.frame(res_qd$result)
get_qd <- function(param, start) {
  v <- qd_tbl$PPORRES[qd_tbl$PPTESTCD == param & qd_tbl$start == start]
  if (length(v) != 1L) stop("no unique ", param, " at start = ", start)
  v
}
auc_tau_d1 <- get_qd("auclast", 0)
auc_tau_d7 <- get_qd("auclast", 144)
accumulation <- auc_tau_d7 / auc_tau_d1

cat(sprintf("AUCtau day 1 = %.0f, day 7 = %.0f ng*h/mL; accumulation = %.2f\n",
            auc_tau_d1, auc_tau_d7, accumulation))
#> AUCtau day 1 = 5837, day 7 = 13935 ng*h/mL; accumulation = 2.39

Bae 2026 reports an accumulation ratio of 2.58-2.79 after multiple once-daily IV dosing (Introduction). The model gives 2.39, modestly below that range; the published figure comes from NCA on the observed sparse day-1 profile, whereas this one uses the model’s own prediction on the same nominal schedule.

Comparison against the published NCA values 

# Bae 2026 published observed values. Cmax after the 200 mg single dose and
# AUCtau after 30 mg once daily are the two safety anchors from Methods
# "Simulation"; Tmax ~2 h and t1/2 27-40 h (SC) / 27-45 h (IV) are from the
# Introduction. The 200 mg Cmax is reported for "Study 1" without naming a
# cohort, and study 1 gave 200 mg in BOTH the 30-min and the 5-min cohort, so
# it is carried on a pooled row (see "Assumptions and deviations").
published <- tibble::tribble(
  ~treatment,                                ~cmax,    ~tmax, ~half.life,
  "40 mg SC (study 3 cohort 1)",                NA,      2.0,       33.5,
  "40 mg IV, 30-min (study 3 cohort 1)",        NA,       NA,       36.0
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = res_single,
  reference = published,
  by = "treatment",
  units = c(cmax = "ng/mL", tmax = "h", half.life = "h"),
  tolerance_pct = 20
)
knitr::kable(
  cmp,
  caption = paste("Simulated (typical-value) vs. published NCA for MIT-001.",
                  "* marks a difference above 20%.",
                  "Published t1/2 references are the midpoints of the",
                  "reported 27-40 h (SC) and 27-45 h (IV) ranges."),
  digits = 3
)
Simulated (typical-value) vs. published NCA for MIT-001. * marks a difference above 20%. Published t1/2 references are the midpoints of the reported 27-40 h (SC) and 27-45 h (IV) ranges.
NCA parameter treatment Reference Simulated % diff
Cmax (ng/mL) 40 mg SC (study 3 cohort 1) — 345 —
Cmax (ng/mL) 40 mg IV, 30-min (study 3 cohort 1) — 1150 —
Tmax (h) 40 mg SC (study 3 cohort 1) 2 2 +0.0%
Tmax (h) 40 mg IV, 30-min (study 3 cohort 1) — 0.5 —
t½ (h) 40 mg SC (study 3 cohort 1) 33.5 33.9 +1.2%
t½ (h) 40 mg IV, 30-min (study 3 cohort 1) 36 34 -5.5%

The two absolute exposure anchors are checked directly, because each needs its own interval and pooling rule.

single_tbl <- as.data.frame(res_single$result)
get_single <- function(trt, param) {
  v <- single_tbl$PPORRES[single_tbl$treatment == trt &
                            single_tbl$PPTESTCD == param]
  if (length(v) != 1L) stop("no unique ", param, " for ", trt)
  v
}

cmax_30min <- get_single("200 mg IV, 30-min (study 1 cohort 1)", "cmax")
cmax_5min <- get_single("200 mg IV, 5-min (study 1 cohort 2)", "cmax")
cmax_pooled <- mean(c(cmax_30min, cmax_5min))

anchors <- tibble::tibble(
  quantity = c("Cmax, 200 mg IV, 30-min cohort (ng/mL)",
               "Cmax, 200 mg IV, 5-min cohort (ng/mL)",
               "Cmax, 200 mg IV, mean of the two cohorts (ng/mL)",
               "AUCtau day 7, 30 mg IV QD (ng*h/mL)"),
  simulated = c(cmax_30min, cmax_5min, cmax_pooled, auc_tau_d7),
  published = c(NA, NA, 10777.2, 16167.4)
) |>
  dplyr::mutate(`difference (%)` = 100 * (simulated / published - 1))

knitr::kable(anchors, digits = 1,
             caption = "Absolute exposure anchors from Bae 2026 Methods.")
Absolute exposure anchors from Bae 2026 Methods.
quantity simulated published difference (%)
Cmax, 200 mg IV, 30-min cohort (ng/mL) 5838.5 NA NA
Cmax, 200 mg IV, 5-min cohort (ng/mL) 14850.1 NA NA
Cmax, 200 mg IV, mean of the two cohorts (ng/mL) 10344.3 10777.2 -4.0
AUCtau day 7, 30 mg IV QD (ng*h/mL) 13935.5 16167.4 -13.8 

stopifnot(
  # The published 200 mg Cmax must lie between the two infusion-duration
  # cohorts, and the pooled mean must be within 20% of it.
  cmax_30min < 10777.2, cmax_5min > 10777.2,
  abs(cmax_pooled / 10777.2 - 1) < 0.20,
  # Steady-state AUCtau within 20% of the observed study 2 value.
  abs(auc_tau_d7 / 16167.4 - 1) < 0.20,
  # Tmax after SC dosing lands on the 2 h sample, as published.
  get_single("40 mg SC (study 3 cohort 1)", "tmax") == 2
)

Reproducing Table 3: probability of target attainment

Bae 2026 simulated 1000 virtual OM patients (60 kg, ALT 40 IU/L) at twice-weekly SC doses of 10-100 mg for 13 weeks and tabulated, for weeks 1, 3 and 13:

  • efficacy – the fraction with Cavg above the 87 ng/mL minimum effective concentration from the mouse OM model;
  • safety – the fraction with AUCtau below 16167.4 ng*h/mL and Cmax below 10777.2 ng/mL.

The safety metric needed reconstruction. The thresholds are per-24-h values taken from once-daily dosing, while the simulated regimen is twice weekly, so the comparable quantity is the regimen’s average daily AUC, Cavg * 24. Reconstructing it that way reproduces all 30 published safety cells (below), which no other reading of “AUCtau” does.

Because both criteria reduce to per-subject thresholds on Cavg, and Cavg at steady state has the closed form verified above, the weeks 3 and 13 rows of Table 3 can be reproduced analytically – no cohort, no seed, no Monte Carlo error. Only IIV on CL and on F enters, because Cavg depends on no other parameter.

om <- ui$omega
sd_log_cavg <- sqrt(om["etalcl", "etalcl"] + om["etalfdepot", "etalfdepot"])

C_EFF <- 87            # ng/mL, mouse minimum effective Cavg
AUCTAU_MAX <- 16167.4  # ng*h/mL, observed after 30 mg IV once daily
CMAX_MAX <- 10777.2    # ng/mL, observed after 200 mg single IV
cavg_safety_limit <- AUCTAU_MAX / 24

doses <- seq(10, 100, by = 10)
cavg_typ <- cavg_closed_form(doses)

pta_cf <- tibble::tibble(
  dose = doses,
  `Cavg (ng/mL)` = cavg_typ,
  efficacy = 100 * pnorm((log(cavg_typ) - log(C_EFF)) / sd_log_cavg),
  safety = 100 * pnorm((log(cavg_safety_limit) - log(cavg_typ)) / sd_log_cavg)
)

# Bae 2026 Table 3, weeks 3 and 13.
pub_tbl <- tibble::tribble(
  ~dose, ~eff_wk3, ~eff_wk13, ~saf_wk3, ~saf_wk13,
     10,      4.9,       5.1,    100.0,     100.0,
     20,     75.9,      75.8,    100.0,     100.0,
     30,     97.4,      98.0,    100.0,     100.0,
     40,    100.0,      99.9,    100.0,     100.0,
     50,     99.9,      99.9,    100.0,     100.0,
     60,    100.0,     100.0,     99.5,      99.3,
     70,    100.0,     100.0,     98.3,      97.7,
     80,    100.0,     100.0,     93.7,      93.6,
     90,    100.0,     100.0,     87.0,      86.6,
    100,    100.0,     100.0,     77.0,      76.1
)

pta_cmp <- pta_cf |>
  dplyr::left_join(pub_tbl, by = "dose") |>
  dplyr::mutate(
    d_eff_wk3 = efficacy - eff_wk3, d_eff_wk13 = efficacy - eff_wk13,
    d_saf_wk3 = safety - saf_wk3,   d_saf_wk13 = safety - saf_wk13
  )

pta_cmp |>
  dplyr::transmute(
    `Dose (mg)` = dose,
    `Cavg (ng/mL)` = `Cavg (ng/mL)`,
    `Efficacy PTA, model (%)` = efficacy,
    `Efficacy PTA, wk 3 (%)` = eff_wk3,
    `Efficacy PTA, wk 13 (%)` = eff_wk13,
    `Safety PTA, model (%)` = safety,
    `Safety PTA, wk 3 (%)` = saf_wk3,
    `Safety PTA, wk 13 (%)` = saf_wk13
  ) |>
  knitr::kable(digits = 1,
               caption = paste("Closed-form steady-state PTA vs. Bae 2026",
                               "Table 3 (weeks 3 and 13)."))
Closed-form steady-state PTA vs. Bae 2026 Table 3 (weeks 3 and 13).
Dose (mg) Cavg (ng/mL) Efficacy PTA, model (%) Efficacy PTA, wk 3 (%) Efficacy PTA, wk 13 (%) Safety PTA, model (%) Safety PTA, wk 3 (%) Safety PTA, wk 13 (%)
10 52.9 4.9 4.9 5.1 100.0 100.0 100.0
20 105.7 74.1 75.9 75.8 100.0 100.0 100.0
30 158.6 97.7 97.4 98.0 100.0 100.0 100.0
40 211.4 99.8 100.0 99.9 100.0 100.0 100.0
50 264.3 100.0 99.9 99.9 99.9 100.0 100.0
60 317.1 100.0 100.0 100.0 99.4 99.5 99.3
70 370.0 100.0 100.0 100.0 97.7 98.3 97.7
80 422.9 100.0 100.0 100.0 93.9 93.7 93.6
90 475.7 100.0 100.0 100.0 87.6 87.0 86.6
100 528.6 100.0 100.0 100.0 79.0 77.0 76.1 

all_diffs <- c(pta_cmp$d_eff_wk3, pta_cmp$d_eff_wk13,
               pta_cmp$d_saf_wk3, pta_cmp$d_saf_wk13)
stopifnot(length(all_diffs) == 40L, !anyNA(all_diffs))
cat(sprintf("40 published cells: max|diff| = %.2f pp, median|diff| = %.2f pp, mean = %+.2f pp\n",
            max(abs(all_diffs)), median(abs(all_diffs)), mean(all_diffs)))
#> 40 published cells: max|diff| = 2.90 pp, median|diff| = 0.04 pp, mean = +0.06 pp
stopifnot(
  max(abs(all_diffs)) < 3.5,
  median(abs(all_diffs)) < 1.0,
  abs(mean(all_diffs)) < 1.0
)

The closed form lands on all 40 published steady-state cells, which validates the structural model, the ALT covariate on CL, the SC bioavailability and the IIV magnitudes on CL and F simultaneously.

The paper’s three headline conclusions follow directly and are asserted here:

pta_at <- function(dose, what) pta_cf[[what]][pta_cf$dose == dose]
stopifnot(
  # "with a 30 mg SC dose, 97.4% of patients were anticipated to reach the
  #  target concentration" and 20 mg does not clear 90%.
  pta_at(30, "efficacy") > 90, pta_at(20, "efficacy") < 90,
  # "up to 50 mg twice weekly SC dosing would not reach the safety thresholds"
  pta_at(50, "safety") > 99.5,
  # "Under the 80 mg twice weekly dosing regimen, exposures remained within the
  #  predefined safety threshold" but 90 mg puts >10% of patients over it.
  pta_at(80, "safety") > 90, pta_at(90, "safety") < 90
)
cat("Recommended range 30-80 mg twice weekly SC is reproduced.\n")
#> Recommended range 30-80 mg twice weekly SC is reproduced.

Monte Carlo cross-check, including week 1

Week 1 is not at steady state, so it needs simulation. This cohort also confirms that the closed form is not hiding a discretisation or accumulation error.

set.seed(20260831)
rxode2::rxSetSeed(20260831)

N_PER_ARM <- 200        # per-arm cap; ample for a PTA cross-check
mc_doses <- c(20, 30, 80, 100)
week_starts <- c(`1` = 0, `3` = 2 * WEEK, `13` = 12 * WEEK)
obs_grid <- sort(unique(as.vector(
  vapply(week_starts, function(t0) seq(t0, t0 + WEEK, by = 1),
         numeric(WEEK + 1))
)))

# One rxSolve call per arm: rxSolve on an rxUi scales poorly with the subject
# count within a single call. With `nSub =` the subject key comes back as
# `sim.id` rather than `id`, so normalise it here.
solve_mc <- function(dose) {
  ev <- rxode2::et(amt = dose, cmt = "depot", time = dose_times) |>
    rxode2::et(obs_grid, cmt = "central") |>
    as.data.frame() |>
    dplyr::mutate(WT = WT_SIM, ALT = ALT_SIM, treatment = paste(dose, "mg"))
  s <- rxode2::rxSolve(ui, events = ev, nSub = N_PER_ARM,
                       keep = c("treatment")) |>
    as.data.frame() |>
    dplyr::filter(!is.na(Cc))
  if (is.null(s$id) && !is.null(s$sim.id)) s$id <- s$sim.id
  # Guard: a cohort that silently collapsed to one subject would make every
  # PTA cell 0% or 100% and the gates below would still "pass".
  stopifnot(!is.null(s$id), dplyr::n_distinct(s$id) == N_PER_ARM)
  s
}

pta_mc <- lapply(mc_doses, function(d) {
  s <- solve_mc(d)
  lapply(names(week_starts), function(wk) {
    t0 <- week_starts[[wk]]
    w <- s |> dplyr::filter(time >= t0, time <= t0 + WEEK)
    per_sub <- w |>
      dplyr::group_by(id) |>
      dplyr::summarise(
        Cavg = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2) / WEEK,
        Cmax = max(Cc), .groups = "drop"
      )
    tibble::tibble(
      dose = d, week = wk,
      efficacy = 100 * mean(per_sub$Cavg > C_EFF),
      safety = 100 * mean(per_sub$Cavg * 24 < AUCTAU_MAX &
                            per_sub$Cmax < CMAX_MAX),
      max_cmax = max(per_sub$Cmax)
    )
  }) |> dplyr::bind_rows()
}) |> dplyr::bind_rows()

pub_long <- pub_tbl |>
  dplyr::select(dose, eff_wk3, eff_wk13, saf_wk3, saf_wk13) |>
  tidyr::pivot_longer(-dose, names_to = "k", values_to = "published") |>
  dplyr::mutate(
    criterion = ifelse(startsWith(k, "eff"), "efficacy", "safety"),
    week = sub("^.*_wk", "", k)
  ) |>
  dplyr::select(-k)

pub_wk1 <- tibble::tribble(
  ~dose, ~criterion, ~published,
     20,  "efficacy",      57.6,
     30,  "efficacy",      95.4,
     80,  "efficacy",     100.0,
    100,  "efficacy",     100.0,
     20,    "safety",     100.0,
     30,    "safety",     100.0,
     80,    "safety",      99.6,
    100,    "safety",      96.7
) |> dplyr::mutate(week = "1")

mc_cmp <- pta_mc |>
  dplyr::select(-max_cmax) |>
  tidyr::pivot_longer(c(efficacy, safety), names_to = "criterion",
                      values_to = "simulated") |>
  dplyr::inner_join(dplyr::bind_rows(pub_long, pub_wk1),
                    by = c("dose", "week", "criterion")) |>
  dplyr::mutate(`difference (pp)` = simulated - published) |>
  dplyr::arrange(criterion, dose, as.numeric(week))

mc_cmp |>
  dplyr::rename("Dose (mg)" = dose, "Week" = week, "Criterion" = criterion,
                "Simulated PTA (%)" = simulated, "Published PTA (%)" = published) |>
  knitr::kable(digits = 1,
               caption = paste("Monte Carlo PTA (", N_PER_ARM,
                               "subjects per arm) vs. Bae 2026 Table 3.",
                               "Week 1 is pre-steady-state."))
Monte Carlo PTA ( 200 subjects per arm) vs. Bae 2026 Table 3. Week 1 is pre-steady-state.
Dose (mg) Week Criterion Simulated PTA (%) Published PTA (%) difference (pp)
20 1 efficacy 68.5 57.6 10.9
20 3 efficacy 77.5 75.9 1.6
20 13 efficacy 77.5 75.8 1.7
30 1 efficacy 97.5 95.4 2.1
30 3 efficacy 98.5 97.4 1.1
30 13 efficacy 98.5 98.0 0.5
80 1 efficacy 100.0 100.0 0.0
80 3 efficacy 100.0 100.0 0.0
80 13 efficacy 100.0 100.0 0.0
100 1 efficacy 100.0 100.0 0.0
100 3 efficacy 100.0 100.0 0.0
100 13 efficacy 100.0 100.0 0.0
20 1 safety 100.0 100.0 0.0
20 3 safety 100.0 100.0 0.0
20 13 safety 100.0 100.0 0.0
30 1 safety 100.0 100.0 0.0
30 3 safety 100.0 100.0 0.0
30 13 safety 100.0 100.0 0.0
80 1 safety 99.0 99.6 -0.6
80 3 safety 95.0 93.7 1.3
80 13 safety 95.0 93.6 1.4
100 1 safety 94.5 96.7 -2.2
100 3 safety 80.0 77.0 3.0
100 13 safety 80.0 76.1 3.9 

# Cohort proportions near a steep threshold are sampling-noisy, so this gate is
# on the aggregate rather than on any single cell (see CLAUDE.md on vignette
# assertions over a simulated cohort).
cat(sprintf("Monte Carlo vs published: median|diff| = %.1f pp, mean = %+.1f pp\n",
            median(abs(mc_cmp$`difference (pp)`)), mean(mc_cmp$`difference (pp)`)))
#> Monte Carlo vs published: median|diff| = 0.2 pp, mean = +1.0 pp
cat(sprintf("highest single-subject Cmax across all arms and weeks: %.0f ng/mL (threshold %.0f)\n",
            max(pta_mc$max_cmax), CMAX_MAX))
#> highest single-subject Cmax across all arms and weeks: 1844 ng/mL (threshold 10777)
stopifnot(
  nrow(mc_cmp) == 24L,
  median(abs(mc_cmp$`difference (pp)`)) < 5,
  abs(mean(mc_cmp$`difference (pp)`)) < 5,
  # The Cmax limb of the safety criterion never binds: even the most extreme
  # subject at the highest dose stays far below the threshold, so safety
  # attainment is governed entirely by AUCtau. This is what lets the closed
  # form above reproduce the published safety column.
  max(pta_mc$max_cmax) < 0.5 * CMAX_MAX
)

Replicating Figure 4

Figure 4 of Bae 2026 shows median and 5th-95th percentile concentration-time profiles for weeks 1, 3 and 13, with the 87 ng/mL mouse effective Cavg marked.

fig_doses <- c(20, 30, 80)
fig_sim <- lapply(fig_doses, solve_mc) |> dplyr::bind_rows()

fig_sim |>
  dplyr::mutate(
    week = dplyr::case_when(
      time <= WEEK ~ "Week 1",
      time >= 2 * WEEK & time <= 3 * WEEK ~ "Week 3",
      time >= 12 * WEEK ~ "Week 13"
    ),
    t_in_week = time - dplyr::case_when(
      time <= WEEK ~ 0,
      time >= 2 * WEEK & time <= 3 * WEEK ~ 2 * WEEK,
      TRUE ~ 12 * WEEK
    )
  ) |>
  dplyr::filter(!is.na(week)) |>
  dplyr::mutate(week = factor(week, levels = c("Week 1", "Week 3", "Week 13"))) |>
  dplyr::group_by(week, treatment, t_in_week) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05), Q50 = median(Cc), Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  ggplot(aes(t_in_week, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue") +
  geom_hline(yintercept = C_EFF, linetype = "dashed", colour = "red") +
  facet_grid(treatment ~ week) +
  scale_y_log10() +
  labs(x = "Time within week (h)", y = "MIT-001 concentration (ng/mL)",
       title = "Simulated twice-weekly SC profiles",
       caption = paste("Replicates Figure 4 of Bae 2026. Dashed red line is the",
                       "87 ng/mL mouse effective Cavg.")) +
  theme_bw()
#> Warning in scale_y_log10(): log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.

Assumptions and deviations

  • Residual error scale. Table 2 reports “Proportional residual error 0.108” without units, which NONMEM users write both as a variance and as an SD. It is read here as a proportional SD (10.8% CV), on the evidence that the whole variability block of Table 2 is reported on the SD / CV scale: the IIV rows are given as percent CV, and their RSEs are on that same scale. The RSE column proves it – an RSE of 7.8% on the CL IIV (and 9.5% on the V1 IIV) is below sqrt(2 / (119 - 1)) = 13.0%, the theoretical floor for the RSE of a variance estimated from 119 subjects, so those RSEs cannot be variance-scale. Reading 0.108 as a variance instead would give a 32.9% proportional error, which also sits well above the assay’s validated precision (accuracy within 15%, CV within 20%). No validation gate in this vignette is sensitive to the choice, because Cmax never approaches its threshold and AUC-based metrics average residual error out.
  • IIV conversion. Percent CVs from Table 2 were converted with omega^2 = log(CV^2 + 1). The naive omega = CV reading changes the reproduced PTA by well under a percentage point, so the published table does not discriminate between them.
  • V2 weight exponent. Table 2 footnote b prints V2 = theta3 * (WT / 68.4) with no exponent, and Table 2 lists a weight exponent (theta5) only for V3. The V2 exponent is therefore encoded as fixed(1) – linear in weight, not estimated.
  • Safety metric reconstruction. Bae 2026 does not state how the once-daily AUCtau threshold was applied to a twice-weekly regimen. It is reconstructed here as the regimen’s average daily AUC, Cavg * 24. This reproduces all 30 published safety cells to within 2.9 percentage points; treating “AUCtau” as the AUC over the actual 72/96 h dosing interval, or as the AUC over the 24 h following a dose, both fail badly (0% and ~15% predicted attainment at 100 mg against a published 77%).
  • Per-study participant counts. Supplementary Table 1 gives 63 participants for study 1 and 26 for study 2; the main-text Table 1 “Male, n (%)” row reports these two the other way round. Supplementary Table 1 is used here, on three independent checks: (i) the published sample counts (1185 and 1038) divided by those participant counts give 19 and 40 samples per participant against scheduled schedules of ~22 and ~44, whereas the Table 1 ordering implies 46 samples per participant in study 1 from a 22-point schedule, which is impossible; (ii) study 1 ran 14 dose levels across two cohorts and study 2 ran 4, consistent with 63 and 26; (iii) the pooled ALT reference of 16 IU/L in Table 2 footnote a is recovered only when study 1 (ALT median 14 IU/L) is the larger study – the Table 1 ordering implies a pooled ALT near 18 IU/L. The total of 119 participants is the same either way, so no model parameter is affected.
  • Sign of the bootstrap CI on the ALT exponent. Table 2 prints the bootstrap 95% CI for theta2 as -0.178 (-0.28-0.079). The upper bound is almost certainly -0.079 with a dropped minus sign: the covariate survived backward elimination at the more stringent dOFV > 6.63 criterion, so its CI must exclude zero. Only the point estimate is used in the model, so nothing downstream depends on this.
  • 200 mg Cmax pooling. The 10777.2 ng/mL safety threshold is attributed to “Study 1” without naming a cohort, but study 1 gave 200 mg in both its 30-min and its 5-min infusion cohort. The model predicts 5839 and 14850 ng/mL for those two, bracketing the published value; their mean is within 4% of it. The gate above asserts the bracketing and the pooled agreement rather than matching a single cohort.
  • Bioavailability. The model estimate of 82% is used, not the 94% AUC-ratio value observed in the 6-participant SC/IV crossover; Bae 2026 reports the former as the model parameter and calls the two “comparable”.
  • Steady-state AUCtau. The model’s day-7 AUCtau for 30 mg once-daily IV is about 14% below the observed 16167.4 ng*h/mL. The model pools three studies and two routes, whereas the observed value comes from study 2 alone (26 of 119 participants), so a difference of this size is expected; it is inside the 20% tolerance and no parameter was adjusted.
  • Virtual cohort. Simulations use the paper’s own assumption of 60 kg and ALT 40 IU/L for OM patients (Methods “Simulation”), and per-study median weight and ALT for the phase 1 NCA replication. Cohorts are capped at 200 subjects per arm rather than the paper’s 1000, so Monte Carlo cells near a steep threshold carry a few percentage points of sampling noise; the tight gates in this vignette are therefore placed on the closed-form reproduction, which has no sampling error.
  • Supplement contents. Supplementary Table 1 (study designs, sampling schedules, sample counts) and Supplementary Figure 1 (pcVPC for the first 12 h after IV dosing) were retrieved and used. The supplement contains no parameter values and no NONMEM control stream; all parameters come from Table 2. No erratum or correction is listed for this article in Europe PMC.