Skip to contents

Emapalumab in macrophage activation syndrome (Brossard 2025)

Source: vignettes/articles/Brossard_2025_emapalumab.Rmd
Brossard_2025_emapalumab.Rmd

Model and source 

mod <- readModelDb("Brossard_2025_emapalumab")()
  • Citation: Brossard P. (2025) Emapalumab in Patients With Macrophage Activation Syndrome Associated With Still’s Disease: A Population Pharmacokinetic/Pharmacodynamic Analysis. Clin Transl Sci 18(2):e70163. doi:10.1111/cts.70163. The structural PK model reproduced in Table S2 is the final model of the upstream analysis Brossard P, Laveille C. (2024) Population Pharmacokinetics of the Anti-Interferon-Gamma Monoclonal Antibody Emapalumab: An Updated Analysis. Rheumatol Ther 11:869-880, doi:10.1007/s40744-024-00669-y, whose Results text supplies the four-point clearance answer key used to recover the unprinted age and total-bilirubin reference values.
  • Description: Two-compartment population PK model for emapalumab with allometric, age and total-bilirubin covariates and an interferon-gamma-dependent clearance component, linked to three parallel indirect-response (turnover) models for CXCL9, soluble IL-2 receptor alpha and ferritin in patients with macrophage activation syndrome associated with Still’s disease
  • Article: https://doi.org/10.1111/cts.70163
  • Supplement (Tables S1 / S2, Figure S1): https://doi.org/10.1111/cts.70163
  • Upstream structural PK analysis: Brossard P, Laveille C. Rheumatol Ther 2024;11:869-880, https://doi.org/10.1007/s40744-024-00669-y

Emapalumab is a fully human monoclonal antibody against interferon-gamma (IFN-gamma). Brossard 2025 re-estimated a two-compartment population PK model, originally built on a pooled primary-haemophagocytic-lymphohistiocytosis (HLH) and macrophage-activation-syndrome (MAS) dataset, on the 14 patients with MAS associated with Still’s disease from NCT03311854, and linked it to three parallel indirect-response (turnover) models for the hyperinflammation biomarkers CXCL9, soluble interleukin-2 receptor alpha (sIL-2Ralpha) and ferritin (Figure 1 of the source).

The paper’s central PK finding is that in MAS the target-mediated (“non-linear”) clearance component is switched off entirely: total IFN-gamma in MAS stays below the roughly 10,000 pg/mL threshold at which IFN-gamma begins to drive emapalumab disposition, so clearance is effectively linear and the half-life is close to that seen in healthy subjects.

Population

The PK/PD layer was estimated in the 14 patients with MAS associated with Still’s disease enrolled in NCT03311854 (Table S1 of the source). They were older and heavier than the primary-HLH patients who contributed the bulk of the structural PK dataset: median age 11.5 years (range 2.1-25.4 years), median body weight 45.5 kg (range 12.0-68.8 kg), 71.4% female, 78.6% White / 14.3% Asian / 7.1% unknown. Median total bilirubin was 12.8 umol/L (range 4.6-76.1) and median creatinine clearance 140.7 mL/min. All received 6 mg/kg IV loading dose, then 3 mg/kg every 3 days until Day 15 and twice weekly until Day 28.

The structural PK model itself was developed on a pooled n = 58 dataset (44 primary-HLH patients from NCT01818492 plus these 14 MAS patients, with long-term follow-up in NCT02069899).

The same information is available programmatically via readModelDb("Brossard_2025_emapalumab")()$population.

Source trace

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

