Skip to contents

Clegg_2024_nirsevimab

Source: vignettes/articles/Clegg_2024_nirsevimab.Rmd
Clegg_2024_nirsevimab.Rmd
library(nlmixr2lib)
library(PKNCA)
#> 
#> Attaching package: 'PKNCA'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)

Nirsevimab population PK simulations

Replicate Figures 4, 5, and 6 from Clegg et al. (2024):

  • Figure 4: Prediction-corrected VPC for the final model, stratified by trial and RSV season (Phase 2b, MELODY, MEDLEY Season 1, MEDLEY Season 2).
  • Figure 5: Prediction-corrected VPC for the MELODY safety cohort.
  • Figure 6: Predicted exposure (AUC₀₋₃₆₅, C_max, and Day-151 concentration) versus body weight by gestational age (GA) group under weight-banded dosing.

Because the original study data are not publicly available, these figures show simulated prediction intervals from virtual infant populations whose covariate distributions approximate the published trial demographics (Table 1, Clegg 2024). Body weight trajectories are generated using WHO weight-for-age growth standards (combined sex, 0–24 months).

WHO weight-for-age growth curve helper

WHO weight-for-age LMS parameters for combined sex, 0–24 months (WHO Multicentre Growth Reference Study Group, 2006). Formula: weight (kg) = M × (1 + L × S × z)^(1/L)

who_lms <- data.frame(
  age_mo = 0:24,
  L = c(0.3487, 0.2297, 0.1970, 0.1738, 0.1553, 0.1395, 0.1257,
    0.1125, 0.0998, 0.0875, 0.0756, 0.0640, 0.0527, 0.0418,
    0.0313, 0.0211, 0.0113, 0.0018, -0.0073, -0.0161, -0.0245,
    -0.0326, -0.0404, -0.0479, -0.0551),
  M = c(3.3464,  4.4709,  5.5675,  6.3762,  7.0023,  7.5105,  7.9340,
    8.2970,  8.6151,  8.9014,  9.1649,  9.4122,  9.6479,  9.8749,
    10.0953, 10.3108, 10.5228, 10.7319, 10.9385, 11.1430, 11.3462,
    11.5480, 11.7478, 11.9459, 12.1424),
  S = c(0.14602, 0.13395, 0.12385, 0.11876, 0.11535, 0.11254, 0.11056,
    0.10947, 0.10868, 0.10814, 0.10765, 0.10722, 0.10706, 0.10695,
    0.10695, 0.10700, 0.10710, 0.10723, 0.10737, 0.10754, 0.10773,
    0.10793, 0.10814, 0.10835, 0.10858)
)

# Weight (kg) for postnatal age (months) and individual z-score
who_weight <- function(pna_mo, z) {
  pna_mo <- pmax(0, pmin(pna_mo, 24))
  L <- approx(who_lms$age_mo, who_lms$L, xout = pna_mo)$y
  M <- approx(who_lms$age_mo, who_lms$M, xout = pna_mo)$y
  S <- approx(who_lms$age_mo, who_lms$S, xout = pna_mo)$y
  M * (1 + L * S * z)^(1 / L)
}

Helper: build a virtual cohort dataset

Builds dose + observation records for a virtual cohort of n infants. ga_range: gestational age range in weeks. pna0_range: postnatal age at dosing range in months. max_day: follow-up duration in days. season2: TRUE for Season 2 (200 mg flat dose, Season 2 CL effect applied). obs_days: observation time points in days. id_offset: integer added to seq_len(n) so multiple cohorts can be bind_rows()-ed without colliding on ID. rxSolve treats ID as the subject key; duplicate IDs across cohorts silently collapse into a single “Frankenstein” subject that receives the summed dose, so a non-zero offset is mandatory for any multi-cohort simulation (see Figure 4).

