Skip to contents

Treosulfan (van der Stoep 2019)

Source: vignettes/articles/vanderStoep_2019_treosulfan.Rmd
vanderStoep_2019_treosulfan.Rmd

Model and source 

ui <- rxode2::rxode(readModelDb("vanderStoep_2019_treosulfan"))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line
  • Citation: van der Stoep MYEC, Zwaveling J, Bertaina A, Locatelli F, Guchelaar HJ, Lankester AC, Moes DJAR. Population pharmacokinetics of treosulfan in paediatric patients undergoing hematopoietic stem cell transplantation. Br J Clin Pharmacol. 2019;85(9):2033-2044. doi:10.1111/bcp.13995
  • Description: Two-compartment IV-infusion population PK model for treosulfan in paediatric patients undergoing haematopoietic stem cell transplantation (van der Stoep 2019). Allometric body-weight scaling normalised to a 70 kg adult, exponents fixed at 0.75 on CL and Q and at 1 on V1 and V2, combined with a sigmoid Emax maturation function of postmenstrual age applied to BOTH CL and Q (TM50 = 38 weeks, Hill = 1.22; clearance reaches 90 percent of adult values at about 4 postnatal years). Full 4x4 correlated IIV block on CL, V1, V2 and Q, plus interoccasion variability on CL across the three conditioning days. Proportional residual error.
  • Article: https://doi.org/10.1111/bcp.13995
  • Supporting information (open access, incl. the final NONMEM control stream as Supporting Information 4): https://www.ebi.ac.uk/europepmc/webservices/rest/PMC6710524/supplementaryFiles

Treosulfan is an alkylating prodrug used in myeloablative conditioning regimens prior to haematopoietic stem cell transplantation (HSCT). van der Stoep 2019 developed a two-compartment IV-infusion population PK model on a multi-institutional paediatric dataset that is substantially larger and spans a far wider age range than the three previously published paediatric treosulfan models, and is the first of them to carry an age-maturation component.

Population

The analysis pooled 91 records (84 unique patients; 7 contributed a second record from a second transplantation more than several months later) from a prospective observational study run at Leiden University Medical Center (n = 63) and IRCCS Bambino Gesu Children’s Hospital, Rome (n = 28) between June 2011 and March 2017. Median age was 4.3 years (range 0.1-18.2) and 33 of 91 (36%) were infants aged 2 years or less – the sub-population the paper set out to characterise. Median body weight was 15.6 kg (range 3.8-75.0) and median body surface area 0.7 m^2 (range 0.3-1.9); 63.7% were male. Underlying diagnoses were haemoglobinopathy (38.5%), primary immune deficiency (28.6%), haematological malignancy (18.7%), bone marrow failure (12.1%) and other (2.2%) (van der Stoep 2019 Table 1).

Patients older than 1 year received 14 g/m^2 per day as a 3 h infusion on three consecutive days (total 42 g/m^2); patients younger than 1 year received 10 or 12 g/m^2 per day. 410 serum concentrations were analysed, 4 of which (1%) were below the 6.8 mg/L HPLC-UV lower limit of quantification and were retained at their measured values. Observed treosulfan AUC(0-inf) was 1658 mg*h/L (range 643-3371).

The same information is available programmatically via the model’s population metadata:

str(ui$population, max.level = 1)
#> List of 13
#>  $ species       : chr "human"
#>  $ n_subjects    : int 91
#>  $ n_studies     : int 1
#>  $ age_range     : chr "0.1-18.2 years (median 4.3); 33 of 91 (36%) were infants aged 2 years or less"
#>  $ weight_range  : chr "3.8-75.0 kg (median 15.6)"
#>  $ bsa_range     : chr "0.3-1.9 m^2 (median 0.7)"
#>  $ sex_female_pct: num 36.3
#>  $ disease_state : chr "Paediatric patients receiving treosulfan-based conditioning prior to allogeneic haematopoietic stem cell transp"| __truncated__
#>  $ dose_range    : chr "Intravenous treosulfan over a 3 h infusion on 3 consecutive days. Patients older than 1 year received 14 g/m^2 "| __truncated__
#>  $ regions       : chr "The Netherlands (Leiden University Medical Center, n = 63) and Italy (IRCCS Bambino Gesu Children's Hospital, Rome, n = 28)"
#>  $ renal_function: chr "eGFR by the revised Schwartz formula, capped at 120 mL/min/1.73 m^2: median 111 (range 16-120). Only 5 patients"| __truncated__
#>  $ n_observations: chr "410 serum treosulfan concentrations across 91 patients (a subset of 35 full profiles from 28 patients was used "| __truncated__
#>  $ notes         : chr "Prospective, observational, multicentre study of patients treated between June 2011 and March 2017. Seven patie"| __truncated__

Source trace

Every ini() entry in inst/modeldb/specificDrugs/vanderStoep_2019_treosulfan.R carries an in-file comment naming its origin. They are collected here for review. “Data S4” is the final NONMEM control stream published as Supporting Information 4; its $THETA block reproduces the Table 2 “Final model” column exactly and its $OMEGA / $SIGMA diagonals reproduce the Table 2 variability percentages as sqrt(variance), which is what identifies it as carrying final rather than initial estimates.