Equation / parameter Value Source location
lcl (CLL) 0.0143 L/h per 70 kg Table 1 “CLL, Typical value”; Table S2 “CLL, L/h/70 kg”
lclnl (CLNL) 0.121 L/h per 70 kg Table 1 “CLNL, Intercept for IFNg = 1e6”; Table S2 “CLNL, L/h/70 kg”
lvc (V1) 3.08 L per 70 kg Table 1 “V1, Typical value”; Table S2 “V1, L/70 kg”
lq (Q) 0.105 L/h per 70 kg Table 1 “Q, Typical value”; Table S2 “Q, L/h/70 kg”
lvp (V2) 4.28 L per 70 kg Table 1 “V2, Typical value”; Table S2 “V2, L/70 kg”
e_wt_cl 0.75 (fixed) Table S2 CLs_BW
e_wt_vc 1 (fixed) Table S2 Vs_BW
e_age_cl 0.188 Table 1 “CL, Q age effect, exponent” (Table S2 CLs_AGE = 0.187)
e_age_vc -0.104 Table 1 “V1, Age effect, exponent”; Table S2 V1_AGE
e_tbili 0.162 Table 1 “CL, V1, V2 bilirubin effect, exponent”; Table S2 CLs_V1_TBIL
e_ifng_clnl 0.542 Table 1 “IFNg effect, exponent”; Table S2 CLNL_IFNg
e_mas_clnl -1 (fixed) Table S2 CLNL_MAS
etalcl 0.361 (SD) Table 1 “CL, Inter-individual variability”
etalvc, etalvp, correlation 0.207, 0.711, 0.569 Table 1 V1 / V2 IIV and “Correlation (V1, V2)”
propSd 0.301 Table 1 “RUV, Additive component, log”
lrbase_cxcl9 8400 ng/L Table 2 CXCL9 Baseline
lkdeg_cxcl9 0.414 /day Table 2 CXCL9 Degradation rate
logitimax_cxcl9 98.3% Table 2 CXCL9 Imax
propSd_cxcl9 46.4% Table 2 CXCL9 RUV
lrbase_sil2ra 6630 ng/L Table 2 sIL-2Ralpha Baseline (see Errata)
lkdeg_sil2ra 0.112 /day Table 2 sIL-2Ralpha Degradation rate
logitimax_sil2ra 87.3% Table 2 sIL-2Ralpha Imax
propSd_sil2ra 29.8% Table 2 sIL-2Ralpha RUV
lrbase_ferritin 15,300 ug/L Table 2 Ferritin Baseline
lkdeg_ferritin 0.207 /day Table 2 Ferritin Degradation rate
logitimax_ferritin 99.6% Table 2 Ferritin Imax
propSd_ferritin 38.9% Table 2 Ferritin RUV
PD IIV block (baselines, degradation rates, Imax) 1.28 / 0.694 / 1.07 / 0.696 / 0.746 / 0.891 / 1.5 / 0.37 / 1.65 Table 2 “Interindividual variability” rows
PD correlations 0.675, 0.805, 0.684 Table 2 “Correlation with baseline ferritin” / “Correlation with Imax”
Two-compartment disposition n/a Methods 2.1; Figure 1
Three parallel turnover models n/a Figure 1; Results 3.2 (“A parallel model … resulted in a better fit”)
Weight reference 70 kg 70 kg Table S2 (“per 70 kg” on every structural typical value)
Upstream clearance answer key 0.00218 / 0.00308 / 0.00623 / 0.01718 L/h Brossard & Laveille 2024 Results text (verified verbatim in the source PDF)
Upstream half-life answer key 19.2 / 13.8 / 7.18 / 3.12 days Brossard & Laveille 2024 Results text, same sentence
Covariate targets (age on CL/Q/V1; TBILI on CLL/V1/V2) n/a Brossard & Laveille 2024 Supplementary Table 1, runs 012-019
Age reference 25 y 25 y Derived - back-solved from the upstream four-point clearance answer key (see Errata)
Total-bilirubin reference 12.8 umol/L 12.8 umol/L Derived - back-solve returns a bilirubin factor of 1.0004; identified as the pooled median of Table S1 (see Errata)
lec50 (potency) 0.0248 ug/mL (fixed) NOT FROM THIS PAPER - Jacqmin 2022 Supplementary Table 7 (see Errata)

Structural falsifier: the upstream four-point clearance answer key

Brossard 2025 prints typical values but never states the age or total-bilirubin reference at which they apply. The upstream Brossard & Laveille 2024 Results text does supply an answer key: for “a 1-year-old patient weighing 10 kg with primary HLH”, emapalumab CL is 0.00218, 0.00308, 0.00623 and 0.01718 L/h at total serum IFN-gamma of 1e3, 1e4, 1e5 and 1e6 pg/mL. Because that reference patient carries the same structural parameters as Table 1 / Table S2, the four values pin the two unprinted reference values exactly.

modZero <- rxode2::zeroRe(mod)

ifng_levels <- c(1e3, 1e4, 1e5, 1e6)

key_cov <- data.frame(
  id      = seq_along(ifng_levels),
  WT      = 10,
  AGE     = 1,
  TBILI   = 12.8,
  IFNG    = ifng_levels,
  DIS_MAS = 0     # the answer-key patient has primary HLH, so CLNL is active
)

# The model declares four endpoints (Cc, cxcl9, sil2ra, ferritin), so rxode2
# needs `dvid` on every observation row to resolve the endpoint mapping; `cmt`
# alone is not sufficient. One dvid is enough -- all four observables plus every
# derived variable come back as columns on the solved rows.
key_ev <- key_cov |>
  dplyr::mutate(time = 0, amt = NA_real_, evid = 0L, cmt = "central", dvid = 1L)

key_sim <- rxode2::rxSolve(
  modZero, events = key_ev, omega = NA, sigma = NA,
  keep = c("IFNG", "DIS_MAS"), returnType = "data.frame"
)

published_cl <- c(0.00218, 0.00308, 0.00623, 0.01718)

key_tab <- key_sim |>
  dplyr::transmute(
    `Total IFN-gamma (pg/mL)` = format(IFNG, scientific = TRUE),
    `Simulated CL (L/h)`      = cl,
    `Published CL (L/h)`      = published_cl,
    `Difference (%)`          = 100 * (cl - published_cl) / published_cl
  )

knitr::kable(
  key_tab, digits = c(0, 5, 5, 2),
  caption = paste(
    "Simulated total clearance for the upstream reference patient",
    "(1 year old, 10 kg, primary HLH, total bilirubin 12.8 umol/L) against the",
    "four values in the Brossard & Laveille 2024 Results text."
  )
)
Simulated total clearance for the upstream reference patient (1 year old, 10 kg, primary HLH, total bilirubin 12.8 umol/L) against the four values in the Brossard & Laveille 2024 Results text.
Total IFN-gamma (pg/mL) Simulated CL (L/h) Published CL (L/h) Difference (%)
1e+03 0.00218 0.00218 -0.12
1e+04 0.00308 0.00308 -0.02
1e+05 0.00622 0.00623 -0.14
1e+06 0.01717 0.01718 -0.08 
# Strict check: every one of the four published clearances is reproduced to
# better than 1%. This is what identifies the 25-year age reference; the
# competing "centre at the dataset median age (1.7 y)" reading is off by 83%.
stopifnot(all(abs(key_sim$cl - published_cl) / published_cl < 0.01))

All four published clearances are reproduced to within 0.3%, which is inside the rounding of the published three-significant-figure values.

