Skip to contents

Trofinetide, DAFFODIL pediatric 2-4 years (Darwish 2025b)

Source: vignettes/articles/Darwish_2025b_trofinetide.Rmd
Darwish_2025b_trofinetide.Rmd

Model and source

  • Citation: Darwish M, Passarell J, Maxwell K, Bradley H, Bishop KM, Youakim JM. Population Pharmacokinetics of Trofinetide in a Pediatric Population Aged 2 to 4 Years with Rett Syndrome. Advances in Therapy. 2025;42(2):1009-1025. doi:10.1007/s12325-024-03058-7
  • Description: Updated population PK model for oral trofinetide in Rett syndrome (Darwish 2025b, DAFFODIL): two-compartment with first-order absorption and linear elimination, re-estimated after adding pediatric data from girls aged 2-4 years to the 13-study pool behind Darwish_2025a_trofinetide.
  • Article: Adv Ther. 2025;42(2):1009-1025

Trofinetide is the first approved treatment for Rett syndrome (RTT), dosed orally twice daily on a body-weight-banded schedule. Acadia and Simulations Plus built a 13-study population pharmacokinetic model to support the weight bands used in the phase 3 LAVENDER study (girls and women aged 5-20 years); that model is packaged separately as Darwish_2025a_trofinetide.

This vignette covers the update: the same model re-estimated after adding 114 whole-blood concentrations from 13 girls aged 2-4 years enrolled in the phase 2/3 DAFFODIL study, used to confirm that the DAFFODIL weight bands (5 g BID for >= 9 to < 12 kg, 6 g BID for >= 12 to < 20 kg) put steady-state exposure inside the 800-1200 ug*h/mL target window.

A transcription hazard specific to this paper

The Table 2 footnote of Darwish 2025b prints a block headed “Model equations:” that gives typical-value expressions for F1, ka, CL, Vc, and Vp. Those equations carry the earlier model’s coefficients, not the ones in the table they annotate. The footnote reads F1 = 0.828 ..., ka = 0.391 ..., CL = 11.8 x (WTKG/58)^0.443 ..., Vc = 24.9 ..., Vp = 35.3 ..., all of which are the Darwish 2025a values; the “Population mean estimate” column beside them reads 0.832, 0.394, 11.7, 0.486, 25.0, and 35.4.

The table wins, and the paper’s own prose proves it – see the cross-model gate below, which reproduces the two published change figures exactly. The model file takes only the functional form from the footnote and every coefficient from the table.

Population

The pooled analysis dataset comprised 5709 trofinetide whole-blood concentrations from 455 participants across 14 clinical studies: eight phase 1 studies, four phase 2 studies in RTT, fragile X syndrome (FXS) and traumatic brain injury (TBI), the phase 3 LAVENDER study, and the phase 2/3 DAFFODIL study. The cohort was 156 healthy volunteers, 198 patients with RTT (female only), 57 with TBI (male only), and 44 with FXS (male only), predominantly female (56.9%), with a mean age of 21.8 years (range 2-64) and mean body weight of 56.5 kg (range 9.8-140).

The DAFFODIL subgroup that motivated the update was 13 girls (92.3% white, 7.7% Asian) with a mean age of 3 years (range 2-4) and mean baseline body weight of 13.4 kg (range 9.8-18.1). Concentrations were measured by LC-MS/MS in lithium-heparinized whole blood, LLOQ 0.100 ug/mL and ULOQ 100 ug/mL.

str(mod_meta$population)
#> List of 12
#>  $ species       : chr "human"
#>  $ n_subjects    : num 455
#>  $ n_observations: num 5709
#>  $ n_studies     : num 14
#>  $ age_range     : chr "2-64 years (mean 21.8)"
#>  $ weight_range  : chr "9.8-140 kg (mean 56.5)"
#>  $ sex_female_pct: num 56.9
#>  $ gfr_reference : chr "124 mL/min/1.73 m^2 (analysis-population median)"
#>  $ disease_state : chr "Pooled analysis of 156 healthy volunteers, 198 patients with Rett syndrome (female only), 57 patients with trau"| __truncated__
#>  $ dose_range    : chr "Oral, gastric-tube, and intravenous (bolus and infusion) trofinetide; oral doses spanned the 2-12 g therapeutic"| __truncated__
#>  $ regions       : chr "Not reported"
#>  $ notes         : chr "Darwish 2025b Results, 'Population Pharmacokinetic Model'. This is the 14-study update of the 13-study model pa"| __truncated__

Source trace

Every value below is transcribed from the “Population mean estimate” column of Darwish 2025b Table 2. The reference values 58 kg, 22.4 years, and 124 mL/min/1.73 m^2 are the analysis-population medians per the Table 2 footnote. The form of each covariate term – a proportional shift (1 + theta * indicator) and power terms on WT, GFR and AGE – comes from the footnote equations; their coefficients do not (see above).