Equation / parameter Value Source location
lcl (CL at 70 kg) 18.8 L/h Table 2 “Final model”; Data S4 $THETA (0, 18.8) TH_CL
lvc (V1 at 70 kg) 20.2 L Table 2 “Final model”; Data S4 $THETA (0, 20.2) TH_V
lq (Q at 70 kg) 21.3 L/h Table 2 “Final model”; Data S4 $THETA (0, 21.3) TH_Q
lvp (V2 at 70 kg) 16.8 L Table 2 “Final model”; Data S4 $THETA (0, 16.8) TH_V2
e_wt_cl, e_wt_q 0.75 (fixed) Methods 2.5 under Eq. 4; Data S4 FSIZE=(WT/70)**0.75
e_wt_vc, e_wt_vp 1 (fixed) Methods 2.5 under Eq. 4; Data S4 (WT/70)
pma_tm50 38 weeks Table 2 “Final model”; Data S4 $THETA (0, 38) TM50
pma_hill 1.22 Data S4 $THETA (0, 1.22) Hill (Table 2 prints 1.2)
IIV block (4x4) see below Data S4 $OMEGA BLOCK(4); diagonals = Table 2 IIV 31.8 / 45.9 / 17.3 / 41.4 %CV
IOV on CL (3 occasions) 0.0194 Data S4 $OMEGA BLOCK(1) + 2x SAME; Table 2 IOV Cl 13.9%
propSd 0.1233 Data S4 $SIGMA 0.0152; Table 2 sigma 12.3%
F_size = (BW/70)^alpha n/a Equation 4
F_mat = 1/(1+(PMA/TM50)^-Hill) n/a Equation 5; Data S4 FMAT=1/(1+(PMA/TM50)**(-HILL))
Cl_tot = Cl_pop * F_size * F_mat, same form for Q n/a Equation 6 + “A similar model was used for Q”; Results 3.3
d/dt(central), d/dt(peripheral1) n/a Data S4 $DES (K=CL/V, K12=Q/V, K21=Q/V2)
Cc ~ prop(propSd) n/a Equation 3 (Y = YPRED * (1 + eps)); Data S4 $ERROR Y = F*(1+ERR(1))
PMA = 40 weeks + postnatal age n/a Methods 2.5
AUC(0-inf) = DOSE * F1 / CL, F = 1 n/a Methods 2.7 and 2.9

The Table 2 variability percentages are sqrt(omega^2), not the log-normal sqrt(exp(omega^2) - 1):

omega_diag <- c(CL = 0.101, V1 = 0.211, V2 = 0.0299, Q = 0.171)
published_cv <- c(CL = 31.8, V1 = 45.9, V2 = 17.3, Q = 41.4)
data.frame(
  parameter = names(omega_diag),
  `omega^2` = unname(omega_diag),
  `sqrt(omega^2) %` = round(100 * sqrt(omega_diag), 1),
  `published Table 2 %` = unname(published_cv),
  check.names = FALSE
)
#>    parameter omega^2 sqrt(omega^2) % published Table 2 %
#> CL        CL  0.1010            31.8                31.8
#> V1        V1  0.2110            45.9                45.9
#> V2        V2  0.0299            17.3                17.3
#> Q          Q  0.1710            41.4                41.4
stopifnot(all(abs(100 * sqrt(omega_diag) - published_cv) < 0.1))

The 4x4 block is positive definite, but only just – a direct consequence of the 70 kg standardisation, which the Discussion notes inflates the V2 and Q IIV relative standard errors to 201% and 153% (against 53% and 46% when standardising to the cohort median 15.6 kg). The published model kept 70 kg for comparability with other publications, so that is what the model file encodes.

omega <- matrix(c(
   0.101,   0.131,  0.0349, -0.0307,
   0.131,   0.211,  0.026,  -0.0354,
   0.0349,  0.026,  0.0299, -0.05,
  -0.0307, -0.0354, -0.05,   0.171
), nrow = 4, byrow = TRUE,
dimnames = list(c("CL", "V1", "V2", "Q"), c("CL", "V1", "V2", "Q")))

round(cov2cor(omega), 3)
#>        CL     V1     V2      Q
#> CL  1.000  0.897  0.635 -0.234
#> V1  0.897  1.000  0.327 -0.186
#> V2  0.635  0.327  1.000 -0.699
#> Q  -0.234 -0.186 -0.699  1.000
ev_omega <- eigen(omega, only.values = TRUE)$values
stopifnot(min(ev_omega) > 0)   # positive definite: rxode2's Cholesky sampler works
signif(min(ev_omega), 3)
#> [1] 8.4e-05

Maturation function (replicates Figure 2)

Equation 5 is a sigmoid Emax function of postmenstrual age driving both CL and Q. Figure 2 of the paper plots it as a percentage of the mature (adult) value against PMA in weeks, with dashed landmarks at birth, 1 year and 4 years.

tm50 <- 38
hill <- 1.22
fmat <- function(pma_weeks) 1 / (1 + (pma_weeks / tm50)^(-hill))
# The paper defines PMA as 40 weeks of assumed gestation plus postnatal age.
pma_of_age <- function(age_years) 40 + age_years * 52.1775

mat_curve <- tibble(pma = seq(31, 390, by = 1)) |>
  mutate(pct = 100 * fmat(pma))

landmarks <- tibble(
  label = c("Birth", "1 year", "4 years"),
  pma = pma_of_age(c(0, 1, 4))
) |>
  mutate(pct = 100 * fmat(pma))