# The paper's central PK finding, as an exact identity: with DIS_MAS = 1 the
# IFN-gamma-dependent clearance component is identically zero at every
# IFN-gamma level, so clearance in MAS is purely linear.
mas_ev <- key_cov |>
  dplyr::mutate(DIS_MAS = 1, time = 0, amt = NA_real_, evid = 0L,
                cmt = "central", dvid = 1L)

mas_sim <- rxode2::rxSolve(
  modZero, events = mas_ev, omega = NA, sigma = NA,
  returnType = "data.frame"
)

stopifnot(all(mas_sim$clnl == 0))
stopifnot(isTRUE(all.equal(mas_sim$cl, mas_sim$cll)))

# ... and CL is then invariant to IFN-gamma across four orders of magnitude.
stopifnot(diff(range(mas_sim$cl)) == 0)

Virtual cohort

Original participant-level data are not public. The cohort below reproduces the MAS column of Table S1: 200 virtual patients whose age, body weight and total bilirubin match the published medians and ranges. Age and weight are rank-coupled so that older patients are heavier, which the published medians (11.5 years / 45.5 kg) imply but Table S1 does not tabulate jointly.

nSub <- 200

rtrunc <- function(n, rfun, lower, upper, ...) {
  out <- rfun(n, ...)
  while (any(bad <- out < lower | out > upper)) {
    out[bad] <- rfun(sum(bad), ...)
  }
  out
}

# Age: Table S1 MAS mean 10.4 y (CV 64.1%), range 2.1-25.4
age <- rtrunc(nSub, rnorm, 2.1, 25.4, mean = 10.4, sd = 0.641 * 10.4)

# Weight: Table S1 MAS mean 38.9 kg (CV 51.7%), range 12.0-68.8, rank-coupled to age
wt_marginal <- sort(rtrunc(nSub, rnorm, 12.0, 68.8, mean = 38.9, sd = 0.517 * 38.9))
wt <- wt_marginal[rank(age, ties.method = "first")]
wt <- pmin(pmax(wt * exp(rnorm(nSub, 0, 0.08)), 12.0), 68.8)

# Total bilirubin: Table S1 MAS median 12.8 umol/L, mean 22.4, range 4.6-76.1.
# Log-normal with meanlog = log(12.8) and sdlog chosen so exp(mu + s^2/2) = 22.4.
sdlog_tbili <- sqrt(2 * (log(22.4) - log(12.8)))
tbili <- rtrunc(nSub, rlnorm, 4.6, 76.1, meanlog = log(12.8), sdlog = sdlog_tbili)

cohort <- data.frame(
  id      = seq_len(nSub),
  WT      = wt,
  AGE     = age,
  TBILI   = tbili,
  IFNG    = 500,   # inert in MAS: DIS_MAS = 1 zeroes CLNL regardless of this value
  DIS_MAS = 1
)

knitr::kable(
  data.frame(
    Covariate  = c("Age (years)", "Body weight (kg)", "Total bilirubin (umol/L)"),
    `Simulated median` = c(median(cohort$AGE), median(cohort$WT), median(cohort$TBILI)),
    `Published median` = c(11.5, 45.5, 12.8),
    `Simulated range`  = c(
      sprintf("%.1f-%.1f", min(cohort$AGE), max(cohort$AGE)),
      sprintf("%.1f-%.1f", min(cohort$WT), max(cohort$WT)),
      sprintf("%.1f-%.1f", min(cohort$TBILI), max(cohort$TBILI))
    ),
    `Published range`  = c("2.1-25.4", "12.0-68.8", "4.6-76.1"),
    check.names = FALSE
  ),
  digits = 1,
  caption = "Virtual cohort covariates against Table S1 (MAS associated with sJIA column)."
)
Virtual cohort covariates against Table S1 (MAS associated with sJIA column).
Covariate Simulated median Published median Simulated range Published range
Age (years) 11.2 11.5 2.3-25.2 2.1-25.4
Body weight (kg) 39.9 45.5 12.1-68.8 12.0-68.8
Total bilirubin (umol/L) 14.4 12.8 4.8-75.7 4.6-76.1

Dosing regimen and simulation

Patients with MAS received emapalumab 6 mg/kg as a 1-hour intravenous infusion, then 3 mg/kg every 3 days until Day 15 and twice weekly until Day 28. The exact twice-weekly calendar is not printed; Days 18, 22, 25 and 28 are used here. The model runs on hours.

dose_days <- c(0, 3, 6, 9, 12, 15, 18, 22, 25, 28)
dose_mgkg <- c(6, rep(3, length(dose_days) - 1L))

doses <- cohort |>
  tidyr::crossing(data.frame(dose_day = dose_days, mgkg = dose_mgkg)) |>
  dplyr::mutate(
    time = dose_day * 24,
    amt  = mgkg * WT,
    rate = amt / 1,     # 1-hour infusion
    evid = 1L,
    cmt  = "central",
    dvid = NA_integer_
  ) |>
  dplyr::select(id, WT, AGE, TBILI, IFNG, DIS_MAS, time, amt, rate, evid, cmt, dvid)

obs_hours <- sort(unique(c(
  c(0, 1, 2, 4, 8, 12),                  # dense over the first infusion
  seq(24, 28 * 24, by = 12),             # every 12 h through the treatment period
  seq(28 * 24 + 24, 140 * 24, by = 24)   # daily through the washout
)))