Equation / parameter Value Source location
lcl (CL) 11.7 L/h Table 2, “CL (L/h)”
e_wt_cl 0.486 Table 2, “Covariate exponent of weight on CL”
e_crcl_cl 0.272 Table 2, “Covariate exponent of GFR on CL”
e_rett_cl -0.146 Table 2, “Shift in CL for RTT = 1”
e_tbi_cl 0.229 Table 2, “Shift in CL for TBI = 1”
lvc (Vc) 25.0 L Table 2, “Vc (L)”
e_age_vc 0.556 Table 2, “Covariate exponent of age on Vc”
e_fxs_vc 1.16 Table 2, “Shift in Vc for FXS = 1”
lq (Q) 1.42 L/h Table 2, “Q (L/h)”
lvp (Vp) 35.4 L Table 2, “Vp (L)”
e_rett_vp 0.805 Table 2, “Shift in Vp for RTT = 1”
e_tbi_vp -0.753 Table 2, “Shift in Vp for TBI = 1”
lka (ka) 0.394 1/h Table 2, “ka (1/h)”
e_fed_ka -0.0969 Table 2, “Shift in ka for FED”
lfdepot (F1) 0.832 Table 2, “F1”
e_fed_f -0.133 Table 2, “Shift in F1 for FED”
e_dose18g_f -0.132 Table 2, “Shift in F1 for 18 g dose”
e_dose24g_f -0.284 Table 2, “Shift in F1 for 24 g dose”
e_diarrhea_f -0.157 Table 2, “Shift in F1 for diarrhea”
etalcl 0.0182 Table 2, “Interindividual variability / CL”; printed as 13.6 %CV
etalvc 0.0906 Table 2, “Interindividual variability / Vc”; printed as 30.8 %CV
etalq 0.355 Table 2, “Interindividual variability / Q”; printed as 65.3 %CV
etalvp 0.0874 Table 2, “Interindividual variability / Vp”; printed as 30.2 %CV
etalfdepot 0.04 Table 2, “Interindividual variability / F1”; printed as 20.2 %CV
expSdHealthy sqrt(0.0788) = 0.281 Table 2, “Residual variability / Healthy subjects”; printed as 28.1 %CV
expSdDisease sqrt(0.140) = 0.374 Table 2, “Residual variability / Patients with RTT, TBI or FXS”; printed as 37.4 %CV
Structural model (2-cmt, first-order absorption, linear elimination) n/a Methods, “Population Pharmacokinetic Model Development”
Proportional-shift form (1 + theta * indicator) n/a Table 2 footnote equations (form only)
Exponential IIV, %CV = sqrt(exp(omega^2) - 1) * 100 n/a Methods Eq. (1)
Log/exponential residual error n/a Methods Eq. (2)

Two transcription points are worth spelling out, because Table 2 reports each variability quantity twice and the two columns are on different scales.

The interindividual-variability block holds variances. The block near the bottom of Table 2 lists 0.0182, 0.0906, 0.355, 0.0874, and 0.04 with no units. Methods Eq. (1) defines %CV = sqrt(exp(omega^2) - 1) * 100; applying it to those five numbers returns 13.6, 30.8, 65.3, 30.2, and 20.2 %CV, which is exactly the “Variability / Estimate” column beside each structural parameter. They are therefore the NONMEM $OMEGA variances and are used directly.

The residual-variability rows hold variances too, but with a different %CV convention. The tabulated 0.0788 and 0.140 map to the printed 28.1 %CV and 37.4 %CV through the plain sqrt(), not through Eq. (1) (sqrt(0.0788) = 0.281, sqrt(0.140) = 0.374). lnorm() takes a log-scale SD, so the model file enters sqrt(0.0788) and sqrt(0.140).

# Reproduce both %CV conventions from the tabulated variances, confirming which
# column is which. Deterministic; no simulation involved.
iiv_var <- c(CL = 0.0182, Vc = 0.0906, Q = 0.355, Vp = 0.0874, F1 = 0.04)
rv_var  <- c(Healthy = 0.0788, Disease = 0.140)

check_var <- tibble::tibble(
  Quantity = c(names(iiv_var), names(rv_var)),
  Block = c(rep("Interindividual", length(iiv_var)), rep("Residual", length(rv_var))),
  `Table 2 variance` = c(iiv_var, rv_var),
  `Recomputed %CV` = 100 * c(sqrt(exp(iiv_var) - 1), sqrt(rv_var)),
  `Table 2 printed %CV` = c(13.6, 30.8, 65.3, 30.2, 20.2, 28.1, 37.4)
)

knitr::kable(check_var, digits = c(0, 0, 4, 1, 1),
             caption = "Both Table 2 variability columns recovered from the tabulated variances.")
Both Table 2 variability columns recovered from the tabulated variances.
Quantity Block Table 2 variance Recomputed %CV Table 2 printed %CV
CL Interindividual 0.0182 13.6 13.6
Vc Interindividual 0.0906 30.8 30.8
Q Interindividual 0.3550 65.3 65.3
Vp Interindividual 0.0874 30.2 30.2
F1 Interindividual 0.0400 20.2 20.2
Healthy Residual 0.0788 28.1 28.1
Disease Residual 0.1400 37.4 37.4 