ggplot(mat_curve, aes(pma, pct)) +
  geom_line(colour = "steelblue", linewidth = 1) +
  geom_vline(data = landmarks, aes(xintercept = pma),
             colour = "red", linetype = "dashed") +
  geom_text(data = landmarks, aes(x = pma, y = 41, label = label),
            hjust = -0.08, fontface = "italic", size = 3) +
  scale_y_continuous(limits = c(40, 100)) +
  labs(x = "Postmenstrual age (weeks)", y = "% of matured clearance",
       caption = "Replicates Figure 2 of van der Stoep 2019.")


landmarks |> mutate(pct = round(pct, 1))
#> # A tibble: 3 × 3
#>   label     pma   pct
#>   <chr>   <dbl> <dbl>
#> 1 Birth    40    51.6
#> 2 1 year   92.2  74.7
#> 3 4 years 249.   90.8

Four claims the paper makes about this curve, checked against the encoded parameters. All are deterministic – no simulation, no random draws – so the tolerances are tight.

# 1. "The maturation of treosulfan clearance reaches 50% of adult values at
#    38 weeks PMA" (Results 3.3). True by construction of TM50.
stopifnot(abs(100 * fmat(38) - 50) < 1e-8)

# 2. "that is 2 weeks prior to birth assuming a full-term gestational age of
#    40 weeks" -- so the curve crosses just above 50% at birth.
stopifnot(100 * fmat(pma_of_age(0)) > 50,
          100 * fmat(pma_of_age(0)) < 55)

# 3. "Clearance reaches 90% of adult values at age approximately 4 years"
#    (Results 3.3, Discussion, and the "What this study adds" box).
stopifnot(abs(100 * fmat(pma_of_age(4)) - 90) < 3)

# 4. The Figure 2 panel runs to ~390 weeks PMA and plateaus near 94-95%,
#    not at 100% -- a property of the Hill coefficient, so this discriminates
#    Hill = 1.22 from a much steeper or shallower value.
stopifnot(abs(100 * fmat(390) - 94.5) < 1.5)

Virtual cohort: the paper’s own dosing simulation (Table 3)

Methods 2.7 describes the paper’s simulation: individual treosulfan doses were computed to target an AUC(0-inf) of 1650 mg*h/L – the daily median exposure in patients receiving the most common 14 g/m^2 dose – for virtual patients whose weights are the 5th, 50th and 95th CDC weight-for-age percentiles for boys. Table 3 publishes the resulting 66 (age, weight, dose) triplets. Because AUC(0-inf) = DOSE / CL with F = 1, reproducing Table 3 is an exact, fully deterministic test of the typical-value clearance model: the population CL, the 0.75 allometric exponent, the maturation function, and the PMA = 40 weeks + postnatal age convention all have to be right simultaneously.

# van der Stoep 2019 Table 3, transcribed verbatim (age category, CDC
# weight-for-age percentile weight in kg, recommended daily dose in mg).
table3 <- tibble::tribble(
  ~age_y, ~wt,  ~dose_pub,
  0,      2.6,  1350,   0,     3.3,  1600,   0,     4.2,  1950,
  0.25,   5.2,  2650,   0.25,  6.4,  3100,   0.25,  7.7,  3550,
  0.5,    6.6,  3500,   0.5,   7.9,  4000,   0.5,   9.5,  4600,
  0.75,   7.4,  4100,   0.75,  8.9,  4700,   0.75, 10.6,  5300,
  1,      8.1,  4600,   1,     9.6,  5200,   1,    11.5,  6000,
  2,     10.1,  6100,   2,    12.2,  7000,   2,    14.7,  8000,
  3,     12.0,  7300,   3,    14.3,  8300,   3,    17.3,  9600,
  4,     13.6,  8250,   4,    16.3,  9450,   4,    20.3, 11100,
  5,     15.2,  9100,   5,    18.5, 10500,   5,    23.5, 12700,
  6,     16.9, 10000,   6,    20.8, 11700,   6,    27.0, 14250,
  7,     18.7, 11000,   7,    23.2, 12800,   7,    30.9, 16000,
  8,     20.7, 11900,   8,    25.8, 14000,   8,    35.3, 17700,
  9,     22.7, 12800,   9,    28.7, 15250,   9,    40.4, 19700,
  10,    24.9, 13800,  10,    32.1, 16700,  10,    46.2, 21900,
  11,    27.5, 14900,  11,    36.1, 18300,  11,    52.6, 24200,
  12,    30.6, 16200,  12,    40.7, 20000,  12,    59.3, 26600,
  13,    34.2, 17600,  13,    45.8, 22000,  13,    66.1, 28900,
  14,    38.5, 19300,  14,    51.2, 24000,  14,    72.7, 31100,
  15,    43.0, 21000,  15,    56.5, 25800,  15,    78.8, 33100,
  16,    47.3, 22600,  16,    61.1, 27400,  16,    84.3, 34900,
  17,    50.8, 23900,  17,    64.7, 28700,  17,    88.8, 36300,
  18,    53.2, 24800,  18,    67.3, 29500,  18,    92.0, 37400
) |>
  mutate(
    id = row_number(),
    PAGE = pma_of_age(age_y),
    WT = wt,
    OCC = 1L,
    arm = "Table 3 target (typical values)"
  )

stopifnot(nrow(table3) == 66L, !anyDuplicated(table3$id))
# Build the event table by hand so the covariate columns ride along. Doses are
# 3 h infusions into `central`; observation rows also name `central` (the ODE
# state), never `Cc` (the algebraic observable).
obs_grid <- seq(0, 24, by = 0.05)