obs <- cohort |>
  tidyr::crossing(data.frame(time = obs_hours)) |>
  dplyr::mutate(amt = NA_real_, rate = NA_real_, evid = 0L, cmt = "central",
                dvid = 1L) |>
  dplyr::select(id, WT, AGE, TBILI, IFNG, DIS_MAS, time, amt, rate, evid, cmt, dvid)

events <- dplyr::bind_rows(doses, obs) |>
  dplyr::arrange(id, time, dplyr::desc(evid)) |>
  as.data.frame()

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("WT", "AGE", "TBILI", "DIS_MAS"),
  returnType = "data.frame"
)

# rxSolve output has no evid column; the dosing rows come back as duplicated
# times, so keep the observation grid explicitly.
sim <- sim |> dplyr::filter(time %in% obs_hours)

stopifnot(dplyr::n_distinct(sim$id) == nSub)   # rxSolve can silently drop subjects

Replicating Figure 3: VPC of the three PD markers

Figure 3 of the source shows visual predictive checks for CXCL9, ferritin and sIL-2Ralpha over 56 days. The model-simulated median and 5th-95th percentile band are reproduced below, with proportional residual error applied at the per-marker RUV values from Table 2.

ruv <- c(cxcl9 = 0.464, sil2ra = 0.298, ferritin = 0.389)

pd_long <- sim |>
  dplyr::filter(time <= 56 * 24) |>
  dplyr::select(id, time, cxcl9, sil2ra, ferritin) |>
  tidyr::pivot_longer(c(cxcl9, sil2ra, ferritin),
                      names_to = "marker", values_to = "ipred") |>
  dplyr::mutate(
    day = time / 24,
    dv  = ipred * (1 + rnorm(dplyr::n(), 0, ruv[marker]))
  ) |>
  dplyr::filter(dv > 0)

pd_pct <- pd_long |>
  dplyr::group_by(marker, day) |>
  dplyr::summarise(
    lo  = quantile(dv, 0.05),
    med = median(dv),
    hi  = quantile(dv, 0.95),
    .groups = "drop"
  ) |>
  dplyr::mutate(marker = factor(
    marker,
    levels = c("cxcl9", "ferritin", "sil2ra"),
    labels = c("CXCL9 (ng/L)", "Ferritin (ug/L)", "sIL-2Ralpha (ng/L)")
  ))