make_cohort <- function(n, ga_range, pna0_range, max_day,
                        season2 = FALSE, obs_days = seq(0, max_day, by = 7),
                        id_offset = 0L) {
  GA <- runif(n, ga_range[1], ga_range[2]) # weeks
  PNA_0 <- runif(n, pna0_range[1], pna0_range[2]) # months at dosing
  wt_z <- pmax(-2, pmin(2, rnorm(n, 0, 1)))
  WT_0 <- who_weight(PNA_0, wt_z)
  AMT <- if (season2) rep(200, n) else ifelse(WT_0 < 5, 50, 100)

  pop <- data.frame(
    ID = id_offset + seq_len(n), GA, PNA_0, wt_z, WT_0, AMT,
    SEASON2 = as.integer(season2),
    ADA_POS = 0L,
    # All White/Nat.Haw reference race for simplicity
    RACE_BLACK_OTH = 0L,
    RACE_ASIAN_AMIND_MULTI = 0L
  )

  # Dose records
  d_dose <- pop |>
    mutate(
      TIME = 0, EVID = 1, CMT = "depot", DV = NA,
      PAGE = GA / 4.35 + PNA_0,
      WT = WT_0
    )

  # Observation records with time-varying WT and PAGE
  d_obs <- pop[rep(seq_len(n), each = length(obs_days)), ] |>
    mutate(
      TIME = rep(obs_days, times = n),
      EVID = 0, CMT = "central", DV = NA, AMT = 0,
      PNA = PNA_0 + TIME / 30.44,
      PAGE = GA / 4.35 + PNA,
      WT = who_weight(pmax(0, PNA), wt_z)
    )

  bind_rows(d_dose, d_obs) |>
    arrange(ID, TIME, desc(EVID)) |>
    select(ID, TIME, AMT, EVID, CMT, DV, WT, PAGE,
      RACE_BLACK_OTH, RACE_ASIAN_AMIND_MULTI, SEASON2, ADA_POS, WT_0)
}

Load model 

mod <- readModelDb("Clegg_2024_nirsevimab")
conc_unit <- rxode2::rxode(mod)$units[["concentration"]]
#> ℹ parameter labels from comments will be replaced by 'label()'

Figure 4 — Prediction-corrected VPC by trial and season

Four panels matching the published trials:

Panel Trial GA range Season Dose
Phase 2b Healthy preterm infants 29–35 wk 1 50/100 mg
MELODY Healthy term/late-preterm ≥35 wk 1 50/100 mg
MEDLEY S1 Preterm + CHD/CLD 24–35 wk 1 50/100 mg
MEDLEY S2 CHD/CLD children ≥35 wk 2 200 mg

The vignette shows simulated prediction intervals only (no observed data; original data are not publicly available).

set.seed(8897) # MEDI8897 = nirsevimab development code

# `id_offset` per cohort is mandatory: rxSolve uses ID as the subject key,
# so duplicate IDs across cohorts silently merge into a single subject that
# receives the summed dose. Without offsets, Figure 4's predictions came out
# ~3-fold too high because IDs 1..100 appeared in all four cohorts and were
# collapsed into Frankenstein subjects with ~400 mg in depot at TIME = 0.
d_p2b <- make_cohort(200, c(29, 35), c(0, 3), 400, id_offset =   0L) |>
  mutate(trial = "Phase 2b\n(healthy infants, 29 to <35 wGA)")
d_melody <- make_cohort(200, c(35, 42), c(0, 3), 500, id_offset = 200L) |>
  mutate(trial = "MELODY\n(healthy infants, \u226535 wGA)")
d_med_s1 <- make_cohort(200, c(24, 35), c(0, 3), 400, id_offset = 400L) |>
  mutate(trial = "MEDLEY Season 1\n(infants \u226435 wGA, CHD/CLD)")
d_med_s2 <- make_cohort(100, c(35, 42), c(12, 24), 400,
  season2 = TRUE, id_offset = 600L) |>
  mutate(trial = "MEDLEY Season 2\n(children with CHD/CLD)")

sim_f4 <- bind_rows(d_p2b, d_melody, d_med_s1, d_med_s2)
# Cheap regression guard: every (ID, TIME, EVID) triple must be unique.
stopifnot(!anyDuplicated(unique(sim_f4[, c("ID", "TIME", "EVID")])))