make_events <- function(subj, grid) {
  dosing <- subj |>
    transmute(id, time = 0, amt = dose_mg, evid = 1L, dur = 3,
              cmt = "central", WT, PAGE, OCC, arm)
  obs <- subj |>
    select(id, WT, PAGE, OCC, arm) |>
    tidyr::crossing(time = grid) |>
    mutate(amt = NA_real_, evid = 0L, dur = NA_real_, cmt = "central")
  bind_rows(dosing, obs) |>
    arrange(id, time, desc(evid))
}

ev_t3 <- table3 |>
  mutate(dose_mg = dose_pub) |>
  make_events(obs_grid)

stopifnot(!anyDuplicated(unique(ev_t3[, c("id", "time", "evid")])))
mod <- readModelDb("vanderStoep_2019_treosulfan")

sim_t3 <- rxode2::rxSolve(
  rxode2::zeroRe(mod),
  events = ev_t3,
  keep = c("WT", "PAGE", "arm")
) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalvp', 'etalq', 'etaiov_cl_1', 'etaiov_cl_2', 'etaiov_cl_3'
#> Warning: multi-subject simulation without without 'omega'

stopifnot(all(sim_t3$Cc >= 0, na.rm = TRUE))

Table 3 exposure gate 

sim_nca_t3 <- sim_t3 |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, arm)

# Guarantee a time-zero record per subject (IV infusion: pre-dose Cc = 0).
sim_nca_t3 <- bind_rows(
  sim_nca_t3,
  sim_nca_t3 |> distinct(id, arm) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, arm, time, .keep_all = TRUE) |>
  arrange(id, arm, time)

conc_t3 <- PKNCA::PKNCAconc(sim_nca_t3, Cc ~ time | arm + id)

dose_t3 <- ev_t3 |>
  filter(evid == 1L) |>
  select(id, time, amt, arm)
dose_obj_t3 <- PKNCA::PKNCAdose(dose_t3, amt ~ time | arm + id)

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, auclast = TRUE,
  aucinf.obs = TRUE, half.life = TRUE
)

nca_t3 <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_t3, dose_obj_t3, intervals = intervals)
)
auc_t3 <- as.data.frame(nca_t3) |>
  filter(PPTESTCD == "aucinf.obs") |>
  select(id, aucinf = PPORRES) |>
  inner_join(table3, by = "id") |>
  mutate(
    # Closed form from Methods 2.7 / 2.9: AUC(0-inf) = DOSE * F1 / CL, F = 1.
    cl_closed = 18.8 * (wt / 70)^0.75 * fmat(PAGE),
    auc_closed = dose_pub / cl_closed,
    pct_vs_closed = 100 * (aucinf - auc_closed) / auc_closed,
    pct_vs_target = 100 * (aucinf - 1650) / 1650
  )

summary(auc_t3$pct_vs_closed)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -0.004693 -0.004206 -0.003874 -0.003853 -0.003492 -0.002964
summary(auc_t3$pct_vs_target)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -1.132322 -0.114930 -0.001419 -0.012471  0.131662  0.743695

Two things are being checked here.

First, the NCA of the solved profile must recover the closed form DOSE / CL exactly. Both sides use the same parameter draw (there is none – zeroRe() zeroes every random effect), so the only difference is trapezoidal / extrapolation error on a finely-sampled grid; a tight bound is the correct assertion, and it is the one that catches a structural transcription error.

stopifnot(max(abs(auc_t3$pct_vs_closed)) < 0.5)

Second, each of the 66 published doses must land on the paper’s own 1650 mg*h/L target. The residual scatter here is the paper’s rounding of Table 3 to the nearest 50 mg, which is largest in absolute terms for the smallest neonatal doses (1350-1950 mg).

stopifnot(
  # Centre: a mis-transcribed CL, exponent or maturation term moves the whole
  # distribution by tens of percent.
  abs(median(auc_t3$pct_vs_target)) < 1,
  # Envelope: the worst single row is the 3.3 kg neonate, where 50 mg of dose
  # rounding is 3% of the dose. Realised max 1.14%; the bound leaves headroom
  # for the grid and for a re-transcription, and still goes red on any real
  # structural error.
  max(abs(auc_t3$pct_vs_target)) < 3
)
auc_t3 |>
  filter(age_y %in% c(0, 1, 4, 10, 18)) |>
  transmute(
    "Age" = ifelse(age_y == 0, "0 mo", paste0(age_y, " y")),
    "Weight (kg)" = wt,
    "Published dose (mg/day)" = dose_pub,
    "Model CL (L/h)" = round(cl_closed, 3),
    "Simulated AUC0-inf (mg*h/L)" = round(aucinf, 0),
    "% vs 1650 target" = round(pct_vs_target, 2)
  ) |>
  knitr::kable(
    caption = "Table 3 dose recommendations re-derived from the packaged model (selected rows; all 66 are gated above)."
  )