stopifnot(max(abs(check_var$`Recomputed %CV` - check_var$`Table 2 printed %CV`)) < 0.06)
ui <- mod_meta
ui$iniDf[, c("name", "est", "label")] |>
  dplyr::rename("Parameter" = name, "Estimate" = est, "Label" = label) |>
  knitr::kable(digits = 4, caption = "Packaged model parameters.")
Packaged model parameters.
Parameter Estimate Label
lcl 2.4596 Clearance at the reference covariate values (L/h)
e_wt_cl 0.4860 Power exponent on (WT/58) for clearance (unitless)
e_crcl_cl 0.2720 Power exponent on (CRCL/124) for clearance (unitless)
e_rett_cl -0.1460 Proportional shift in clearance for Rett syndrome (fraction)
e_tbi_cl 0.2290 Proportional shift in clearance for traumatic brain injury (fraction)
lvc 3.2189 Central volume of distribution at the reference age (L)
e_age_vc 0.5560 Power exponent on (AGE/22.4) for central volume (unitless)
e_fxs_vc 1.1600 Proportional shift in central volume for fragile X syndrome (fraction)
lq 0.3507 Intercompartmental clearance (L/h)
lvp 3.5667 Peripheral volume of distribution (L)
e_rett_vp 0.8050 Proportional shift in peripheral volume for Rett syndrome (fraction)
e_tbi_vp -0.7530 Proportional shift in peripheral volume for traumatic brain injury (fraction)
lka -0.9314 First-order absorption rate constant in the fasted state (1/h)
e_fed_ka -0.0969 Proportional shift in ka for the fed state (fraction)
lfdepot -0.1839 Oral bioavailability in the fasted, therapeutic-dose, diarrhea-free reference state (fraction)
e_fed_f -0.1330 Proportional shift in bioavailability for the fed state (fraction)
e_dose18g_f -0.1320 Proportional shift in bioavailability for the 18 g dose level (fraction)
e_dose24g_f -0.2840 Proportional shift in bioavailability for the 24 g dose level (fraction)
e_diarrhea_f -0.1570 Proportional shift in bioavailability during diarrhea (fraction)
expSdHealthy 0.2807 Log-scale residual SD, healthy subjects (log units)
expSdDisease 0.3742 Log-scale residual SD, subjects with Rett syndrome, TBI, or FXS (log units)
etalcl 0.0182 Table 2 IIV block: omega^2 for CL = 0.0182, printed as 13.6 %CV (RSE 17.4%; eta-shrinkage 37.5%)
etalvc 0.0906 Table 2 IIV block: omega^2 for Vc = 0.0906, printed as 30.8 %CV (RSE 12.5%; eta-shrinkage 31.9%)
etalq 0.3550 Table 2 IIV block: omega^2 for Q = 0.355, printed as 65.3 %CV (RSE 17.7%; eta-shrinkage 9.8%)
etalvp 0.0874 Table 2 IIV block: omega^2 for Vp = 0.0874, printed as 30.2 %CV (RSE 14.3%; eta-shrinkage 41.0%)
etalfdepot 0.0400 Table 2 IIV block: omega^2 for F1 = 0.04, printed as 20.2 %CV (RSE 16.0%; eta-shrinkage 54.8%)

Cross-model gate: reproducing the paper’s own update summary

This is the sharpest available check on the transcription, and it needs no simulation at all. Darwish 2025b Results makes two quantitative claims about how the parameters moved when the DAFFODIL data were added to the 13-study pool:

  1. “CL, Q, Vc, Vp, first-order absorption rate constant (ka), and F1 differing by < 1.5% between the two models.”
  2. “Covariate effect estimates were generally minimally affected (< 10% change) …; notable exceptions were the shift in CL and Vp with positive RTT disease status, which increased by 13.6% and 30.7%, respectively, compared with the previous model.”

Both models are packaged here, so the claims can be evaluated directly against the two ini() blocks. If the update had been transcribed from the stale footnote equations instead of from Table 2, every one of these deltas would be exactly zero and the gate would fail.

est_of <- function(m) setNames(m$iniDf$est, m$iniDf$name)
new <- est_of(mod_meta)
old <- est_of(prev_meta)

pct_change <- function(nm, backtransform = FALSE) {
  a <- old[[nm]]; b <- new[[nm]]
  if (backtransform) { a <- exp(a); b <- exp(b) }
  100 * (b - a) / abs(a)
}

primary <- c(CL = "lcl", Q = "lq", Vc = "lvc", Vp = "lvp", ka = "lka", F1 = "lfdepot")
cov_eff <- c("e_wt_cl", "e_crcl_cl", "e_rett_cl", "e_tbi_cl", "e_age_vc",
             "e_fxs_vc", "e_rett_vp", "e_tbi_vp", "e_fed_ka", "e_fed_f",
             "e_dose18g_f", "e_dose24g_f", "e_diarrhea_f")

primary_tbl <- tibble::tibble(
  Parameter = names(primary),
  `Darwish 2025a` = exp(old[primary]),
  `Darwish 2025b` = exp(new[primary]),
  `% change` = vapply(primary, pct_change, numeric(1), backtransform = TRUE)
)

cov_tbl <- tibble::tibble(
  `Covariate effect` = cov_eff,
  `Darwish 2025a` = old[cov_eff],
  `Darwish 2025b` = new[cov_eff],
  `% change` = vapply(cov_eff, pct_change, numeric(1))
)