# Carry `trial` through rxSolve via `keep =` rather than a post-hoc left_join.
# A `left_join(distinct(ID, trial))` fans out every row when an ID happens to
# repeat across trials, silently re-pooling cohorts; `keep = "trial"` aligns
# the trial label per output row.
out_f4 <- rxode2::rxSolve(mod, events = sim_f4, keep = "trial") |>
  as.data.frame()
d_f4 <- out_f4 |>
  group_by(trial, time) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(trial = factor(trial, levels = c(
    "Phase 2b\n(healthy infants, 29 to <35 wGA)",
    "MELODY\n(healthy infants, \u226535 wGA)",
    "MEDLEY Season 1\n(infants \u226435 wGA, CHD/CLD)",
    "MEDLEY Season 2\n(children with CHD/CLD)"
  )))
ggplot(d_f4, aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), fill = "#4682b4", alpha = 0.25) +
  geom_line(colour = "#4682b4", linewidth = 0.8) +
  facet_wrap(~trial, scales = "free_x", nrow = 2) +
  scale_y_log10(
    limits = c(0.1, NA),
    labels = scales::label_number(drop0trailing = TRUE)
  ) +
  scale_x_continuous(breaks = seq(0, 500, by = 100)) +
  labs(
    x = "Time after dose (days)",
    y = paste0("Nirsevimab serum concentration (", conc_unit, ")\nPrediction corrected"),
    title = "Figure 4 \u2014 Prediction-corrected VPC by trial and season",
    caption = paste0(
      "Simulated median (line) and 90% prediction interval (shaded).\n",
      "Original observed data not available; reference race = White, ADA negative."
    )
  ) +
  theme_bw() +
  theme(strip.text = element_text(size = 9))
#> Warning in scale_y_log10(limits = c(0.1, NA), labels = scales::label_number(drop0trailing = TRUE)): log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.

Figure 5 — Prediction-corrected VPC for the MELODY safety cohort

Healthy term and late-preterm infants (≥35 weeks GA), Season 1, followed for 360 days. The paper uses 10th/50th/90th percentiles (rather than the 5th/95th used in Figure 4).

set.seed(3979) # NCT03979313 = MELODY trial
d_f5 <- make_cohort(500, c(35, 42), c(0, 3), 360,
  obs_days = seq(0, 360, by = 7))

out_f5 <- rxode2::rxSolve(mod, events = d_f5) |> as.data.frame()
d_f5_plot <- out_f5 |>
  group_by(time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  )
ggplot(d_f5_plot, aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), fill = "#4682b4", alpha = 0.25) +
  geom_line(colour = "#4682b4", linewidth = 0.8) +
  geom_line(aes(y = Q10), linetype = "dashed", linewidth = 0.5) +
  geom_line(aes(y = Q90), linetype = "dashed", linewidth = 0.5) +
  scale_y_log10(
    limits = c(0.5, NA),
    labels = scales::label_number(drop0trailing = TRUE)
  ) +
  scale_x_continuous(breaks = seq(0, 360, by = 30)) +
  labs(
    x = "Time after dose (days)",
    y = paste0("Nirsevimab serum concentration (", conc_unit, ")\nPrediction corrected"),
    title = "Figure 5 \u2014 Prediction-corrected VPC: MELODY safety cohort",
    caption = paste0(
      "Simulated 10th/50th/90th percentiles (N = 500 virtual infants, \u226535 wGA).\n",
      "Original observed data not available."
    )
  ) +
  theme_bw()
#> Warning in scale_y_log10(limits = c(0.5, NA), labels = scales::label_number(drop0trailing = TRUE)): log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.

Figure 6 — Predicted exposure versus body weight by GA group

Simulate 3 GA groups (GA <29 weeks, ≥29 to <35 weeks, ≥35 weeks), dosing at birth, with weight-banded dosing (50 mg ≤5 kg, 100 mg >5 kg). Compute AUC₀₋₃₆₅ (trapezoidal), C_max, and Day-151 concentration. Use daily time steps for accurate NCA.

Adult reference lines from the paper (3000 mg IV dose): AUC₀₋₃₆₅ = 57.0 day·mg/mL; C_max = 1140 μg/mL. EC₉₀ (preclinical) = 6.8 μg/mL (for Day-151 concentration panel).