Table 3 dose recommendations re-derived from the packaged model (selected rows; all 66 are gated above).
Age Weight (kg) Published dose (mg/day) Model CL (L/h) Simulated AUC0-inf (mg*h/L) % vs 1650 target
0 mo 2.6 1350 0.820 1646 -0.25
0 mo 3.3 1600 0.981 1631 -1.13
0 mo 4.2 1950 1.175 1659 0.56
1 y 8.1 4600 2.785 1652 0.09
1 y 9.6 5200 3.164 1644 -0.39
1 y 11.5 6000 3.622 1656 0.38
4 y 13.6 8250 4.997 1651 0.06
4 y 16.3 9450 5.724 1651 0.06
4 y 20.3 11100 6.748 1645 -0.30
10 y 24.9 13800 8.347 1653 0.19
10 y 32.1 16700 10.099 1654 0.22
10 y 46.2 21900 13.270 1650 0.02
18 y 53.2 24800 15.018 1651 0.08
18 y 67.3 29500 17.913 1647 -0.20
18 y 92.0 37400 22.647 1651 0.08

Dose-vs-age relationship (replicates Figure 5)

Figure 5 plots the daily dose required to reach the 1650 mg*h/L target against age, panel A in mg/kg and panel B as the absolute dose, with the 50th weight percentile as a solid line and the 5th and 95th as dashed. The paper’s reading of it: “The recommended treosulfan dose per kg is lower in early years of life and reaches a maximum at approximately 4 years accounting for maturation of clearance and because dose per kg is higher in younger children based on allometric theory.”

fig5 <- auc_t3 |>
  arrange(age_y, wt) |>
  group_by(age_y) |>
  mutate(percentile = c("5th", "50th", "95th")[rank(wt)]) |>
  ungroup() |>
  mutate(dose_per_kg = dose_pub / wt)

p_a <- ggplot(fig5, aes(age_y, dose_per_kg, group = percentile,
                        linetype = percentile)) +
  geom_line() +
  scale_linetype_manual(values = c("5th" = "dashed", "50th" = "solid",
                                   "95th" = "dashed")) +
  labs(x = "Age (years)", y = "Treosulfan dose (mg/kg per day)",
       title = "A", linetype = "Weight percentile")

p_b <- ggplot(fig5, aes(age_y, dose_pub, group = percentile,
                        linetype = percentile)) +
  geom_line() +
  scale_linetype_manual(values = c("5th" = "dashed", "50th" = "solid",
                                   "95th" = "dashed")) +
  labs(x = "Age (years)", y = "Treosulfan dose (mg per day)",
       title = "B", linetype = "Weight percentile")

p_a

p_b

# The paper's claim: dose per kg peaks at about 4 years on the 50th-percentile
# curve. Deterministic (Table 3 is a printed table), so this is exact.
peak_age <- fig5 |>
  filter(percentile == "50th") |>
  slice_max(dose_per_kg, n = 1) |>
  pull(age_y)
peak_age
#> [1] 3
stopifnot(peak_age >= 3, peak_age <= 5)

# Panel B: absolute dose rises monotonically with age on every percentile curve.
stopifnot(
  fig5 |>
    arrange(percentile, age_y) |>
    group_by(percentile) |>
    summarise(mono = all(diff(dose_pub) > 0), .groups = "drop") |>
    pull(mono) |>
    all()
)

Study-protocol cohort with between-subject variability (replicates Figure 1)

Figure 1 shows the observed full concentration-time profiles of the 27 patients sampled at 1.5, 3.5, 4, 5, 7 and 9 h after the start of a 3 h infusion of 14 g/m^2. Reproducing it, and the published cohort exposure, needs three things the paper reports only as marginal summaries: an age distribution, a weight distribution, and a dose in mg (i.e. a body surface area per subject).

The cohort below is built to reproduce the published marginals exactly where they are stated:

  • Age is drawn from a piecewise-linear quantile function pinned to the three published facts – minimum 0.1 y, 36% aged 2 y or less, median 4.3 y, maximum 18.2 y (Table 1).
  • Weight is the Table 3 50th-percentile CDC weight-for-age at that age, multiplied by 0.88 and by log-normal scatter, which brings the cohort weight median onto the published 15.6 kg. The sub-unit factor is physiologically expected: these are children sick enough to need a transplant, so they sit below the healthy weight-for-age median.
  • BSA is scaled from weight as BSA = 0.7 * (WT / 15.6)^(2/3), anchored on the Table 1 medians (15.6 kg, 0.7 m^2) using the standard two-thirds-power geometric relationship; BSA itself was computed by the Mosteller formula in the paper (Supporting Information 2), which needs a height the paper does not publish. See “Assumptions and deviations” for the check that this reproduces the published BSA range from the published weight range.
  • Dose follows the study protocol, not a flat 14 g/m^2. Methods 2.1: patients older than 1 year received 14 g/m^2 per day, while patients younger than 1 year received 10 or 12 g/m^2 (16 and 2 of 91 patients respectively). This matters: infants are exactly the subgroup whose maturation fraction is lowest, so dosing them at the older-child rate would inflate their exposure by about 40% and bias the whole cohort median upward. Reproducing the published AUC requires reproducing the published dose rule.
# `set.seed()` seeds R's RNG, not rxode2's; rxode2 partitions its streams per
# solver thread, so the cohort below differs between a 2-thread CI runner and a
# 16-thread workstation. Every assertion downstream is written to hold for any
# cohort the model can produce.
set.seed(20190901)

n_cohort <- 150L # <= 200 per arm

# Quantile function pinned to the Table 1 age marginals.
age_quantile <- approxfun(
  x = c(0, 0.36, 0.50, 1.00),
  y = c(0.1, 2.0, 4.3, 18.2)
)

wt50 <- fig5 |> filter(percentile == "50th") |> arrange(age_y)

