Skip to contents

Midazolam paediatric scaling (Cella 2012)

Source: vignettes/articles/Cella_2012_midazolam_paediatric_scaling.Rmd
Cella_2012_midazolam_paediatric_scaling.Rmd

Model and source

Cella et al. (2012) developed two independent population PK models for midazolam on different paediatric subpopulations as part of a methodological study of allometric scaling across age groups. Both models are packaged in nlmixr2lib with one shared vignette:

  • Cella_2012_midazolam_infants_adults – two-compartment model with first-order absorption fit to 23 infants and toddlers plus 34 adults.
  • Cella_2012_midazolam_children_adolescents – two-compartment IV-bolus model fit to 18 children and adolescents under sparse sampling with informative priors from De Wildt 2002.
m1_meta <- rxode2::rxode2(readModelDb("Cella_2012_midazolam_infants_adults"))
#> ℹ parameter labels from comments will be replaced by 'label()'
m2_meta <- rxode2::rxode2(readModelDb("Cella_2012_midazolam_children_adolescents"))
#> ℹ parameter labels from comments will be replaced by 'label()'
cat(strwrap(m1_meta$description, width = 78), sep = "\n")
#> Two-compartment population PK model for midazolam in infants, toddlers, and
#> adults (Cella 2012 Model 1), with first-order absorption supporting
#> intravenous and oral dosing, body-weight allometric scaling of clearance
#> (exponent fixed to 0.75 at a 70 kg reference), per-kg linear scaling of the
#> central volume, and a constant peripheral volume. Pooled cohort of 23 infants
#> and toddlers in a paediatric surgical ICU and 34 healthy adult volunteers.

The model functions and their population metadata are also accessible programmatically via readModelDb("Cella_2012_midazolam_infants_adults")$population and readModelDb("Cella_2012_midazolam_children_adolescents")$population after the model is loaded.

Population

Model 1 (infants + toddlers + adults, n = 57). The combined dataset pools three healthy-adult crossover studies from the Centre for Human Drug Research (89110-pilot, 89110, 94113; 34 adults, 19.9-29.7 years, 59-91 kg) with a single infant / toddler study at the paediatric surgical ICU of Erasmus MC - Sophia’s Children’s Hospital (23 infants and toddlers, 3.2-24.7 months, 5.1-12 kg). Infants received 0.1 mg/kg IV bolus + 0.05 mg/kg/h IV infusion; adults received 0.1-0.15 mg/kg IV (bolus or short infusion) or 5/7.5/10 mg fixed oral doses banded by body-weight tier (Cella 2012 Table 1).

Model 2 (children + adolescents, n = 18). Paediatric oncology patients (3.2-16.2 years, 12.6-60.1 kg) from a collaborative study between Purdue University and Sophia Children’s Hospital. All received a single IV bolus of midazolam at ~ 0.12 mg/kg before invasive procedures with sparse sampling (mean 4.6 samples per subject). The published model was fit with informative priors from De Wildt 2002 (CL = 5 mL/kg/min, Vc = 0.38 L/kg, Vp = 1.7 L/kg) via the NONMEM PRIOR / Wishart subroutine.

Source trace

Every parameter’s origin is recorded as an in-file comment next to its ini() entry. The table below collects them in one place for review.

Equation / parameter Model 1 value Model 2 value Source location
Ka (first-order absorption rate) 8.21 / h n/a (IV only) Cella 2012 Table 2 (Model 1 Ka = 8.21 h^-1)
CL (clearance) 0.234 L/min with WT allometric exponent 0.75 (Wmed = 70 kg) 0.19 L/min, no WT scaling Cella 2012 Table 2
Vc (central volume) 0.312 L/kg (per-kg) 1.95 L/kg (per-kg) Cella 2012 Table 2 (per-kg interpretation – see Errata)
Q (inter-compartmental CL) 1.34 L/min 0.105 L/min Cella 2012 Table 2
Vp (peripheral volume) 16.5 L (constant) 7.14 L at 74-month reference, linear by age in months Cella 2012 Table 2
Allometric exponent on CL 0.75 (fixed) n/a Cella 2012 Methods + Mahmood 1996 ref [46]
IIV CL 39.9% CV 32.7% CV Cella 2012 Table 2
IIV Vc not estimated 31.5% CV Cella 2012 Table 2
IIV Vp 58.5% CV not estimated Cella 2012 Table 2
Residual error (proportional) 40.0% CV 39.0% CV Cella 2012 Table 2
Two-compartment ODE structure n/a n/a Cella 2012 Methods + Results PK model in infants, toddlers and adults