set.seed(2024)
n_per_ga <- 200
obs_daily <- seq(0, 365, by = 1)

ga_groups <- list(
  list(label = "GA < 29 weeks", ga = c(24, 29)),
  list(label = "GA \u226529 to <35 weeks", ga = c(29, 35)),
  list(label = "GA \u226535 weeks", ga = c(35, 42))
)

sim_f6_list <- lapply(ga_groups, function(g) {
  GA <- runif(n_per_ga, g$ga[1], g$ga[2])
  PNA_0 <- rep(0, n_per_ga) # dose at birth
  wt_z <- pmax(-2, pmin(2, rnorm(n_per_ga, 0, 1)))
  WT_0 <- who_weight(PNA_0, wt_z)
  AMT <- ifelse(WT_0 <= 5, 50, 100)

  pop <- data.frame(
    ID = seq_len(n_per_ga), GA, PNA_0, wt_z, WT_0, AMT,
    SEASON2 = 0L, ADA_POS = 0L,
    RACE_BLACK_OTH = 0L, RACE_ASIAN_AMIND_MULTI = 0L
  )

  d_dose <- pop |>
    mutate(TIME = 0, EVID = 1, CMT = "depot", DV = NA,
      PAGE = GA / 4.35, WT = WT_0)

  d_obs <- pop[rep(seq_len(n_per_ga), each = length(obs_daily)), ] |>
    mutate(
      TIME = rep(obs_daily, times = n_per_ga),
      EVID = 0, CMT = "central", DV = NA, AMT = 0,
      PNA = PNA_0 + TIME / 30.44,
      PAGE = GA / 4.35 + PNA,
      WT = who_weight(pmax(0, PNA), wt_z)
    )

  d_events <- bind_rows(d_dose, d_obs) |>
    arrange(ID, TIME, desc(EVID)) |>
    select(ID, TIME, AMT, EVID, CMT, DV, WT, PAGE,
      RACE_BLACK_OTH, RACE_ASIAN_AMIND_MULTI, SEASON2, ADA_POS, WT_0)

  out <- rxode2::rxSolve(mod, events = d_events) |>
    as.data.frame()

  # NCA by individual: AUC0_365 and Cmax via PKNCA; C_day151 as direct lookup
  out_nca <- out |> filter(time >= 0, !is.na(Cc)) |> mutate(ga_group = g$label)

  # Build PKNCA dose data from the per-subject dose amounts
  dose_df <- pop |>
    mutate(time = 0, ga_group = g$label) |>
    select(id = ID, time, amt = AMT, ga_group)

  conc_obj <- PKNCAconc(out_nca, Cc ~ time | ga_group + id)
  dose_obj <- PKNCAdose(dose_df, amt ~ time | ga_group + id)
  nca_res  <- pk.nca(PKNCAdata(conc_obj, dose_obj,
    intervals = data.frame(start = 0, end = 365,
                           auclast = TRUE, cmax = TRUE)))

  nca_pknca <- as.data.frame(nca_res$result) |>
    filter(PPTESTCD %in% c("auclast", "cmax")) |>
    select(id, PPTESTCD, PPORRES) |>
    tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
    rename(AUC0_365 = auclast, Cmax = cmax)

  # Day-151 concentration: direct lookup (not a standard NCA parameter)
  c_day151 <- out |>
    group_by(id) |>
    summarise(C_day151 = Cc[which.min(abs(time - 151))],
              WT_base  = WT_0[1], .groups = "drop")

  nca <- left_join(nca_pknca, c_day151, by = "id") |>
    mutate(ga_group = g$label)
})

d_f6_nca <- bind_rows(sim_f6_list) |>
  mutate(ga_group = factor(ga_group, levels = sapply(ga_groups, `[[`, "label")))
# 90% prediction intervals (5th-95th) by baseline weight bin
wt_bins <- seq(1, 11, by = 0.5)