knitr::kable(primary_tbl, digits = c(0, 4, 4, 2),
             caption = "Primary parameters: the paper states all six move by < 1.5%.")
Primary parameters: the paper states all six move by < 1.5%.
Parameter Darwish 2025a Darwish 2025b % change
CL 11.800 11.700 -0.85
Q 1.440 1.420 -1.39
Vc 24.900 25.000 0.40
Vp 35.300 35.400 0.28
ka 0.391 0.394 0.77
F1 0.828 0.832 0.48 
knitr::kable(cov_tbl, digits = c(0, 4, 4, 2),
             caption = "Covariate effects: the paper states all move by < 10% except the two RTT shifts.")
Covariate effects: the paper states all move by < 10% except the two RTT shifts.
Covariate effect Darwish 2025a Darwish 2025b % change
e_wt_cl 0.4430 0.4860 9.71
e_crcl_cl 0.2730 0.2720 -0.37
e_rett_cl -0.1690 -0.1460 13.61
e_tbi_cl 0.2350 0.2290 -2.55
e_age_vc 0.5490 0.5560 1.28
e_fxs_vc 1.1500 1.1600 0.87
e_rett_vp 0.6160 0.8050 30.68
e_tbi_vp -0.7520 -0.7530 -0.13
e_fed_ka -0.0949 -0.0969 -2.11
e_fed_f -0.1330 -0.1330 0.00
e_dose18g_f -0.1320 -0.1320 0.00
e_dose24g_f -0.2840 -0.2840 0.00
e_diarrhea_f -0.1480 -0.1570 -6.08 
# Claim 1: all six primary parameters within 1.5%.
stopifnot(nrow(primary_tbl) == 6L)
stopifnot(all(abs(primary_tbl$`% change`) < 1.5))

# Claim 2: exactly two covariate effects exceed 10%, and they are the RTT
# shifts on CL and Vp.
exceptions <- cov_tbl$`Covariate effect`[abs(cov_tbl$`% change`) >= 10]
stopifnot(identical(sort(exceptions), c("e_rett_cl", "e_rett_vp")))

# Claim 2, quantitative: those two increased by 13.6% and 30.7%.
delta <- setNames(cov_tbl$`% change`, cov_tbl$`Covariate effect`)
stopifnot(abs(delta[["e_rett_cl"]] - 13.6) < 0.1)
stopifnot(abs(delta[["e_rett_vp"]] - 30.7) < 0.1)

cat(sprintf(
  "Reproduced: RTT shift on CL %+.1f%% (paper: +13.6%%); RTT shift on Vp %+.1f%% (paper: +30.7%%).\n",
  delta[["e_rett_cl"]], delta[["e_rett_vp"]]
))
#> Reproduced: RTT shift on CL +13.6% (paper: +13.6%); RTT shift on Vp +30.7% (paper: +30.7%).

Both published change figures come back to three significant figures. That is only possible if the coefficients were read from Table 2, which settles the footnote conflict flagged at the top of this vignette.

Virtual cohort

Individual-level data are not public (Darwish 2025b Data Availability). The cohorts below approximate the two studies compared in the paper’s Figure 3: the DAFFODIL bands for girls aged 2-4 years, and the four LAVENDER bands for ages 5-20 years. Weight-band boundaries and doses are from Supplementary Material Table S2.

Following the paper’s own exposure simulations, all subjects have RTT, are fasted, are free of diarrhea, and receive therapeutic (not supratherapeutic) doses.

tau <- 12 # h, BID dosing interval

bands <- tibble::tribble(
  ~study,      ~band,          ~wt_lo, ~wt_hi, ~dose_g, ~age_lo, ~age_hi, ~gfr_mean,
  "DAFFODIL",  "9-12 kg",         9,     12,      5,       2,       4,      150,
  "DAFFODIL",  "12-20 kg",       12,     20,      6,       2,       4,      150,
  "LAVENDER",  "12-20 kg (L)",   12,     20,      6,       5,       8,      162,
  "LAVENDER",  "20-35 kg",       20,     35,      8,       7,      13,      162,
  "LAVENDER",  "35-50 kg",       35,     50,     10,      11,      17,      162,
  "LAVENDER",  ">50 kg",         50,     70,     12,      14,      20,      162
) |>
  dplyr::mutate(amt_mg = dose_g * 1000)

# Keep `band` a plain character column everywhere (rxSolve `keep=`, PKNCA
# grouping, and the joins below are all type-sensitive); apply the display
# ordering only at plotting time.
band_levels <- bands$band
daffodil_bands <- bands$band[bands$study == "DAFFODIL"]
set.seed(20250907)
n_per_band <- 100L # well under the 200/arm cap

