---
title: "Nirsevimab (Clegg 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Nirsevimab (Clegg 2024)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>"
)
```

```{r setup}
library(nlmixr2lib)
library(PKNCA)
library(dplyr)
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)

```{r who-lms}
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).

```{r make-cohort}
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

```{r load-model}
mod <- readModelDb("Clegg_2024_nirsevimab")
conc_unit <- rxode2::rxode(mod)$units[["concentration"]]
```

---

## 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).

```{r fig4-sim, message=FALSE}
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.
# Cohort sizes downsampled from 200/200/200/100 for vignette build budget;
# VPC ribbons remain smooth at n=30-40 per panel.
d_p2b <- make_cohort(40, 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(40, c(35, 42), c(0, 3), 500, id_offset = 200L) |>
  mutate(trial = "MELODY\n(healthy infants, \u226535 wGA)")
d_med_s1 <- make_cohort(40, 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(30, 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()
```

```{r fig4-summarise}
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)"
  )))
```

```{r fig4-plot, fig.width=8, fig.height=7}
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))
```

---

## 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).

```{r fig5-sim, message=FALSE}
set.seed(3979) # NCT03979313 = MELODY trial
# Downsampled from 500 for vignette build budget; 10-50-90 VPC band shape preserved.
d_f5 <- make_cohort(80, 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()
```

```{r fig5-summarise}
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"
  )
```

```{r fig5-plot, fig.width=6.5, fig.height=4}
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 = 80 virtual infants, \u226535 wGA).\n",
      "Original observed data not available."
    )
  ) +
  theme_bw()
```

---

## 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).

```{r fig6-sim, message=FALSE}
set.seed(2024)
# Downsampled from 200 per GA group for vignette build budget; daily obs grid
# preserved (PKNCA AUC0-365 needs dense sampling). Weight-binned PI bands remain
# the right shape at n=25/group.
n_per_ga <- 25
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(!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")))
```

```{r fig6-summarise}
# 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)"
)
```

```{r fig6-plot, fig.width=9, fig.height=8}
# 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)
  )
```

## 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

* `r rxode2::rxode(readModelDb("Clegg_2024_nirsevimab"))$reference`