d_f6_long <- d_f6_nca |>
  pivot_longer(c(AUC0_365, Cmax, C_day151), names_to = "metric", values_to = "value") |>
  mutate(wt_bin = cut(WT_base, breaks = wt_bins, right = FALSE))

d_f6_pi <- d_f6_long |>
  group_by(ga_group, metric, wt_bin) |>
  summarise(
    wt_mid = mean(WT_base),
    Q05 = quantile(value, 0.05, na.rm = TRUE),
    Q50 = quantile(value, 0.50, na.rm = TRUE),
    Q95 = quantile(value, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  filter(!is.na(wt_bin))

metric_labels <- c(
  AUC0_365 = "AUC[0-365]~(day%.%mg/mL)",
  Cmax     = "C[max]~(mu*g/mL)",
  C_day151 = "C[day~151]~(mu*g/mL)"
)
# Adult reference values (3000 mg IV)
adult_refs <- data.frame(
  metric = c("AUC0_365", "Cmax"),
  yref   = c(57.0, 1140)
)
ec90_val <- 6.8 # preclinical EC90, μg/mL

ggplot(d_f6_pi, aes(x = wt_mid)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), fill = "#87ceeb", alpha = 0.6) +
  geom_line(aes(y = Q50), linewidth = 0.8) +
  # Individual posthoc predictions
  geom_point(
    data = d_f6_nca |>
      pivot_longer(c(AUC0_365, Cmax, C_day151), names_to = "metric", values_to = "value"),
    aes(x = WT_base, y = value),
    colour = "grey50", size = 0.5, alpha = 0.4
  ) +
  # Adult reference lines (AUC and Cmax panels only)
  geom_hline(
    data = adult_refs,
    aes(yintercept = yref),
    linetype = "dotted", linewidth = 0.7
  ) +
  # EC90 line (Day-151 panel only)
  geom_hline(
    data = data.frame(metric = "C_day151", yref = ec90_val),
    aes(yintercept = yref),
    linetype = "dashed", linewidth = 0.7
  ) +
  facet_grid(
    metric ~ ga_group,
    scales = "free_y",
    labeller = labeller(metric = as_labeller(metric_labels, label_parsed))
  ) +
  scale_x_continuous(breaks = seq(1, 11, by = 2)) +
  labs(
    x = "Baseline weight (kg)",
    y = NULL,
    title = "Figure 6 \u2014 Predicted exposure vs body weight by GA group",
    caption = paste0(
      "Blue band: 90% prediction interval (5th\u201395th percentile). ",
      "Black line: median. Grey dots: individual predictions.\n",
      "Dotted line: adult reference (3000 mg IV). ",
      "Dashed line: preclinical EC\u2090 = 6.8 \u03bcg/mL (Day-151 panel).\n",
      "White race, ADA negative, Season 1 assumed."
    )
  ) +
  theme_bw() +
  theme(
    strip.text.x = element_text(size = 9),
    strip.text.y = element_text(size = 8)
  )
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?

Notes on the simulation

  • Model: 2-compartment with first-order IM absorption; published parameters are for a 70 kg adult reference, scaled allometrically to infant weights.
  • Virtual populations: Covariate distributions match published trial demographics (Table 1). Race is set to the White/Native Hawaiian reference and ADA to negative for all simulated subjects to isolate weight and maturation effects.
  • Weight trajectories: WHO weight-for-age LMS growth standards (combined sex, 0–24 months), with each infant’s growth percentile (z-score) held constant.
  • IIV: Simulated using published omega² values: CL (26% CV), V2 (43% CV), Ka (44% CV; high shrinkage 83%), with correlation r = 0.785 between CL and V2. Q and V3 have no IIV in the final model.
  • Residual error: Proportional 21%.
  • Figure 6 NCA: AUC₀₋₃₆₅ and C_max computed via PKNCA on daily simulated concentrations; Day-151 concentration derived by direct lookup from the same densely sampled time course.

Reference

  • Clegg L, Freshwater E, Leach A, Villafana T, Wahlby Hamren U. Population Pharmacokinetics of Nirsevimab in Preterm and Term Infants. J Clin Pharmacol. 2024;64(5):555-567. doi:10.1002/jcph.2401