make_cohort <- function(band_row, n, id_offset) {
  subj <- tibble::tibble(
    id = id_offset + seq_len(n),
    study = band_row$study,
    band = band_row$band,
    amt_mg = band_row$amt_mg,
    WT = stats::runif(n, band_row$wt_lo, band_row$wt_hi),
    AGE = stats::runif(n, band_row$age_lo, band_row$age_hi),
    CRCL = pmin(pmax(stats::rnorm(n, band_row$gfr_mean, 27), 80), 220),
    DIS_RETT = 1, DIS_TBI = 0, DIS_FXS = 0,
    FED = 0, AE_DIARRHEA = 0, DOSE_18G = 0, DOSE_24G = 0
  )

  # Steady-state dose record. `ss = 1` with `ii = tau` establishes steady state
  # directly; trofinetide's disposition is slow relative to the 12 h interval,
  # so an explicit multiple-dose run would need weeks of simulated doses.
  dosing <- subj |>
    dplyr::mutate(time = 0, amt = amt_mg, evid = 1L, cmt = "depot",
                  ii = tau, ss = 1L)

  # Observations are placed on the `central` ODE state; rxode2 returns the
  # algebraic observable `Cc` as a column at those records.
  obs <- subj |>
    tidyr::crossing(time = seq(0, tau, by = 0.25)) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = "central",
                  ii = 0, ss = 0L)

  dplyr::bind_rows(dosing, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(
  lapply(seq_len(nrow(bands)), function(i) {
    make_cohort(bands[i, ], n_per_band, id_offset = (i - 1L) * n_per_band)
  })
)

stopifnot(!anyDuplicated(events[, c("id", "time", "evid")]))
stopifnot(nrow(dplyr::distinct(events, id)) == n_per_band * nrow(bands))

Simulation 

mod <- readModelDb("Darwish_2025b_trofinetide")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep = c("study", "band", "amt_mg")
) |>
  as.data.frame() |>
  # rxSolve returns observation records only (there is no `evid` column in the
  # output) and may return character keeps as factors; force them back to
  # character so the joins against the PKNCA results below are type-stable.
  dplyr::mutate(band = as.character(band), study = as.character(study))
#> ℹ parameter labels from comments will be replaced by 'label()'

stopifnot(!anyNA(sim$Cc), nrow(sim) > 0)

Replicating Figure 4: individual steady-state profiles vs the target window

Figure 4 of Darwish 2025b overlays the 13 DAFFODIL participants’ model-predicted steady-state profiles on the average curves corresponding to AUC0-12 values of 800 and 1200 ug*h/mL, with a median curve for the subjects inside that range. Those reference curves are, by construction, flat average concentrations of AUC0-12 / 12.

target_lines <- tibble::tibble(
  label = c("AUC0-12 = 800", "AUC0-12 = 1000", "AUC0-12 = 1200"),
  cav = c(800, 1000, 1200) / tau
)

sim |>
  dplyr::filter(band %in% daffodil_bands) |>
  dplyr::mutate(band = factor(band, levels = daffodil_bands)) |>
  dplyr::group_by(band, time) |>
  dplyr::summarise(
    Q10 = stats::quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = stats::quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  geom_hline(data = target_lines, aes(yintercept = cav, colour = label),
             linetype = "dashed") +
  facet_wrap(~band) +
  labs(
    x = "Time since dose (h)", y = "Trofinetide whole blood (ug/mL)",
    colour = NULL,
    title = "DAFFODIL steady-state profiles by weight band",
    caption = "Replicates Figure 4 of Darwish 2025b (median with 10th-90th percentile band)."
  ) +
  theme(legend.position = "bottom")

PKNCA validation 

sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, band)

# Guarantee a time = 0 record per (id, band). The ss = 1 dose means the
# time = 0 concentration is the steady-state trough, which the solve already
# produces; this bind_rows is a defensive no-op that `distinct()` collapses.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, band) |> dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, band, time, .keep_all = TRUE) |>
  dplyr::arrange(id, band, time)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | band + id,
  concu = "ug/mL", timeu = "h"
)

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

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | band + id, doseu = "mg")

intervals <- data.frame(
  start = 0, end = tau,
  cmax = TRUE, tmax = TRUE, ctrough = TRUE, auclast = TRUE, cav = TRUE
)

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

nca_auc <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD == "auclast") |>
  dplyr::select(id, band, auc_pknca = PPORRES) |>
  dplyr::mutate(band = as.character(band))
stopifnot(nrow(nca_auc) == n_per_band * nrow(bands))

Closed-form clearance gate

At steady state the exposure over a dosing interval reduces to AUC0-tau = F1 * Dose / CL, independent of the distribution parameters. This exercises the weight, GFR and disease-state covariate terms on CL and the whole bioavailability chain, and it must agree with the numerically integrated PKNCA result. Because both CL and F1 carry log-normal IIV, the median of the simulated AUC converges on the typical-value closed form.

th <- setNames(ui$iniDf$est, ui$iniDf$name)

cl_typical <- function(WT, CRCL, DIS_RETT = 1, DIS_TBI = 0) {
  exp(th[["lcl"]]) * (WT / 58)^th[["e_wt_cl"]] * (CRCL / 124)^th[["e_crcl_cl"]] *
    (1 + th[["e_tbi_cl"]] * DIS_TBI) * (1 + th[["e_rett_cl"]] * DIS_RETT)
}