cohort14 <- tibble(
  id = 1000L + seq_len(n_cohort),
  age_y = age_quantile(runif(n_cohort))
) |>
  mutate(
    WT = pmin(
      pmax(
        approx(wt50$age_y, wt50$wt, xout = age_y, rule = 2)$y *
          0.88 * exp(rnorm(n_cohort, 0, 0.25)),
        3.8
      ),
      75.0
    ),
    PAGE = pma_of_age(age_y),
    OCC = 1L,
    bsa = 0.7 * (WT / 15.6)^(2 / 3),
    # Study protocol (Methods 2.1): 10 g/m^2 under 1 year, 14 g/m^2 at or above.
    dose_g_m2 = ifelse(age_y < 1, 10, 14),
    dose_mg = dose_g_m2 * 1000 * bsa,
    arm = "Study protocol cohort (with IIV)"
  )

summary(cohort14[, c("age_y", "WT", "bsa", "dose_mg")])
#>      age_y               WT              bsa            dose_mg     
#>  Min.   : 0.1016   Min.   : 3.800   Min.   :0.2730   Min.   : 2730  
#>  1st Qu.: 1.7666   1st Qu.: 9.676   1st Qu.:0.5091   1st Qu.: 6943  
#>  Median : 4.8205   Median :16.995   Median :0.7411   Median :10376  
#>  Mean   : 6.8153   Mean   :23.136   Mean   :0.8576   Mean   :11789  
#>  3rd Qu.:11.6387   3rd Qu.:33.303   3rd Qu.:1.1606   3rd Qu.:16248  
#>  Max.   :17.9953   Max.   :75.000   Max.   :1.9940   Max.   :27916

The reconstructed cohort against the published marginals (Table 1):

tibble::tribble(
  ~Quantity, ~Published, ~Simulated,
  "Median age (y)", "4.3", sprintf("%.1f", median(cohort14$age_y)),
  "Age range (y)", "0.1-18.2",
  sprintf("%.1f-%.1f", min(cohort14$age_y), max(cohort14$age_y)),
  "Aged 2 y or less", "36%",
  sprintf("%.0f%%", 100 * mean(cohort14$age_y <= 2)),
  "Median weight (kg)", "15.6", sprintf("%.1f", median(cohort14$WT)),
  "Weight range (kg)", "3.8-75.0",
  sprintf("%.1f-%.1f", min(cohort14$WT), max(cohort14$WT)),
  "Median BSA (m^2)", "0.7", sprintf("%.1f", median(cohort14$bsa)),
  "Receiving 14 g/m^2", "80%",
  sprintf("%.0f%%", 100 * mean(cohort14$dose_g_m2 == 14))
) |>
  knitr::kable(caption = "Reconstructed cohort vs van der Stoep 2019 Table 1.")
Reconstructed cohort vs van der Stoep 2019 Table 1.
Quantity Published Simulated
Median age (y) 4.3 4.8
Age range (y) 0.1-18.2 0.1-18.0
Aged 2 y or less 36% 32%
Median weight (kg) 15.6 17.0
Weight range (kg) 3.8-75.0 3.8-75.0
Median BSA (m^2) 0.7 0.7
Receiving 14 g/m^2 80% 86% 
ev_14 <- make_events(cohort14, obs_grid)
stopifnot(!anyDuplicated(unique(ev_14[, c("id", "time", "evid")])))

rxode2::rxSetSeed(20190901)
sim_14 <- rxode2::rxSolve(
  mod,
  events = ev_14,
  keep = c("WT", "PAGE", "arm")
) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line

stopifnot(all(sim_14$Cc >= 0, na.rm = TRUE))
# Replicates Figure 1 of van der Stoep 2019: full concentration-time profiles at
# 14 g/m^2, shown on the paper's sampling times and log concentration axis. The
# published panel holds 27 patients, all of whom received 14 g/m^2, so the
# subset is drawn from the 14 g/m^2 subjects only.
sample_times <- c(1.5, 3.5, 4, 5, 7, 9)

ids_14 <- cohort14$id[cohort14$dose_g_m2 == 14]

fig1 <- sim_14 |>
  filter(time %in% sample_times) |>
  filter(id %in% sample(ids_14, 27))

ggplot(fig1, aes(time, Cc, group = id)) +
  geom_line(linetype = "dotted", linewidth = 0.3) +
  geom_point(size = 1.2) +
  scale_y_log10(limits = c(1, 1000),
                breaks = c(1, 10, 100, 1000),
                labels = c("1", "10", "100", "1000")) +
  scale_x_continuous(limits = c(0, 14), breaks = seq(0, 14, by = 2)) +
  labs(x = "Time (h)", y = "Treosulfan blood concentration (mg/L)",
       caption = "Replicates Figure 1 of van der Stoep 2019 (27 simulated subjects at 14 g/m2).")

The published panel spans roughly 190-620 mg/L at 1.5 h and falls to roughly 5-80 mg/L at 9 h. The simulated envelope is compared on the same two anchor times. These are stochastic quantities, so the assertions are on robust quantiles and are set outside the range seen across repeated renders.

envelope <- sim_14 |>
  filter(time %in% c(1.5, 9), id %in% ids_14) |>
  group_by(time) |>
  summarise(
    q10 = quantile(Cc, 0.10),
    q50 = quantile(Cc, 0.50),
    q90 = quantile(Cc, 0.90),
    .groups = "drop"
  )