ggplot(pd_pct, aes(day, med)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey70", alpha = 0.5) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~marker, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 56, by = 14)) +
  scale_y_log10() +
  labs(x = "Time after first dose (days)", y = "Absolute concentration") +
  theme_bw()
Replicates Figure 3 of Brossard 2025: simulated median (line) and 5th-95th percentile band (shaded) for each PD marker over 56 days.

Replicates Figure 3 of Brossard 2025: simulated median (line) and 5th-95th percentile band (shaded) for each PD marker over 56 days.

Structural check: the suppressed plateau equals baseline x (1 - Imax)

Because emapalumab concentrations in this regimen sit far above the potency term, inhibition is essentially saturated and each marker decays from its baseline towards rbase * (1 - Imax). That plateau is an exact structural identity of the turnover model and is independently readable off Figure 3.

typ_ev <- data.frame(
  id = 1L, WT = 45.5, AGE = 11.5, TBILI = 12.8, IFNG = 500, DIS_MAS = 1
) |>
  tidyr::crossing(data.frame(dose_day = dose_days, mgkg = dose_mgkg)) |>
  dplyr::mutate(time = dose_day * 24, amt = mgkg * 45.5, rate = amt / 1,
                evid = 1L, cmt = "central", dvid = NA_integer_) |>
  dplyr::select(id, WT, AGE, TBILI, IFNG, DIS_MAS, time, amt, rate, evid, cmt, dvid) |>
  dplyr::bind_rows(
    data.frame(id = 1L, WT = 45.5, AGE = 11.5, TBILI = 12.8, IFNG = 500,
               DIS_MAS = 1, time = obs_hours, amt = NA_real_, rate = NA_real_,
               evid = 0L, cmt = "central", dvid = 1L)
  ) |>
  dplyr::arrange(time, dplyr::desc(evid)) |>
  as.data.frame()

typ <- rxode2::rxSolve(
  modZero, events = typ_ev, omega = NA, sigma = NA,
  returnType = "data.frame"
)

plateau_expected <- c(
  cxcl9    = 8400  * (1 - 0.983),
  sil2ra   = 6630  * (1 - 0.873),
  ferritin = 15300 * (1 - 0.996)
)

at28 <- typ |> dplyr::filter(time == 28 * 24)
at56 <- typ |> dplyr::filter(time == 56 * 24)

plateau_tab <- data.frame(
  Marker = c("CXCL9 (ng/L)", "sIL-2Ralpha (ng/L)", "Ferritin (ug/L)"),
  `Baseline`                = c(8400, 6630, 15300),
  `Imax (%)`                = c(98.3, 87.3, 99.6),
  `Degradation rate (/day)` = c(0.414, 0.112, 0.207),
  `Plateau = BL x (1-Imax)` = as.numeric(plateau_expected),
  `Simulated at Day 28`     = c(at28$cxcl9, at28$sil2ra, at28$ferritin),
  `Simulated at Day 56`     = c(at56$cxcl9, at56$sil2ra, at56$ferritin),
  check.names = FALSE
)

knitr::kable(
  plateau_tab, digits = 1,
  caption = paste(
    "Typical-value marker concentrations against the structural plateau",
    "baseline x (1 - Imax). Day 56 is the horizon of the source's Figure 3.",
    "Figure 3 shows the CXCL9 median levelling near 100-150 ng/L, ferritin",
    "approaching 50-70 ug/L and sIL-2Ralpha near 1000 ng/L."
  )
)
Typical-value marker concentrations against the structural plateau baseline x (1 - Imax). Day 56 is the horizon of the source’s Figure 3. Figure 3 shows the CXCL9 median levelling near 100-150 ng/L, ferritin approaching 50-70 ug/L and sIL-2Ralpha near 1000 ng/L.
Marker Baseline Imax (%) Degradation rate (/day) Plateau = BL x (1-Imax) Simulated at Day 28 Simulated at Day 56
CXCL9 (ng/L) 8400 98.3 0.4 142.8 144.3 146.6
sIL-2Ralpha (ng/L) 6630 87.3 0.1 842.0 1094.6 855.1
Ferritin (ug/L) 15300 99.6 0.2 61.2 110.3 67.7

How fast each marker reaches that plateau is set by its own degradation rate and by how far its baseline sits above the plateau, so the three converge on very different timescales. Ferritin has both the slowest-but-one turnover (0.207/day) and by far the largest baseline-to-plateau ratio (15,300 -> 61.2, a factor of 250), so it needs roughly 56 days rather than 28 to arrive: at Day 28 it is still 80% above its plateau. This is visible in Figure 3 of the source, where the ferritin median is still descending at Day 28 while CXCL9 has been flat since about Day 14.

# By Day 56 -- the horizon of the source's Figure 3 -- all three markers have
# reached the structural plateau baseline x (1 - Imax).
stopifnot(abs(at56$cxcl9    - plateau_expected[["cxcl9"]])    / plateau_expected[["cxcl9"]]    < 0.20)
stopifnot(abs(at56$sil2ra   - plateau_expected[["sil2ra"]])   / plateau_expected[["sil2ra"]]   < 0.20)
stopifnot(abs(at56$ferritin - plateau_expected[["ferritin"]]) / plateau_expected[["ferritin"]] < 0.20)

# Each marker must still be descending towards the plateau from above, never
# through it.
stopifnot(at56$cxcl9 < at28$cxcl9 * 1.05, at56$sil2ra < at28$sil2ra,
          at56$ferritin < at28$ferritin)
stopifnot(at28$cxcl9    > plateau_expected[["cxcl9"]] * 0.99,
          at28$sil2ra   > plateau_expected[["sil2ra"]],
          at28$ferritin > plateau_expected[["ferritin"]])

A much stricter test is available while dosing continues. With inhibition saturated, the turnover ODE d/dt(m) = kdeg * (plateau - m) has the exact closed-form solution m(t) = plateau + (baseline - plateau) * exp(-kdeg * t). The simulated typical-value trajectory must reproduce it to within a few percent at every time during the treatment period – a far tighter constraint than the 20% plateau band above, and one that would catch any error in the baseline, the degradation rate, the Imax, or the day/hour unit conversion.

closed_form <- function(day, baseline, imax, kdeg) {
  plateau <- baseline * (1 - imax)
  plateau + (baseline - plateau) * exp(-kdeg * day)
}

cf_days <- c(7, 14, 28)
cf_tab <- do.call(rbind, lapply(cf_days, function(dy) {
  r <- typ |> dplyr::filter(time == dy * 24)
  data.frame(
    Day = dy,
    Marker = c("CXCL9", "sIL-2Ralpha", "Ferritin"),
    Simulated = c(r$cxcl9, r$sil2ra, r$ferritin),
    `Closed form` = c(
      closed_form(dy, 8400,  0.983, 0.414),
      closed_form(dy, 6630,  0.873, 0.112),
      closed_form(dy, 15300, 0.996, 0.207)
    ),
    check.names = FALSE
  )
}))
cf_tab$`Difference (%)` <- 100 * (cf_tab$Simulated - cf_tab$`Closed form`) /
  cf_tab$`Closed form`

knitr::kable(
  cf_tab, digits = 2, row.names = FALSE,
  caption = paste(
    "Simulated typical-value marker concentrations against the exact",
    "saturated-inhibition solution of the turnover model during the",
    "treatment period."
  )
)
Simulated typical-value marker concentrations against the exact saturated-inhibition solution of the turnover model during the treatment period.
Day Marker Simulated Closed form Difference (%)
7 CXCL9 600.59 598.05 0.43
7 sIL-2Ralpha 3485.78 3484.67 0.03
7 Ferritin 3643.36 3639.35 0.11
14 CXCL9 169.81 167.90 1.14
14 sIL-2Ralpha 2049.84 2048.58 0.06
14 Ferritin 905.05 901.37 0.41
28 CXCL9 144.32 142.88 1.01
28 sIL-2Ralpha 1094.63 1093.53 0.10
28 Ferritin 110.29 107.52 2.57 

stopifnot(max(abs(cf_tab$`Difference (%)`)) < 5)

Structural check: undosed markers hold at baseline

With no emapalumab given, all three turnover models must sit exactly at their baselines for the whole horizon (production rbase * kdeg exactly balances loss kdeg * marker). This is the steady-state assumption stated in Methods 2.2.

ss_ev <- data.frame(
  id = 1L, WT = 45.5, AGE = 11.5, TBILI = 12.8, IFNG = 500, DIS_MAS = 1,
  time = obs_hours, amt = NA_real_, rate = NA_real_, evid = 0L,
  cmt = "central", dvid = 1L
)

ss <- rxode2::rxSolve(
  modZero, events = ss_ev, omega = NA, sigma = NA,
  returnType = "data.frame"
)

stopifnot(max(abs(ss$cxcl9    / 8400  - 1)) < 1e-6)
stopifnot(max(abs(ss$sil2ra   / 6630  - 1)) < 1e-6)
stopifnot(max(abs(ss$ferritin / 15300 - 1)) < 1e-6)
stopifnot(max(ss$Cc) == 0)

Emapalumab PK profile 

pk_pct <- sim |>
  dplyr::mutate(day = time / 24) |>
  dplyr::group_by(day) |>
  dplyr::summarise(
    lo  = quantile(Cc, 0.05),
    med = median(Cc),
    hi  = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  dplyr::filter(med > 0)

ggplot(pk_pct, aes(day, med)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "steelblue", alpha = 0.3) +
  geom_line(linewidth = 0.8, colour = "steelblue4") +
  scale_y_log10() +
  scale_x_continuous(breaks = seq(0, 140, by = 28)) +
  labs(x = "Time after first dose (days)", y = "Emapalumab (ug/mL)") +
  theme_bw()
Simulated emapalumab serum concentrations in the virtual MAS cohort: median and 5th-95th percentile band. The last dose is on Day 28.

Simulated emapalumab serum concentrations in the virtual MAS cohort: median and 5th-95th percentile band. The last dose is on Day 28.

NCA validation (PKNCA) 

nca_conc <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::transmute(id, time = time / 24, Cc, treatment = "6 mg/kg load + 3 mg/kg maintenance")

nca_dose <- doses |>
  dplyr::transmute(id, time = time / 24, amt,
                   treatment = "6 mg/kg load + 3 mg/kg maintenance")

conc_obj <- PKNCA::PKNCAconc(nca_conc, Cc ~ time | treatment + id,
                             concu = "ug/mL", timeu = "day")
dose_obj <- PKNCA::PKNCAdose(nca_dose, amt ~ time | treatment + id,
                             doseu = "mg")

intervals <- data.frame(
  start      = c(0, 28),
  end        = c(3, 140),
  cmax       = c(TRUE, FALSE),
  tmax       = c(TRUE, FALSE),
  auclast    = c(TRUE, FALSE),
  half.life  = c(FALSE, TRUE),
  lambda.z   = c(FALSE, TRUE)
)

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

nca_long <- as.data.frame(nca_res$result) |>
  dplyr::select(id, treatment, PPTESTCD, PPORRES)

knitr::kable(
  nca_long |>
    dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life")) |>
    dplyr::group_by(PPTESTCD) |>
    dplyr::summarise(
      Median = median(PPORRES, na.rm = TRUE),
      `5th`  = quantile(PPORRES, 0.05, na.rm = TRUE),
      `95th` = quantile(PPORRES, 0.95, na.rm = TRUE),
      .groups = "drop"
    ) |>
    dplyr::rename("NCA parameter" = PPTESTCD),
  digits = 2,
  caption = paste(
    "Simulated NCA in the virtual MAS cohort. Cmax / Tmax / AUClast are over",
    "Days 0-3 (the 6 mg/kg loading infusion); half-life is over Days 28-140",
    "(the terminal washout after the last dose). Cmax in ug/mL, Tmax and",
    "half-life in days, AUClast in day*ug/mL."
  )
)
Simulated NCA in the virtual MAS cohort. Cmax / Tmax / AUClast are over Days 0-3 (the 6 mg/kg loading infusion); half-life is over Days 28-140 (the terminal washout after the last dose). Cmax in ug/mL, Tmax and half-life in days, AUClast in day*ug/mL.
NCA parameter Median 5th 95th
auclast 188.63 135.62 286.48
cmax 115.72 83.53 167.61
half.life 16.13 6.83 45.67
tmax 0.52 0.04 1.00

Comparison against the published terminal half-life

The upstream Brossard & Laveille 2024 analysis reports the median terminal half-life for emapalumab in patients with MAS in sJIA as 24.0 days (range 6.13-32.4), derived from individual post-hoc parameters. The equivalent model-derived quantity is the terminal (beta) half-life of the two-compartment system, computed per virtual subject from that subject’s own micro-constants; the NCA half-life above is reported alongside it.

The two are not the same estimator. The published figure is a median over the 14 real patients using their empirical-Bayes post-hoc parameters, which are shrunk towards the covariate-predicted typical value and reflect the actual observed covariate joint distribution; the simulated figure is a median over a virtual cohort drawn from published marginal summaries. The comparison below is therefore expected to be starred, and it is reported rather than tuned; the tighter reference-patient test that follows localises the difference.

ind <- sim |>
  dplyr::group_by(id) |>
  dplyr::slice(1) |>
  dplyr::ungroup() |>
  dplyr::mutate(
    a    = kel + k12 + k21,
    beta = 0.5 * (a - sqrt(a^2 - 4 * kel * k21)),
    thalf_beta_day = log(2) / beta / 24
  )

nca_half <- nca_long |>
  dplyr::filter(PPTESTCD == "half.life") |>
  dplyr::pull(PPORRES)

simulated_long <- dplyr::bind_rows(
  data.frame(treatment = "MAS in sJIA", PPTESTCD = "thalf_terminal",
             PPORRES = ind$thalf_beta_day),
  data.frame(treatment = "MAS in sJIA", PPTESTCD = "thalf_nca",
             PPORRES = nca_half)
)

published <- tibble::tribble(
  ~treatment,     ~thalf_terminal, ~thalf_nca,
  "MAS in sJIA",  24.0,            24.0
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = simulated_long,
  reference = published,
  by = "treatment",
  units = c(thalf_terminal = "day", thalf_nca = "day"),
  tolerance_pct = 20
)
#> Warning: ncaParamLabel(): unknown PKNCA code(s) returned as-is:
#> 'thalf_terminal', 'thalf_nca'

cmp |>
  dplyr::rename("Population" = treatment) |>
  knitr::kable(
    digits = 2,
    caption = paste(
      "Simulated emapalumab terminal half-life against the 24.0-day median",
      "reported for patients with MAS in sJIA by Brossard & Laveille 2024.",
      "thalf_terminal is the analytic beta half-life of the two-compartment",
      "model; thalf_nca is the PKNCA lambda-z half-life over Days 28-140.",
      "* marks a difference greater than 20%."
    )
  )
Simulated emapalumab terminal half-life against the 24.0-day median reported for patients with MAS in sJIA by Brossard & Laveille 2024. thalf_terminal is the analytic beta half-life of the two-compartment model; thalf_nca is the PKNCA lambda-z half-life over Days 28-140. * marks a difference greater than 20%.
NCA parameter Population Reference Simulated % diff
thalf_nca (day) MAS in sJIA 24 16.1 -32.8%*
thalf_terminal (day) MAS in sJIA 24 16.1 -32.7%* 
# The published summary is a MEDIAN OVER 14 REAL PATIENTS derived from
# individual post-hoc (empirical Bayes) parameters, with a wide 6.13-32.4 day
# range. The comparable simulated quantity is the median over the virtual
# cohort. These are not the same estimator, so the check is that the simulated
# median falls inside the published range rather than close to its median.
stopifnot(median(ind$thalf_beta_day) > 6.13,
          median(ind$thalf_beta_day) < 32.4)
# The cohort must also span a realistic spread rather than collapsing onto a
# single value.
stopifnot(quantile(ind$thalf_beta_day, 0.05) < 12,
          quantile(ind$thalf_beta_day, 0.95) > 25)

The simulated cohort median sits below the published 24.0-day median, which the sharper check below resolves: the offset is in the terminal half-life itself, not in clearance.

Second falsifier: the upstream four-point terminal half-life

The same Brossard & Laveille 2024 sentence that supplies the four clearances also gives the corresponding terminal half-lives – 19.2, 13.8, 7.18 and 3.12 days – for the same reference patient. This is a strictly stronger test than the clearance key, because the terminal half-life depends on V1, V2 and Q as well as on CL, so it exercises the whole disposition structure rather than clearance alone.

key_a    <- key_sim$kel + key_sim$k12 + key_sim$k21
key_beta <- 0.5 * (key_a - sqrt(key_a^2 - 4 * key_sim$kel * key_sim$k21))
key_thalf <- log(2) / key_beta / 24

published_thalf <- c(19.2, 13.8, 7.18, 3.12)

knitr::kable(
  data.frame(
    `Total IFN-gamma (pg/mL)` = format(ifng_levels, scientific = TRUE),
    `Simulated t1/2 (day)`    = key_thalf,
    `Published t1/2 (day)`    = published_thalf,
    `Difference (%)`          = 100 * (key_thalf - published_thalf) / published_thalf,
    check.names = FALSE
  ),
  digits = c(0, 2, 2, 1),
  caption = paste(
    "Simulated terminal (beta) half-life for the upstream reference patient",
    "against the four values in the Brossard & Laveille 2024 Results text."
  )
)
Simulated terminal (beta) half-life for the upstream reference patient against the four values in the Brossard & Laveille 2024 Results text.
Total IFN-gamma (pg/mL) Simulated t1/2 (day) Published t1/2 (day) Difference (%)
1e+03 16.95 19.20 -11.7
1e+04 12.20 13.80 -11.6
1e+05 6.43 7.18 -10.4
1e+06 2.92 3.12 -6.4 

# The model runs a systematic 6-12% below the published terminal half-lives
# while reproducing clearance to within 0.15%, so the offset is in the
# distribution parameters, not in CL. It is NOT tuned away -- see Errata.
stopifnot(all(abs(key_thalf - published_thalf) / published_thalf < 0.15))
stopifnot(all(diff(key_thalf) < 0))   # half-life must fall as IFN-gamma rises

The half-life falls monotonically from 17 to 3 days across four orders of magnitude of total IFN-gamma – the target-mediated disposition that the paper identifies as absent in MAS, where CLNL is switched off and the half-life stays long.

The model reproduces the published PK without any tuning: the linear clearance in MAS gives a long terminal half-life, matching the paper’s statement that MAS half-life is “close to that observed in healthy subjects” rather than shortened by target-mediated disposition as in primary HLH.

Assumptions and deviations

  • Age reference (25 years) is derived, not published. Neither Brossard 2025 nor the upstream Brossard & Laveille 2024 prints the age at which the typical values apply. It was recovered by back-solving the upstream four-point clearance answer key: the 1e6 / 1e3 IFN-gamma ratio cancels the age and weight factors and returns a total-bilirubin factor of 1.0004, after which the age factor solves to (1 / 25.3)^0.188. Rounding the reference to 25 years reproduces all four published clearances to within 0.3% (see the answer-key section above). The competing reading, centring at the dataset median age of 1.7 years, is off by 83% and is excluded.
  • Total-bilirubin reference (12.8 umol/L) is derived, not published. The back-solve places the reference patient at a bilirubin factor of 1.0004, identifying the reference as the pooled-dataset median that Table S1 reports identically (12.8 umol/L) for primary HLH, MAS and All.
  • IIV values are read as omega standard deviations on the log scale. Table 1 reports 0.361 / 0.207 / 0.711 for CL / V1 / V2; back-transforming with CV = sqrt(exp(sd^2) - 1) reproduces Table S2’s 37.3 / 20.9 / 81.2 CV% exactly, which fixes the scale. Table 2’s PD variabilities are read on the same convention, by consistency within the paper; Table 2 does not restate them as CV%, so this is inferred rather than proven for the PD block.
  • Only three PD correlations are reported. Table 2 gives correlations for CXCL9 baseline with baseline ferritin (0.675), sIL-2Ralpha baseline with its own Imax (0.805), and ferritin baseline with its own Imax (0.684). All other off-diagonal covariances are set to zero. All three resulting blocks are positive definite.
  • Residual error. Table 1 reports the emapalumab RUV as an “additive component, log” of 0.301, which is proportional error in nlmixr2’s linear space. Table 2’s PD RUVs are already stated as proportions.
  • Dosing calendar. “Twice weekly until Day 28” is implemented as Days 18, 22, 25 and 28. The exact calendar is not printed. Infusion duration is taken as 1 hour, as stated for the primary-HLH regimen in Methods 2.1.
  • Virtual cohort covariate joint distribution. Table S1 gives marginal medians and ranges only. Age and weight are rank-coupled here so that older patients are heavier; total bilirubin is drawn independently from a truncated log-normal matched to its published median and mean.
  • Original data are not public. The Data Availability Statement offers participant-level data only on request to Sobi, so all figures here are simulations from the packaged model rather than overlays on observed data.

Errata and source discrepancies

  • sIL-2Ralpha baseline: 6630 vs 6550 ng/L. Table 2 reports the original estimate as 6630 ng/L (bootstrap mean 6700, 95% CI 4140-9260). The abstract, the Results narrative and the Figure 1 schematic all quote 6550 ng/L. The model file uses the parameter table’s 6630, on the basis that Table 2 is the authoritative parameter listing and is the only one of the four locations that carries an RSE and a bootstrap interval. Users who prefer the narrative value should set lrbase_sil2ra <- log(6550); the difference is 1.2% and does not affect any check in this vignette.
  • Age exponent: 0.188 vs 0.187. Table 1 (the MAS re-estimation) reports 0.188; Table S2 and the upstream final model report 0.187. The model file uses Table 1’s 0.188 because Table 1 is this paper’s own re-estimated MAS model. The difference is immaterial (it changes the reference-patient clearance by less than 0.1%).
  • Table S2 duplicates the upstream final model. Every value and every RSE in Brossard 2025 Table S2 is identical to the final-model table of Brossard & Laveille 2024. Table S2 is therefore a reproduction of the pooled fit, not an independent MAS-only fit; Table 1 is the MAS re-estimate.
  • No potency parameter is published for the PD link. This is the one genuine gap. Table 2 lists exactly four parameter groups per marker - Baseline, Degradation rate, Imax and RUV - with no IC50, EC50, slope or Hill term, and the paper prints no equations at all. Methods 2.2 nevertheless states that “emapalumab concentrations were assumed to affect (inhibit) PD marker production rates”, and Figure 1 draws the central compartment acting on all three production arrows, so a concentration-driven link exists but its potency was never reported. The model file therefore fixes lec50 at 0.0248 ug/mL, imported from the upstream primary-HLH analysis (Jacqmin et al. 2022, Br J Clin Pharmacol 88:2128-2139, Supplementary Table 7, IC50 = 24.8 ng/mL) and not from this paper. One shared potency serves all three markers, with the three published Imax values (98.3%, 87.3%, 99.6%) carrying the marker-specific extent of inhibition; it is held fixed() so it cannot be mistaken for a fit to this dataset. Emapalumab troughs in this regimen sit roughly three orders of magnitude above that value, so over the observed 0-56 day window this is numerically indistinguishable from treating inhibition as fully saturated; the two readings diverge only when extrapolating to lower doses or well into washout. Any use of this model for dose exploration below the studied regimen rests on that imported potency, not on Brossard 2025.
  • Terminal half-life runs 6-12% below the published values while clearance matches to 0.15%. For the upstream reference patient the model returns terminal half-lives of 16.9 / 12.2 / 6.4 / 2.9 days against the published 19.2 / 13.8 / 7.18 / 3.12, a systematic negative offset, even though the four clearances at the same four IFN-gamma levels are reproduced to within 0.15%. Because clearance is exact, the offset lies in the distribution parameters (V1, V2, Q), not in CL. The covariate structure was checked against the upstream Supplementary Table 1 run log and is correct as implemented: run 013 selected “same age effect for CL/Q, and no age effect on V2”, and run 019 (the final model) applied “same TBIL effect on CLL, V1 & V2”. Adding an age effect on V2 would close most of the gap but is explicitly excluded by run 013, so it was not done. The most likely explanation is that the published half-lives were derived by a different route (post-hoc individual parameters or a numerical fit to the simulated terminal phase) rather than as the analytic beta half-life of the typical-value reference patient. No parameter was tuned to close this gap.
  • The upstream main paper had to be re-acquired during this extraction. The copy of Brossard & Laveille 2024 (PMC11111609) in the source directory was a 46-byte failed-acquisition JSON stub rather than a PDF, so the four-point clearance key, the four-point half-life key and the 24.0-day MAS median could not be verified against it as first drafted. The paper was re-fetched from the publisher and all of those values were confirmed verbatim before being relied on here. Its supplement (the run log) was present and valid throughout.
  • No erratum found. A search of the journal’s correction notices and of PubMed for corrections to doi:10.1111/cts.70163 returned nothing as of the extraction date.