gate <- sim |>
  dplyr::distinct(id, study, band, WT, CRCL, amt_mg) |>
  dplyr::mutate(auc_closed_form = exp(th[["lfdepot"]]) * amt_mg / cl_typical(WT, CRCL)) |>
  dplyr::inner_join(nca_auc, by = c("id", "band")) |>
  dplyr::group_by(study, band) |>
  dplyr::summarise(
    `Closed form F1*Dose/CL` = stats::median(auc_closed_form),
    `PKNCA AUC0-12` = stats::median(auc_pknca),
    .groups = "drop"
  ) |>
  dplyr::mutate(`% diff` = 100 * (`PKNCA AUC0-12` - `Closed form F1*Dose/CL`) /
                  `Closed form F1*Dose/CL`)

knitr::kable(gate, digits = 1,
             caption = "Median closed-form vs. numerically integrated steady-state AUC0-12.")
Median closed-form vs. numerically integrated steady-state AUC0-12.
study band Closed form F1*Dose/CL PKNCA AUC0-12 % diff
DAFFODIL 12-20 kg 917.7 902.9 -1.6
DAFFODIL 9-12 kg 905.1 866.4 -4.3
LAVENDER 12-20 kg (L) 864.4 846.2 -2.1
LAVENDER 20-35 kg 905.7 863.1 -4.7
LAVENDER 35-50 kg 914.5 936.7 2.4
LAVENDER >50 kg 915.6 940.4 2.7 

# Both sides use the same drawn parameters, so the difference here is numerical
# (solver tolerance plus the median-of-a-ratio-of-lognormals offset), not
# physical. A tight bound is the correct assertion for this comparison.
stopifnot(nrow(gate) == nrow(bands))
stopifnot(all(abs(gate$`% diff`) < 10))

Comparison against the published target exposure window

Darwish 2025b does not tabulate per-band AUC values; they appear only in the Figure 3 box plots. What the paper does state numerically is its exposure-matching criterion and its verdict: the median steady-state AUC0-12 for each weight band must fall inside 800-1200 ugh/mL, with 1000 ugh/mL drawn as the median target, and the DAFFODIL medians did. Expressing the criterion as “within 20% of 1000 ug*h/mL” makes it exactly the tolerance the comparison table already applies.

published <- tibble::tibble(
  band = daffodil_bands,
  auclast = 1000 # ug*h/mL, midpoint of the 800-1200 target window
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = published,
  by = "band",
  params = "auclast",
  units = c(auclast = "ug*h/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated median steady-state AUC0-12 for the DAFFODIL bands vs. the 1000 ug*h/mL median target (Darwish 2025b Methods). * flags a deviation beyond the 800-1200 ug*h/mL target window."
)
Simulated median steady-state AUC0-12 for the DAFFODIL bands vs. the 1000 ugh/mL median target (Darwish 2025b Methods). flags a deviation beyond the 800-1200 ug*h/mL target window.
NCA parameter band Reference Simulated % diff
AUClast (ug*h/mL) 9-12 kg 1000 866 -13.4%
AUClast (ug*h/mL) 12-20 kg 1000 903 -9.7% 
attr(cmp, "footnote")
#> NULL

The two assertions below are evaluated on the typical subject of each band rather than on the simulated cohort median. With a single log-normal eta on each of CL and F1 and an AUC that is monotone in both, the population median AUC is the typical-value AUC, so nothing is lost – and the check becomes exactly reproducible instead of depending on which cohort the RNG happened to draw. That matters for the second claim in particular: the band ordering is a genuine but small (~2%) effect that a 100-subject cohort median cannot resolve, because the Monte Carlo standard error of the median at ~24% combined CV is itself ~3%.

band_typical <- bands |>
  dplyr::filter(study == "DAFFODIL") |>
  dplyr::mutate(
    wt_median = (wt_lo + wt_hi) / 2,
    cl = cl_typical(wt_median, gfr_mean),
    auc_typical = exp(th[["lfdepot"]]) * amt_mg / cl
  )

knitr::kable(
  band_typical[, c("band", "dose_g", "wt_median", "gfr_mean", "cl", "auc_typical")] |>
    dplyr::rename("Weight band" = band, "Dose (g)" = dose_g,
                  "Median WT (kg)" = wt_median, "GFR" = gfr_mean,
                  "Typical CL (L/h)" = cl, "Typical AUC0-12 (ug*h/mL)" = auc_typical),
  digits = c(0, 0, 1, 0, 2, 0),
  caption = "Typical-value steady-state exposure for the two DAFFODIL weight bands."
)
Typical-value steady-state exposure for the two DAFFODIL weight bands.
Weight band Dose (g) Median WT (kg) GFR Typical CL (L/h) Typical AUC0-12 (ug*h/mL)
9-12 kg 5 10.5 150 4.59 907
12-20 kg 6 16.0 150 5.63 887 

# Darwish 2025b Results / Conclusion: "Median steady-state AUC0-12 values for
# both body-weight bands fell within the target exposure range."
stopifnot(nrow(band_typical) == 2L)
stopifnot(all(band_typical$auc_typical >= 800), all(band_typical$auc_typical <= 1200))

# Darwish 2025b Discussion: "participants in the lower body weight group
# (>= 9 to < 12 kg) having slightly higher exposures compared with those in the
# higher body weight group". Structural under this model: the 5 g -> 6 g dose
# step (x1.20) is smaller than the clearance step implied by the 0.486 weight
# exponent over the same band midpoints ((16/10.5)^0.486 = x1.23).
lo <- band_typical$auc_typical[band_typical$band == "9-12 kg"]
hi <- band_typical$auc_typical[band_typical$band == "12-20 kg"]
stopifnot(lo > hi)

# The simulated cohort medians must land in the same place, but are asserted
# only loosely: they carry RNG-dependent Monte Carlo error.
band_median <- nca_auc |>
  dplyr::filter(band %in% daffodil_bands) |>
  dplyr::group_by(band) |>
  dplyr::summarise(median_auc = stats::median(auc_pknca), .groups = "drop") |>
  # Join by name rather than by position: PKNCA returns groups alphabetically,
  # so "12-20 kg" precedes "9-12 kg" and a positional compare would transpose.
  dplyr::inner_join(band_typical[, c("band", "auc_typical")], by = "band")
stopifnot(nrow(band_median) == 2L)
stopifnot(all(abs(band_median$median_auc - band_median$auc_typical) < 150))

cat(sprintf(
  "Typical AUC0-12,ss: 9-12 kg = %.0f, 12-20 kg = %.0f ug*h/mL (both inside 800-1200; lower band higher, ratio %.3f).\n",
  lo, hi, lo / hi
))
#> Typical AUC0-12,ss: 9-12 kg = 907, 12-20 kg = 887 ug*h/mL (both inside 800-1200; lower band higher, ratio 1.023).

Replicating Figure 3: distribution against the target window 

auc_by_subject <- nca_auc |>
  dplyr::left_join(dplyr::distinct(bands, study, band), by = "band") |>
  dplyr::mutate(band = factor(band, levels = band_levels))

ggplot(auc_by_subject, aes(band, auc_pknca, fill = study)) +
  geom_boxplot(outlier.alpha = 0.3) +
  geom_hline(yintercept = c(800, 1200), linetype = "dashed") +
  geom_hline(yintercept = 1000, linetype = "dotted") +
  labs(
    x = "Weight-banded dosage regimen", y = "Steady-state AUC0-12 (ug*h/mL)",
    fill = NULL,
    title = "Simulated steady-state exposure by weight band",
    caption = "Replicates Figure 3b of Darwish 2025b. Dashed lines: 800-1200 ug*h/mL target range; dotted line: 1000 ug*h/mL median target."
  ) +
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 20, hjust = 1))