envelope
#> # A tibble: 2 × 4
#>    time    q10   q50   q90
#>   <dbl>  <dbl> <dbl> <dbl>
#> 1   1.5 240.   352.  546. 
#> 2   9     7.03  17.6  39.8

# Figure 1 medians read approximately 350 mg/L at 1.5 h and 20 mg/L at 9 h. A
# mis-transcribed volume, clearance or dose unit moves these by a factor, not by
# tens of percent, so a 2-fold band is a gate that still goes red.
med_15 <- envelope$q50[envelope$time == 1.5]
med_9 <- envelope$q50[envelope$time == 9]
stopifnot(med_15 > 175, med_15 < 700)
stopifnot(med_9 > 10, med_9 < 40)

PKNCA validation 

sim_nca_14 <- sim_14 |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, arm)

sim_nca_14 <- bind_rows(
  sim_nca_14,
  sim_nca_14 |> distinct(id, arm) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, arm, time, .keep_all = TRUE) |>
  arrange(id, arm, time)

conc_14 <- PKNCA::PKNCAconc(sim_nca_14, Cc ~ time | arm + id)

dose_14 <- ev_14 |>
  filter(evid == 1L) |>
  select(id, time, amt, arm)
dose_obj_14 <- PKNCA::PKNCAdose(dose_14, amt ~ time | arm + id)

nca_14 <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_14, dose_obj_14, intervals = intervals)
)
nca_all <- bind_rows(as.data.frame(nca_t3), as.data.frame(nca_14))

nca_all |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life")) |>
  group_by(arm, PPTESTCD) |>
  summarise(
    median = median(PPORRES),
    p10 = quantile(PPORRES, 0.10),
    p90 = quantile(PPORRES, 0.90),
    .groups = "drop"
  ) |>
  mutate(across(c(median, p10, p90), \(x) signif(x, 3))) |>
  rename(
    "Arm" = arm,
    "NCA parameter" = PPTESTCD,
    "Median" = median,
    "10th pctile" = p10,
    "90th pctile" = p90
  ) |>
  knitr::kable(caption = "Simulated NCA summary by arm.")
Simulated NCA summary by arm.
Arm NCA parameter Median 10th pctile 90th pctile
Study protocol cohort (with IIV) aucinf.obs 1730.00 1140.00 2770.00
Study protocol cohort (with IIV) cmax 465.00 311.00 732.00
Study protocol cohort (with IIV) half.life 1.44 1.10 1.94
Study protocol cohort (with IIV) tmax 3.00 3.00 3.00
Table 3 target (typical values) aucinf.obs 1650.00 1650.00 1650.00
Table 3 target (typical values) cmax 448.00 425.00 465.00
Table 3 target (typical values) half.life 1.43 1.27 1.67
Table 3 target (typical values) tmax 3.00 3.00 3.00

Comparison against published exposure

van der Stoep 2019 does not publish Cmax, Tmax or half-life. It publishes two AUC(0-inf) figures, both of which this model should land on: the 1650 mgh/L target that Table 3 was constructed around (Methods 2.7), and the 1658 mgh/L median observed in the analysis cohort (Table 1).

published <- tibble::tribble(
  ~arm, ~aucinf.obs,
  "Table 3 target (typical values)", 1650,
  "Study protocol cohort (with IIV)", 1658
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_all,
  reference = published,
  by = "arm",
  units = c(aucinf.obs = "mg*h/L"),
  tolerance_pct = 25
)

knitr::kable(
  cmp,
  caption = "Simulated vs. published treosulfan AUC(0-inf). * differs from reference by >25%.",
  align = c("l", "l", "r", "r", "r")
)
Simulated vs. published treosulfan AUC(0-inf). * differs from reference by >25%.
NCA parameter arm Reference Simulated % diff
AUC0-∞ (obs) (mg*h/L) Table 3 target (typical values) 1650 1650 -0.0%
AUC0-∞ (obs) (mg*h/L) Study protocol cohort (with IIV) 1660 1730 +4.3% 
attr(cmp, "footnote")
#> NULL
pct <- as.numeric(gsub("[^0-9.+-]", "", cmp$`% diff`))
names(pct) <- cmp$arm

# The typical-value arm is deterministic -- no random effects, no reconstructed
# covariates -- and lands on the paper's own target to well under a percent. It
# carries the tight bound and is the load-bearing gate on the structural model.
stopifnot(abs(pct[["Table 3 target (typical values)"]]) < 1)

# The IIV arm is a cohort median from 150 subjects whose age, weight and BSA are
# all reconstructed, so it carries two irreducible sources of scatter. Over 400
# replicate cohorts of this design the median AUC realised 1748 +/- 61 (SD),
# i.e. -6.2% to +17.3% against the published 1658 -- roughly +5% of structural
# offset from the reconstruction, on top of +/-7% of pure sampling noise at
# n = 150. The bound is set outside that measured range and must not be
# tightened back: a transcription error in a dose, volume or clearance moves
# this by a factor, not by tens of percent, so 25% still goes red on anything
# real.
stopifnot(abs(pct[["Study protocol cohort (with IIV)"]]) < 25)