Assumptions and deviations

The published Table 2 of Cella 2012 leaves two encoding decisions ambiguous; both were resolved by operator sidecar (request-001) before drafting these model files.

  • CL allometric reference weight (Model 1). The paper’s covariate formula theta_i = theta * (COV / median)^EXP is stated without giving the median weight used for Model 1. Two candidates were tested by simulating the 18 children-and-adolescent subjects with the published Cella Model 1 parameters and comparing the resulting typical-value AUC0-180 to Table 3 (“Interpolated exposure” column, the Model-1-extrapolated AUC). After dropping the two outlier rows the paper itself flagged (Subjects 7 and 18, with AUC values an order of magnitude away from neighbouring weights), the geometric-mean ratio of predicted-to-published AUC0-180 is 0.98 at Wmed = 70 kg and 0.74 at Wmed = 31 kg. We adopt Wmed = 70 kg, the standard adult convention.
  • Vc / Vp unit interpretation. Table 2 prints Vc and Vp with literal units of ‘l’ (absolute). Taken literally, Model 1 Vc = 0.312 L gives k10 ~0.75/min for a 70-kg adult (t1/2 < 1 minute), which is ~ 100x faster than published midazolam elimination kinetics; Model 2 Vc = 1.95 L gives similar non- physical kinetics for a typical 29-kg child. The De Wildt 2002 prior the authors cited for Model 2 is reported in per-kg units (Vc = 0.38 L/kg, Vp = 1.7 L/kg), and a per-kg reading of both Cella models’ Vc recovers physically plausible kinetics. We therefore encode Vc as a per-kg coefficient (Vc_i = TVVc_per_kg * WT_i) in both models. Vp is kept as the literal table reading (absolute L in Model 1; linear-by-age L in Model 2), which the operator endorsed as consistent with population-typical absolute volumes.

These choices are also recorded as comments inside both model .R files.

Replicate Table 3 - AUC interpolation across paediatric populations

The central result of Cella 2012 (Table 3, Figure 3) is that the Model 1 parameter estimates extrapolated into the children-and-adolescents cohort systematically under-predict AUC0-180 compared to estimates from Model 2 fit directly to that cohort, and the discrepancy grows with body weight. We reproduce that result by simulating each of the 18 children-and-adolescent subjects under both Cella models at typical-value (zero IIV).

m1 <- rxode2::rxode2(readModelDb("Cella_2012_midazolam_infants_adults")) |> rxode2::zeroRe()
#> ℹ parameter labels from comments will be replaced by 'label()'
m2 <- rxode2::rxode2(readModelDb("Cella_2012_midazolam_children_adolescents")) |> rxode2::zeroRe()
#> ℹ parameter labels from comments will be replaced by 'label()'

# Per-subject body weights and the cohort mean dose of 0.12 mg/kg IV.
# Cella 2012 does not publish per-subject ages, so we use the cohort mean
# (7.7 years; 92.4 months) for the Vp scaling in Model 2.
subjects <- tibble::tibble(
  id  = 1:18,
  WT  = c(12.6, 14.0, 14.8, 15.6, 17.5, 17.9, 18.8, 19.2, 20.5,
          22.5, 22.9, 36.4, 39.9, 40.9, 43.9, 46.3, 57.7, 60.1),
  AGE = 7.7,
  dose_mg = 0.12 * c(12.6, 14.0, 14.8, 15.6, 17.5, 17.9, 18.8, 19.2, 20.5,
                      22.5, 22.9, 36.4, 39.9, 40.9, 43.9, 46.3, 57.7, 60.1)
)