attainment <- auc_by_subject |>
  dplyr::group_by(study, band) |>
  dplyr::summarise(
    N = dplyr::n(),
    `Median AUC0-12` = stats::median(auc_pknca),
    Q1 = stats::quantile(auc_pknca, 0.25),
    Q3 = stats::quantile(auc_pknca, 0.75),
    `% within 800-1200` = 100 * mean(auc_pknca >= 800 & auc_pknca <= 1200),
    .groups = "drop"
  )

knitr::kable(attainment, digits = 1,
             caption = "Simulated steady-state AUC0-12 by weight band.")
Simulated steady-state AUC0-12 by weight band.
study band N Median AUC0-12 Q1 Q3 % within 800-1200
DAFFODIL 9-12 kg 100 866.4 749.1 1023.2 53
DAFFODIL 12-20 kg 100 902.9 785.2 1060.5 60
LAVENDER 12-20 kg (L) 100 846.2 725.8 1022.0 51
LAVENDER 20-35 kg 100 863.1 718.0 1010.7 46
LAVENDER 35-50 kg 100 936.7 806.6 1082.3 59
LAVENDER >50 kg 100 940.4 803.0 1132.3 60

The paper reports the corresponding attainment percentages as 67% (>= 9 to < 12 kg) and 80% (>= 12 to < 20 kg) for DAFFODIL, and a range of 57% (>= 50 kg) to 62% (>= 35 to < 50 kg) for LAVENDER. Those DAFFODIL figures come from only 13 participants – 67% and 80% are 4/6 and 4/5 – so they are reported here for context but are deliberately not asserted on: a percentage with a denominator of five or six cannot discriminate between a correct and a mildly incorrect transcription. The gate that does the work is the median-within-target check above, which is the paper’s own stated criterion.

Covariate impact on steady-state exposure

Steady-state AUC over a dosing interval depends only on F1 and CL, so the model-implied exposure ratios for every covariate are available in closed form. FXS scales the central volume alone, so it moves Cmax without touching AUC.

gmr <- tibble::tribble(
  ~Covariate,                            ~`Exposure ratio (AUC0-12)`,
  "Rett syndrome vs. healthy",           1 / (1 + th[["e_rett_cl"]]),
  "TBI vs. healthy",                     1 / (1 + th[["e_tbi_cl"]]),
  "Fragile X syndrome vs. healthy",      1,
  "Diarrhea vs. none",                   1 + th[["e_diarrhea_f"]],
  "Fed vs. fasted",                      1 + th[["e_fed_f"]],
  "18 g vs. therapeutic dose",           1 + th[["e_dose18g_f"]],
  "24 g vs. therapeutic dose",           1 + th[["e_dose24g_f"]],
  "WT 13.4 vs. 58 kg",                   (13.4 / 58)^-th[["e_wt_cl"]],
  "GFR 150 vs. 124 mL/min/1.73 m^2",     (150 / 124)^-th[["e_crcl_cl"]]
)