The IIV arm sits slightly above the observed 1658 mg*h/L. The residual is attributable to the reconstruction rather than to the model: the patients’ own ages, weights and body surface areas are not published, and the largest single driver of individual exposure in this model – the maturation fraction, which ranges from about 0.52 at birth to 0.98 in adolescence – is a function of exactly the covariate that has to be reconstructed. Dosing every subject at 14 g/m^2 instead of following the protocol’s reduced infant dose raises this arm’s median by a further 8% or so, which is what makes the protocol dose rule worth encoding. The typical-value arm carries no reconstruction at all and is the gate that would catch a structural error.

Assumptions and deviations

  • Hill coefficient 1.22 rather than the printed 1.2. Table 2 prints the final Hill estimate as 1.2; the Data S4 control stream carries 1.22. Every other final $THETA in that stream matches its Table 2 entry exactly, which identifies it as carrying the final estimates, so the unrounded value is used. It is also marginally the better reproduction of Table 3 (root-mean- square deviation across the 66 rows 0.40% at 1.22 against 0.49% at 1.20), and the two are indistinguishable for practical purposes.

  • Gestational age is assumed, not observed. Methods 2.5: “PMA was estimated by adding a gestational age of 40 weeks to postnatal age.” The cohort carries no recorded gestational age, so PAGE in this model is always 40 + postnatal age in weeks. A user with genuine gestational ages should supply PAGE = GA_weeks + postnatal age in weeks instead; the maturation term is poorly informed below term either way, since the youngest patient was 0.1 years old.

  • PAGE is carried in weeks, not the register-default months. The PAGE canonical in inst/references/covariate-columns.md defaults to months but explicitly provides for models whose source equations are written in weeks. TM50 = 38 weeks only makes sense on the week scale. The model file declares units = "weeks" in its covariateData$PAGE entry.

  • Body surface area is reconstructed for the IIV arm. The paper doses by BSA but publishes BSA only as a cohort median and range. The scaling BSA = 0.7 * (WT / 15.6)^(2/3) is anchored on the Table 1 medians (15.6 kg, 0.7 m^2) and uses the standard two-thirds-power geometric relationship. It reproduces the published BSA range from the published weight range to within the printed precision:

    bsa_from_wt <- function(wt) 0.7 * (wt / 15.6)^(2 / 3)
c(min_wt_3.8kg = bsa_from_wt(3.8), max_wt_75kg = bsa_from_wt(75))
#> min_wt_3.8kg  max_wt_75kg 
#>     0.273026     1.993988
# Table 1 reports BSA 0.7 m^2 (range 0.3-1.9).
stopifnot(
  abs(bsa_from_wt(3.8) - 0.3) < 0.1,
  abs(bsa_from_wt(75) - 1.9) < 0.15
)

    The Table 3 arm – which carries every gate on the structural model – uses the paper’s own published doses in mg and needs no BSA at all.

  • Age, weight and their pairing in the IIV arm are reconstructed. Ages come from a piecewise-linear quantile function pinned to the four published age facts (minimum, 36th percentile, median, maximum); weights are the Table 3 50th-percentile CDC weight-for-age curve scaled by 0.88 with log-normal scatter, the sub-unit factor bringing the cohort weight median onto the published 15.6 kg and reflecting that transplant candidates sit below the healthy weight-for-age median. The interpolation between the published anchors is an assumption; no per-patient covariate data are published. This affects the IIV arm’s exposure distribution only – no model parameter depends on it.

  • The IIV arm follows the study’s protocol dose, not a flat 14 g/m^2. Methods 2.1 gives 14 g/m^2 per day to patients older than 1 year and 10 or 12 g/m^2 to those younger (16 and 2 of 91 patients); the vignette uses 10 g/m^2 below 1 year, the dominant of the two. This is load-bearing for the exposure comparison rather than cosmetic: infants have the lowest maturation fraction, so dosing them at the older-child rate inflates their AUC by about 40% and the cohort median by about 8%.

  • eGFR is deliberately absent. Results 3.3 reports that eGFR was a statistically significant covariate on CL in the stepwise screen (dOFV = -16.72), but the prediction-corrected VPC worsened and the IIV of the PK parameters increased, so the authors excluded it from the final model. The model file therefore has no renal covariate; the screen is recorded in population$renal_function.

  • Interoccasion variability is encoded for three occasions. Methods 2.4 defines each dose and its subsequent sampling as a separate occasion and the regimen is three consecutive days; Data S4 multiplexes three etas off the DAY column with a shared variance. Simulations here use a single day, so OCC = 1 throughout and only the first IOV eta is active. A user simulating the full three-day conditioning course should supply OCC = 1, 2, 3 on the respective days.

  • Residual error is proportional only. Methods 2.4 states that a combined proportional-plus-additive model was examined and the proportional model was retained. Four of 410 samples fell below the 6.8 mg/L LLOQ and were kept at their measured values (the Hecht 2011 approach), so no censoring is encoded.

  • NCA is run on uncensored simulated profiles. The paper’s AUC(0-inf) is model-derived (DOSE * F1 / CL, Methods 2.9), not NCA-derived, so an uncensored NCA is the correct comparator for it. Half-life is reported in the summary table for information only and is not gated: lambda-z estimated from uncensored simulated concentrations that decay far below the 6.8 mg/L assay limit is biased relative to any half-life an analyst would obtain from real assay data, and the paper publishes no half-life to compare against.

  • No non-paper-derived parameter values. Every ini() entry comes from van der Stoep 2019 Table 2 or from the Supporting Information 4 NONMEM control stream. No value was digitised from a figure, supplied by correspondence, or carried from an upstream model.