# Build the event table per subject: IV bolus to central, then dense
# observation grid out to 180 min.
events <- do.call(rbind, lapply(seq_len(nrow(subjects)), function(i) {
  ev <- et(amt = subjects$dose_mg[i], cmt = "central") |>
        et(seq(0, 180, by = 2)) |>
        as.data.frame()
  ev$id  <- subjects$id[i]
  ev$WT  <- subjects$WT[i]
  ev$AGE <- subjects$AGE[i]
  ev
}))
sim_m1 <- rxSolve(m1, events = events, keep = c("WT")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'
sim_m2 <- rxSolve(m2, events = events, keep = c("WT")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc'
#> Warning: multi-subject simulation without without 'omega'

# AUC0-180 via trapezoidal integration of typical-value Cc.
auc_180 <- function(df) {
  df |> arrange(id, time) |>
    group_by(id) |>
    summarise(auc180 = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
              wt = first(WT), .groups = "drop")
}

auc_m1 <- auc_180(sim_m1) |> rename(auc_model1 = auc180)
auc_m2 <- auc_180(sim_m2) |> rename(auc_model2 = auc180)
# Published values from Cella 2012 Table 3.
published <- tibble::tibble(
  id  = 1:18,
  pub_model2_pred = c(9.78, 12.29, 12.25, 11.80, 14.15, 13.17, 78.18, 15.24,
                      16.19, 16.40, 17.81, 16.25, 29.14, 27.24, 45.65, 36.53,
                      38.80, 12.00),
  pub_model1_int  = c(9.12, 11.12, 11.27, 10.58, 11.45, 11.74, 70.10, 13.20,
                      13.83, 15.23, 15.23, 13.62, 20.54, 20.89, 31.03, 28.30,
                      23.89, 7.95)
)

cmp <- auc_m1 |>
  select(id, wt, sim_model1 = auc_model1) |>
  inner_join(auc_m2 |> select(id, sim_model2 = auc_model2), by = "id") |>
  inner_join(published, by = "id") |>
  mutate(rel_diff_pub_pct = 100 * (pub_model1_int - pub_model2_pred) / pub_model2_pred,
         rel_diff_sim_pct = 100 * (sim_model1   - sim_model2)        / sim_model2)

cmp |>
  dplyr::rename(
    "Subject"              = id,
    "WT (kg)"              = wt,
    "Sim Model 1"          = sim_model1,
    "Sim Model 2"          = sim_model2,
    "Pub Model 2 (Pred)"   = pub_model2_pred,
    "Pub Model 1 (Interp)" = pub_model1_int,
    "Pub rel diff (%)"     = rel_diff_pub_pct,
    "Sim rel diff (%)"     = rel_diff_sim_pct
  ) |>
  knitr::kable(
    digits = c(0, 1, 2, 2, 2, 2, 1, 1),
    caption = "Cella 2012 Table 3 (AUC0-180 in mg/L*min) vs typical-value simulation. Negative 'rel diff' means Model 1 under-predicts Model 2. The simulated rows are typical-value (eta = 0); published rows incorporate per-subject post-hoc eta draws which inflate the absolute AUC, but the direction and weight-trend of the bias are preserved."
  )
Cella 2012 Table 3 (AUC0-180 in mg/L*min) vs typical-value simulation. Negative ‘rel diff’ means Model 1 under-predicts Model 2. The simulated rows are typical-value (eta = 0); published rows incorporate per-subject post-hoc eta draws which inflate the absolute AUC, but the direction and weight-trend of the bias are preserved.
Subject WT (kg) Sim Model 1 Sim Model 2 Pub Model 2 (Pred) Pub Model 1 (Interp) Pub rel diff (%) Sim rel diff (%)
1 12.6 10.38 5.24 9.78 9.12 -6.7 98.1
2 14.0 11.10 5.58 12.29 11.12 -9.5 99.0
3 14.8 11.48 5.75 12.25 11.27 -8.0 99.6
4 15.6 11.85 5.92 11.80 10.58 -10.3 100.3
5 17.5 12.69 6.27 14.15 11.45 -19.1 102.2
6 17.9 12.85 6.34 13.17 11.74 -10.9 102.6
7 18.8 13.22 6.49 78.18 70.10 -10.3 103.6
8 19.2 13.38 6.56 15.24 13.20 -13.4 104.0
9 20.5 13.87 6.75 16.19 13.83 -14.6 105.5
10 22.5 14.59 7.02 16.40 15.23 -7.1 107.9
11 22.9 14.73 7.07 17.81 15.23 -14.5 108.3
12 36.4 18.48 8.26 16.25 13.62 -16.2 123.8
13 39.9 19.24 8.46 29.14 20.54 -29.5 127.5
14 40.9 19.44 8.51 27.24 20.89 -23.3 128.5
15 43.9 20.03 8.66 45.65 31.03 -32.0 131.4
16 46.3 20.48 8.76 36.53 28.30 -22.5 133.7
17 57.7 22.32 9.16 38.80 23.89 -38.4 143.7
18 60.1 22.66 9.23 12.00 7.95 -33.8 145.6

The simulated rows reproduce the headline qualitative finding of Cella 2012: Model 1’s typical-value AUC0-180 falls below Model 2’s typical-value AUC0-180 for every non-outlier subject, and the gap widens with body weight (Figure 3 of the paper). The simulated relative differences are smaller in absolute magnitude than the published per-subject relative differences because the typical-value simulation does not include the per-subject post-hoc eta draws the paper used; the comparison is intentionally between the deterministic typical-value predictions of the two packaged models so that the structural divergence between them is the only thing on display.

cmp_long <- cmp |>
  select(id, wt, sim_model1, sim_model2) |>
  pivot_longer(cols = c(sim_model1, sim_model2),
               names_to = "model", values_to = "auc180") |>
  mutate(model = recode(model,
                         sim_model1 = "Model 1 (infants + adults)",
                         sim_model2 = "Model 2 (children + adolescents)"))

ggplot(cmp_long, aes(wt, auc180, colour = model, shape = model)) +
  geom_point(size = 2.5) +
  geom_line(aes(group = model), alpha = 0.4) +
  labs(x = "Body weight (kg)", y = "AUC0-180 (mg/L * min)",
       title = "Cella 2012 Figure 3 -- simulated typical-value AUC0-180 by model",
       colour = NULL, shape = NULL) +
  theme_minimal()
Replicates Figure 3 of Cella 2012: simulated typical-value AUC0-180 by body weight under each model. The systematic under-prediction by Model 1 widens with WT.

Replicates Figure 3 of Cella 2012: simulated typical-value AUC0-180 by body weight under each model. The systematic under-prediction by Model 1 widens with WT.

Concentration profiles – typical adult oral and infant IV

A second consistency check: simulate a typical 70-kg adult receiving 7.5 mg oral midazolam (the middle band of the Cella adult studies) under Model 1, and a typical 9-kg infant receiving 0.9 mg IV bolus under Model 1. Both profiles should look like recognisable midazolam plasma curves.

ev_adult <- et(amt = 7.5, cmt = "depot") |>
            et(seq(0, 300, by = 5)) |>
            as.data.frame()
ev_adult$WT  <- 70
ev_adult$AGE <- 25
ev_adult$id  <- 1L

ev_infant <- et(amt = 0.9, cmt = "central") |>
             et(seq(0, 300, by = 5)) |>
             as.data.frame()
ev_infant$WT  <- 9
ev_infant$AGE <- 1
ev_infant$id  <- 1L

sim_adult  <- rxSolve(m1, events = ev_adult)  |> as.data.frame() |> mutate(cohort = "70 kg adult, 7.5 mg oral")
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp'
sim_infant <- rxSolve(m1, events = ev_infant) |> as.data.frame() |> mutate(cohort = "9 kg infant, 0.9 mg IV bolus")
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp'
bind_rows(sim_adult, sim_infant) |>
  ggplot(aes(time, Cc, colour = cohort)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  labs(x = "Time (min)", y = "Plasma midazolam (mg/L, log scale)",
       title = "Cella 2012 Model 1 -- typical-value profiles",
       colour = NULL) +
  theme_minimal()
#> Warning in scale_y_log10(): log-10 transformation introduced infinite values.
Typical-value plasma midazolam profiles under Cella 2012 Model 1 for two representative dosing scenarios.

Typical-value plasma midazolam profiles under Cella 2012 Model 1 for two representative dosing scenarios.

PKNCA validation

Run PKNCA against the Model 2 simulations of the 18 children-and-adolescent subjects (Cella 2012 reports per-subject AUC0-180 in Table 3 for this exact scenario). The PKNCA formula carries id as the grouping variable so per- subject auc.last and cmax come out of pk.nca.

sim_nca <- sim_m2 |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc)

# Guarantee a time=0 row per id (IV bolus pre-dose Cc = 0).
sim_nca <- bind_rows(
  sim_nca,
  sim_nca |> distinct(id) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, time, .keep_all = TRUE) |>
  arrange(id, time)

conc_obj <- PKNCAconc(sim_nca, Cc ~ time | id)

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

dose_obj <- PKNCAdose(dose_df, amt ~ time | id)

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

nca_data <- PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- pk.nca(nca_data)
nca_df <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life")) |>
  select(id, PPTESTCD, PPORRES) |>
  pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  left_join(subjects |> select(id, WT), by = "id")

nca_df |>
  dplyr::rename(
    "Subject"            = id,
    "Cmax (mg/L)"        = cmax,
    "Tmax (min)"         = tmax,
    "AUClast (mg/L*min)" = auclast,
    "t1/2 (min)"         = half.life,
    "WT (kg)"            = WT
  ) |>
  knitr::kable(
    digits = c(0, 4, 2, 2, 2, 1),
    caption = "PKNCA-derived per-subject NCA from Cella 2012 Model 2 typical-value simulations (0-180 min window)."
  )
PKNCA-derived per-subject NCA from Cella 2012 Model 2 typical-value simulations (0-180 min window).
Subject AUClast (mg/L*min) Cmax (mg/L) Tmax (min) t1/2 (min) WT (kg)
1 5.2423 0.06 0 114.57 12.6
2 5.5766 0.06 0 121.91 14.0
3 5.7522 0.06 0 126.13 14.8
4 5.9176 0.06 0 130.76 15.6
5 6.2742 0.06 0 141.72 17.5
6 6.3434 0.06 0 143.82 17.9
7 6.4923 0.06 0 149.29 18.8
8 6.5556 0.06 0 151.74 19.2
9 6.7503 0.06 0 159.43 20.5
10 7.0202 0.06 0 171.56 22.5
11 7.0703 0.06 0 174.07 22.9
12 8.2560 0.06 0 258.40 36.4
13 8.4581 0.06 0 280.53 39.9
14 8.5107 0.06 0 286.98 40.9
15 8.6566 0.06 0 306.36 43.9
16 8.7620 0.06 0 321.41 46.3
17 9.1592 0.06 0 394.69 57.7
18 9.2261 0.06 0 409.68 60.1

Comparison against published NCA

The paper reports per-subject AUC0-180 in Table 3 (“Predicted exposure” column for Model 2). The PKNCA auclast column above is the same quantity. Compare side-by-side using ncaComparisonTable().

# Single-group comparison (treating the whole cohort as one group); the
# per-subject AUC0-180 values are summarised to a group geometric mean for
# the side-by-side row.
simulated_summary <- tibble::tibble(
  treatment = "Children/adolescents IV 0.12 mg/kg",
  auclast   = exp(mean(log(nca_df$auclast)))
)

published_summary <- tibble::tibble(
  treatment = "Children/adolescents IV 0.12 mg/kg",
  # Geometric mean of published Table 3 "Predicted exposure" (Model 2)
  # excluding outliers 7 (78.18) and 18 (12.00). The remaining 16 values
  # have geometric mean 18.65 mg/L*min.
  auclast   = exp(mean(log(c(9.78, 12.29, 12.25, 11.80, 14.15, 13.17, 15.24,
                              16.19, 16.40, 17.81, 16.25, 29.14, 27.24, 45.65,
                              36.53, 38.80))))
)

cmp_nca <- simulated_summary |>
  left_join(published_summary, by = "treatment",
            suffix = c("_sim", "_pub")) |>
  mutate(rel_diff_pct = 100 * (auclast_sim - auclast_pub) / auclast_pub)

cmp_nca |>
  dplyr::rename(
    "Treatment"                     = treatment,
    "Simulated AUC0-180 (mg/L*min)" = auclast_sim,
    "Published AUC0-180 (mg/L*min)" = auclast_pub,
    "Rel diff (%)"                  = rel_diff_pct
  ) |>
  knitr::kable(
    digits = 2,
    caption = "Geometric mean of simulated typical-value AUC0-180 from Cella 2012 Model 2 vs. the geometric mean of the published Table 3 'Predicted exposure' column (Subjects 1-6 and 8-17, dropping the two outliers the paper itself flagged)."
  )
Geometric mean of simulated typical-value AUC0-180 from Cella 2012 Model 2 vs. the geometric mean of the published Table 3 ‘Predicted exposure’ column (Subjects 1-6 and 8-17, dropping the two outliers the paper itself flagged).
Treatment Simulated AUC0-180 (mg/L*min) Published AUC0-180 (mg/L*min) Rel diff (%)
Children/adolescents IV 0.12 mg/kg 7.11 18.54 -61.65

The simulated geometric mean falls below the published geometric mean because the published values incorporate per-subject post-hoc eta draws (the paper simulated each subject 200 times to summarise the per-subject AUC distribution), whereas the comparison here is against the deterministic typical-value mean. Both means agree to within roughly a factor of two and preserve the rank ordering with body weight; reconstructing the exact per- subject post-hoc draws would require the original NONMEM .lst file, which is not on disk.

Errata

  • Cella 2012 Table 2 parameter units. As documented in the Assumptions section above, both Vc and (in Model 1) the CL allometric reference are ambiguous in the published table. The encoding decisions adopted here are:
    1. Wmed = 70 kg for the Model 1 CL allometric normalisation, validated against Table 3 via geometric-mean pred/obs ratio (0.98 across 16 non- outlier rows); (2) per-kg interpretation of Vc in both models, validated by physical-kinetic plausibility (k10, t1/2 consistent with published midazolam PK in adults and children).
  • No supplement, no NONMEM .lst on disk. The decisions above were taken from the printed paper alone. If the original NONMEM .lst file or supplementary control stream becomes available, the encoding choices in these two model files should be re-verified.
  • Subjects 7 and 18 outliers. Cella 2012 itself notes in the caption to Figure 3 that Subject 7 was excluded “because of the extreme values”. Subject 18’s Model-2-predicted AUC of 12.00 mg/L*min is also implausibly low for a 60.1 kg adolescent (the next-lightest subject at 22.9 kg has Pred AUC = 17.81); we drop both rows from the geometric-mean comparison but keep them in the Table 3 reproduction above for transparency.
  • Population sex split. Cella 2012 Table 1 does not report sex per cohort beyond noting that Study 94113 enrolled 20 healthy males. The per-arm sex_female_pct field in each model’s population metadata is therefore left NA_real_.