knitr::kable(gmr, digits = 3,
             caption = "Model-implied steady-state AUC0-12 ratios relative to the reference covariate state.")
Model-implied steady-state AUC0-12 ratios relative to the reference covariate state.
Covariate Exposure ratio (AUC0-12)
Rett syndrome vs. healthy 1.171
TBI vs. healthy 0.814
Fragile X syndrome vs. healthy 1.000
Diarrhea vs. none 0.843
Fed vs. fasted 0.867
18 g vs. therapeutic dose 0.868
24 g vs. therapeutic dose 0.716
WT 13.4 vs. 58 kg 2.038
GFR 150 vs. 124 mL/min/1.73 m^2 0.950 

# Structural sanity: the two supratherapeutic dose levels must reduce exposure,
# 24 g more than 18 g, and low body weight must raise it (the whole reason the
# regimen is weight-banded).
stopifnot(nrow(gmr) == 9L)
r <- setNames(gmr$`Exposure ratio (AUC0-12)`, gmr$Covariate)
stopifnot(r[["24 g vs. therapeutic dose"]] < r[["18 g vs. therapeutic dose"]])
stopifnot(r[["18 g vs. therapeutic dose"]] < 1)
stopifnot(r[["WT 13.4 vs. 58 kg"]] > 1.5)

The 2 kg-scale DAFFODIL participant clears trofinetide roughly half as fast as the 58 kg reference subject, which is why a 5 g dose in a 10 kg child lands in the same exposure window as a 12 g dose in a 50 kg adolescent.

Assumptions and deviations

  • Coefficients come from Table 2, not from the Table 2 footnote equations. The footnote’s “Model equations:” block reproduces the earlier 13-study model’s coefficients verbatim and contradicts the table it annotates. Table 2 is used because it is the only source that reproduces the paper’s own stated update deltas (< 1.5% on primaries; +13.6% and +30.7% on the two RTT shifts), verified in the cross-model gate above. Only the functional form of the covariate terms is taken from the footnote.
  • Age-for-weight mapping. Darwish 2025b reports marginal covariate distributions per study (Table 1) but not the joint weight/age/GFR distribution within a weight band. Ages are drawn uniformly over a plausible window per band: 2-4 years for both DAFFODIL bands (the study’s own eligibility range) and an increasing window across the LAVENDER bands. Age enters only the central volume, so this affects Cmax and the profile shape but not steady-state AUC or any assertion in this vignette.
  • GFR distribution. Simulated as Normal(150, 27) for DAFFODIL and Normal(162, 27) for LAVENDER, truncated to 80-220 mL/min/1.73 m^2, matching the per-study means and SDs in Table 1 (ACP-2566-009: 150 (26.9); Neu-2566-Rett-003: 162 (36.9)). Note these paediatric GFR values sit well above the 124 mL/min/1.73 m^2 analysis-population median, which raises CL and partly offsets the weight effect.
  • Weight within a band. Sampled uniformly across the band; the open-ended “>50 kg” LAVENDER band is capped at 70 kg. Real weights are not uniform, so the per-band medians here should be read as band-representative rather than as reproductions of the study’s realised medians.
  • Fasted, diarrhea-free, therapeutic dose. All simulated subjects have FED = 0, AE_DIARRHEA = 0, DOSE_18G = 0, and DOSE_24G = 0. DAFFODIL and LAVENDER dosing is not described as fed or fasted in Table 1 (both the “Fed” and “Fasted” counts are zero for ACP-2566-009), so the fasted reference state is used.
  • Steady state via ss = 1. The dosing records use ss = 1 with ii = 12 rather than an explicit multi-week dosing history, which would add simulation time for no additional fidelity. The DAFFODIL titration through 2 g and 4 g BID in weeks 1-3 is therefore not simulated; the paper’s exposure comparison is likewise made at the steady state of the final weight-banded dose.
  • Attainment percentages are not asserted. The published 67% / 80% DAFFODIL figures rest on 13 participants and are too coarse to gate on; see the note under the Figure 3 replication.
  • Residual error is not applied to the figures. The plots and NCA use Cc (the individual prediction) rather than a residual-error-perturbed observation, so they correspond to the paper’s model-predicted profiles rather than to observed concentrations. Both published residual magnitudes are packaged and selected by the disease indicators; the RTT cohort simulated here exercises only the disease branch.
  • Bioavailability may exceed 1 for some simulated subjects. F1 has a typical value of 0.832 with 20.2 %CV log-normal IIV, so the upper tail crosses 1. That is a property of the published parameterization, not of this implementation.
  • All parameter values come from the paper’s text, tables, and supplement. No value was digitized from a figure, obtained by correspondence, or carried from an upstream model. Weight bands, dose levels, and PK sampling schedules come from Supplementary Material Tables S1 and S2.