---
title: "Propranolol (Del Frari 2018)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Propranolol (Del Frari 2018)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("DelFrari_2018_propranolol"))$reference`
* Description: `r rxode2::rxode2(readModelDb("DelFrari_2018_propranolol"))$description`
* Article: <https://doi.org/10.1002/prp2.399>

# Population

Del Frari 2018 developed a one-compartment population PK model for oral
propranolol in 22 infants with proliferative Infantile Hemangiomas
recruited from four French hospitals (Bordeaux, Lyon, Nantes, Nice;
EudraCT 2009-018102-22). Infants were stratified into two groups by age
at inclusion: Group 1 (n = 10) was 35-90 days old at inclusion with
the steady-state PK assessment performed after 4 weeks of treatment
(D28); Group 2 (n = 12) was 91-150 days old at inclusion with the
steady-state PK assessment performed after 12 weeks of treatment (D84).
167 plasma concentrations contributed to the final analysis: 21
predose troughs at D7 and 18 at D14 during titration, then 60
intensive samples in Group 1 at D28 and 68 in Group 2 at D84 (predose,
1, 2, 4, 6, and 9 h post-dose).

The pooled cohort had median age 104 days (range 50-151) and median
weight 6.3 kg across all visits D1-D84. Mean weight at the PK
assessment was 5.62 kg in Group 1 (D28) and 7.78 kg in Group 2 (D84).
Treatment was a pediatric oral solution of propranolol hydrochloride
titrated over 2 weeks (1 mg/kg/day, then 2 mg/kg/day) to the target
3 mg/kg/day BID dose for weeks 3-12. Sample baseline demographics are
in Del Frari 2018 Table 2.

The same information is available programmatically:
`readModelDb("DelFrari_2018_propranolol")$population`.

# Source trace

Per-parameter origin (also recorded as in-file comments next to each
`ini()` entry of
`inst/modeldb/specificDrugs/DelFrari_2018_propranolol.R`):

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl`                | `log(19.3)` | Del Frari 2018 Table 5: CL/F = 19.3 L/h at WT = 6.3 kg (95% CI 15.4-23.2, RSE 10.3%) |
| `lvc`                | `log(122)`  | Del Frari 2018 Table 5: V/F = 122 L (95% CI 104-140, RSE 7.32%); no covariate retained |
| `lka`                | `log(0.993)`| Del Frari 2018 Table 5: Ka = 0.993 1/h (95% CI 0.558-1.43, RSE 22.4%) |
| `e_wt_cl`            | `fixed(0.75)` | Del Frari 2018 Results 3.4: allometric exponent on CL/F fixed to 0.75 on AIC grounds; reference WT = 6.3 kg per Equation 5 |
| `etalcl`             | `0.195` | Del Frari 2018 Table 5: omega^2 on CL = 0.195 (44.2% CV; 95% CI 0.0492-0.341, RSE 38.2%) |
| `etalka`             | `1.75`  | Del Frari 2018 Table 5: omega^2 on Ka = 1.75 (132% CV; 95% CI 0.370-3.13, RSE 40.2%) |
| `propSd`             | `sqrt(0.0953) = 0.309` | Del Frari 2018 Table 5: proportional sigma^2 = 0.0953 (30.9% CV; 95% CI 0.0675-0.123, RSE 14.9%) |
| `(WT / 6.3)^0.75` on CL | n/a | Del Frari 2018 Equation 5; reference weight 6.3 kg = pooled median across D1-D84 (paragraph below Eq 5 in Section 3.4) |
| `d/dt(depot)`, `d/dt(central)` | n/a | Del Frari 2018 Results 3.3: best structural model is a one-compartment disposition with first-order absorption and first-order elimination |
| `Cc <- central / vc * 1000` | n/a | Dose mg / volume L -> mg/L; ng/mL = mg/L * 1000 to match the paper's reporting units |
| `Cc ~ prop(propSd)` | n/a | Del Frari 2018 Results 3.3: residual variability described by a proportional error model |

# Virtual cohort

Original observed concentrations are not publicly available. The
simulation below replicates the Del Frari 2018 study cohort: 10 Group 1
infants at the D28 median weight (5.6 kg) and 12 Group 2 infants at the
D84 median weight (7.65 kg). Each infant is replicated 50 times by
Monte Carlo to give 1100 simulated trajectories per dosing regimen
(the source paper replicated the n = 22 cohort 1000 times to derive
Table 6 percentiles; 50 replicates keeps the vignette inside the
pkgdown 5-minute per-vignette wall-time budget while preserving the
qualitative Cmin / Cmax features).

```{r cohort}
set.seed(20260529L)

n_replicates <- 50L
tau          <- 12      # h; primary BID dosing interval
n_days_ss    <- 7L      # 7 days to reach steady state (~12 elimination half-lives)
total_dur    <- n_days_ss * 24

# Per-subject weights for the published cohort. Group 1 median weight at
# D28 (5.6 kg) and Group 2 median weight at D84 (7.65 kg) per Del Frari
# 2018 Table 2 D28 and D84 rows.
subject_specs <- dplyr::bind_rows(
  tibble::tibble(group = "Group 1", WT = 5.6,  n = 10L),
  tibble::tibble(group = "Group 2", WT = 7.65, n = 12L)
)

# Build the per-replicate subject set with disjoint IDs.
build_subjects <- function() {
  ids <- integer(0)
  wt  <- numeric(0)
  grp <- character(0)
  off <- 0L
  for (i in seq_len(nrow(subject_specs))) {
    n_i <- subject_specs$n[i]
    ids <- c(ids, off + seq_len(n_i))
    wt  <- c(wt,  rep(subject_specs$WT[i], n_i))
    grp <- c(grp, rep(subject_specs$group[i], n_i))
    off <- off + n_i
  }
  tibble::tibble(id = ids, WT = wt, group = grp)
}

# Build dosing + observation events for a single regimen.
# `dose_times_h` is the vector of clock-times-since-midnight at which doses
# are taken on a single day; this is repeated for `n_days_ss` days, with
# dose amount = total_daily_mgkg / length(dose_times_h) * WT per dose.
make_cohort <- function(regimen_label, dose_times_h,
                        total_daily_mgkg = 3, id_offset = 0L) {
  subj <- build_subjects() |>
    dplyr::mutate(rep_id = id_offset + id) |>
    dplyr::select(id = rep_id, WT, group)

  per_dose_mgkg <- total_daily_mgkg / length(dose_times_h)
  dose_offsets  <- tidyr::expand_grid(
    day  = 0:(n_days_ss - 1L),
    hour = dose_times_h
  ) |>
    dplyr::mutate(time = day * 24 + hour) |>
    dplyr::pull(time)

  dose_rows <- tidyr::expand_grid(
    id   = subj$id,
    time = dose_offsets
  ) |>
    dplyr::left_join(subj, by = "id") |>
    dplyr::mutate(
      amt  = per_dose_mgkg * WT,
      evid = 1L,
      cmt  = 1L
    )

  obs_grid <- seq(0, total_dur, by = 0.5)
  obs_rows <- tidyr::expand_grid(
    id   = subj$id,
    time = obs_grid
  ) |>
    dplyr::left_join(subj, by = "id") |>
    dplyr::mutate(
      amt  = 0,
      evid = 0L,
      cmt  = NA_integer_
    )

  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::mutate(regimen = regimen_label) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

regimens <- list(
  list(label = "Regular BID (8h, 20h)",   doses = c(8, 20)),
  list(label = "Irregular BID (8h, 17h)", doses = c(8, 17)),
  list(label = "TID (8h, 14h, 20h)",      doses = c(8, 14, 20)),
  list(label = "TID (8h, 12h, 20h)",      doses = c(8, 12, 20))
)

# Build all 50 replicates of each regimen with disjoint IDs across the
# whole simulation.
n_subj_per_rep <- sum(subject_specs$n)
events <- do.call(
  dplyr::bind_rows,
  lapply(seq_along(regimens), function(r_idx) {
    reg <- regimens[[r_idx]]
    do.call(
      dplyr::bind_rows,
      lapply(seq_len(n_replicates), function(rep_idx) {
        off <- ((r_idx - 1L) * n_replicates + (rep_idx - 1L)) * n_subj_per_rep
        make_cohort(
          regimen_label    = reg$label,
          dose_times_h     = reg$doses,
          total_daily_mgkg = 3,
          id_offset        = off
        )
      })
    )
  })
)

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode2(readModelDb("DelFrari_2018_propranolol"))

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("WT", "group", "regimen")
) |>
  as.data.frame()
```

# Typical-clearance check (Section 3.4 anchor values)

Del Frari 2018 Section 3.4 reports the typical CL/F for the median
infant in each group: 17.2 L/h at 5.6 kg (Group 1, D28) and 22.3 L/h
at 7.65 kg (Group 2, D84). The packaged model uses the published
CL/F = 19.3 L/h at the pooled reference weight 6.3 kg with a fixed
allometric exponent of 0.75 on body weight. The table below reproduces
the per-group typical clearance directly from the model's `lcl` and
`e_wt_cl` parameters.

```{r typical-cl}
ini_vals <- mod$ini

cl_typical <- function(wt, lcl_val, e_wt_val, ref_wt = 6.3) {
  exp(lcl_val) * (wt / ref_wt) ^ e_wt_val
}

lcl_value  <- ini_vals$est[ini_vals$name == "lcl"]
ewt_value  <- ini_vals$est[ini_vals$name == "e_wt_cl"]

typical_cl_tbl <- tibble::tribble(
  ~group,     ~wt_kg, ~published_cl_Lh,
  "Group 1",     5.6,             17.2,
  "Group 2",    7.65,             22.3
) |>
  dplyr::mutate(
    sim_cl_Lh = round(cl_typical(wt_kg, lcl_value, ewt_value), 2),
    pct_diff  = round(100 * (sim_cl_Lh - published_cl_Lh) / published_cl_Lh, 1)
  )

knitr::kable(
  typical_cl_tbl,
  caption = paste(
    "Typical apparent oral clearance CL/F (L/h) at the per-group median",
    "weight, compared to the values reported in Del Frari 2018 Section 3.4."
  )
)
```

# Replicate Figure 7: simulated steady-state profiles under 4 dosing regimens

Del Frari 2018 Figure 7 plots the 5th, 50th, and 95th percentile
envelopes of simulated propranolol concentration over a 24-hour
steady-state window for four dosing regimens at a total daily dose of
3 mg/kg/day: (A) regular BID 8:00 + 20:00, (B) irregular BID 8:00 +
17:00, (C) TID 8:00 + 14:00 + 20:00, and (D) TID 8:00 + 12:00 + 20:00.
The figure below reproduces those envelopes from the packaged model.

```{r figure-7, fig.width = 9, fig.height = 6}
last_day_start <- (n_days_ss - 1L) * 24
sim_ss <- sim |>
  dplyr::filter(time >= last_day_start, !is.na(Cc)) |>
  dplyr::mutate(time_of_day = time - last_day_start)

ss_pi <- sim_ss |>
  dplyr::group_by(regimen, time_of_day) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(ss_pi, aes(time_of_day, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(linewidth = 0.6, colour = "steelblue4") +
  facet_wrap(~ regimen) +
  labs(
    x = "Time of day (h since 00:00 on day 7)",
    y = "Propranolol plasma conc. (ng/mL)",
    title = "Replicates Del Frari 2018 Figure 7",
    caption = paste(
      "Median + 5-95% PI at steady state on day 7,",
      "1100 simulated trajectories per regimen (n = 22 cohort x 50 replicates)."
    )
  )
```

## Comparison against Table 6 simulated Cmin / Cmax percentiles

Del Frari 2018 Table 6 reports the 5th / 50th / 95th percentiles of
simulated Cmin and Cmax over a 24-hour steady-state window for the
four dosing regimens above. The table below reproduces those
percentiles from the packaged model and renders them alongside the
published anchors. Per Del Frari 2018 Section 3.6, simulated median
Cmax differed by less than 10% across regimens and median Cmin by
less than 20%; the packaged-model simulations should reproduce that
qualitative finding.

```{r table-6-compare}
cmin_cmax_sim <- sim_ss |>
  dplyr::group_by(regimen, id) |>
  dplyr::summarise(
    Cmin = min(Cc, na.rm = TRUE),
    Cmax = max(Cc, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::group_by(regimen) |>
  dplyr::summarise(
    sim_cmin_p05 = round(stats::quantile(Cmin, 0.05, na.rm = TRUE), 2),
    sim_cmin_p50 = round(stats::quantile(Cmin, 0.50, na.rm = TRUE), 2),
    sim_cmin_p95 = round(stats::quantile(Cmin, 0.95, na.rm = TRUE), 2),
    sim_cmax_p05 = round(stats::quantile(Cmax, 0.05, na.rm = TRUE), 2),
    sim_cmax_p50 = round(stats::quantile(Cmax, 0.50, na.rm = TRUE), 2),
    sim_cmax_p95 = round(stats::quantile(Cmax, 0.95, na.rm = TRUE), 2),
    .groups = "drop"
  )

published_t6 <- tibble::tribble(
  ~regimen,                       ~pub_cmin_p05, ~pub_cmin_p50, ~pub_cmin_p95, ~pub_cmax_p05, ~pub_cmax_p50, ~pub_cmax_p95,
  "Regular BID (8h, 20h)",                 1.99,          14.0,          39.4,          24.0,          61.8,           131,
  "Irregular BID (8h, 17h)",              0.803,          11.6,          47.0,          26.4,          67.1,           141,
  "TID (8h, 14h, 20h)",                    1.43,          14.2,          50.2,          24.7,          62.2,           130,
  "TID (8h, 12h, 20h)",                    1.53,          15.2,          43.4,          23.6,          60.1,           128
)

t6_compare <- dplyr::left_join(cmin_cmax_sim, published_t6, by = "regimen")

knitr::kable(
  t6_compare,
  caption = paste(
    "Simulated steady-state Cmin and Cmax percentiles (5 / 50 / 95)",
    "by dosing regimen, alongside Del Frari 2018 Table 6 published anchors.",
    "All concentrations in ng/mL."
  )
)
```

# PKNCA validation

Steady-state NCA on the regular-BID regimen, computing Cmax, Tmax,
Cmin (trough), and AUC over the 12-hour dosing interval at the end of
day 7.

```{r pknca}
ss_bid <- sim |>
  dplyr::filter(regimen == "Regular BID (8h, 20h)", time >= (n_days_ss - 1L) * 24 + 8, time <= (n_days_ss - 1L) * 24 + 20, !is.na(Cc)) |>
  dplyr::mutate(
    tad       = time - ((n_days_ss - 1L) * 24 + 8),
    treatment = group
  )

ss_dose <- ss_bid |>
  dplyr::group_by(id, treatment, WT) |>
  dplyr::summarise(.groups = "drop") |>
  dplyr::mutate(
    time = 0,
    amt  = 1.5 * WT     # 1.5 mg/kg per BID dose
  )

conc_obj <- PKNCA::PKNCAconc(
  ss_bid |> dplyr::transmute(id, time = tad, Cc, treatment),
  Cc ~ time | treatment + id
)

dose_obj <- PKNCA::PKNCAdose(
  ss_dose |> dplyr::transmute(id, time, amt, treatment),
  amt ~ time | treatment + id,
  route = "extravascular"
)

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

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

nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "cmin", "tmax", "auclast", "half.life")) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(
    median = round(stats::median(PPORRES, na.rm = TRUE), 2),
    p05    = round(stats::quantile(PPORRES, 0.05, na.rm = TRUE), 2),
    p95    = round(stats::quantile(PPORRES, 0.95, na.rm = TRUE), 2),
    .groups = "drop"
  ) |>
  tidyr::pivot_wider(
    names_from  = PPTESTCD,
    values_from = c(median, p05, p95),
    names_glue  = "{PPTESTCD}_{.value}"
  ) |>
  dplyr::arrange(treatment)

knitr::kable(
  nca_summary,
  caption = paste(
    "Steady-state NCA over the 12-hour BID dosing interval on day 7 by",
    "study group (median; 5-95% PI). Cmax / Cmin in ng/mL,",
    "Tmax / half.life in h, auclast in ng/mL*h."
  )
)
```

## Comparison against Table 3 observed steady-state concentrations

Del Frari 2018 Table 3 reports the observed median (range) propranolol
plasma concentration at predose and 1, 2, 4, 6, and 9 h after the
morning dose at D28 (Group 1) and D84 (Group 2). The table below
extracts the corresponding times from the simulated steady-state Day 7
profile and compares to the published observed-data medians.

```{r table-3-compare}
key_times <- c(0, 1, 2, 4, 6, 9)
sim_t3 <- sim_ss |>
  dplyr::filter(
    regimen == "Regular BID (8h, 20h)",
    time_of_day >= 8,
    (time_of_day - 8) %in% key_times
  ) |>
  dplyr::mutate(time_after_dose = time_of_day - 8) |>
  dplyr::group_by(group, time_after_dose) |>
  dplyr::summarise(
    sim_median = round(stats::median(Cc, na.rm = TRUE), 2),
    .groups = "drop"
  )

published_t3 <- tibble::tribble(
  ~group,    ~time_after_dose, ~pub_median,
  "Group 1",                0,        17.5,
  "Group 1",                1,        62.2,
  "Group 1",                2,        61.8,
  "Group 1",                4,        47.0,
  "Group 1",                6,        37.2,
  "Group 1",                9,        33.8,
  "Group 2",                0,        9.74,
  "Group 2",                1,        41.9,
  "Group 2",                2,        74.3,
  "Group 2",                4,        44.2,
  "Group 2",                6,        31.4,
  "Group 2",                9,        21.7
)

t3_compare <- dplyr::left_join(sim_t3, published_t3,
                               by = c("group", "time_after_dose")) |>
  dplyr::mutate(
    pct_diff = round(100 * (sim_median - pub_median) / pub_median, 1)
  )

knitr::kable(
  t3_compare,
  caption = paste(
    "Simulated steady-state median propranolol concentration at the Del Frari",
    "2018 Table 3 sampling times, compared against the published observed",
    "medians from Group 1 (D28) and Group 2 (D84). All values in ng/mL."
  )
)
```

# Assumptions and deviations

- **Cohort weight stratification at the per-group D28 / D84 median.**
  The packaged validation vignette simulates 10 Group 1 infants at the
  D28 median weight (5.6 kg) and 12 Group 2 infants at the D84 median
  weight (7.65 kg), matching the sample sizes and per-group steady-
  state assessment timing of the source clinical study (Del Frari 2018
  Table 2 D28 and D84 rows). Within-group weight variability is held
  at the median rather than sampled from the observed range to keep
  the residual variability dominated by the packaged-model
  between-subject IIV on CL/F and Ka, which is what the Figure 7 and
  Table 6 percentile envelopes reflect.
- **Allometric reference weight 6.3 kg.** Del Frari 2018 Equation 5
  and the paragraph immediately below it define the allometric
  reference as 6.3 kg (the median weight in the pooled population
  from D1 to D84). The same value is hard-coded in the model file's
  `(WT / 6.3)^e_wt_cl` expression.
- **Allometric exponent fixed at 0.75.** Del Frari 2018 Results
  Section 3.4 paragraph 2 reports that the fixed-0.75 exponent yielded
  a lower AIC than the estimated-exponent variant; the packaged model
  encodes the exponent as `fixed(0.75)` to preserve this status. A
  downstream user re-fitting the model can release this constraint
  but the packaged simulation should not.
- **No BSV on V/F.** Del Frari 2018 Results Section 3.4 paragraph 4
  removed BSV on V/F from the final model because the 95% CI of
  omega^2 included 0 and the eta-shrinkage was 62.3%. The packaged
  model has no `etalvc` term, consistent with the final Table 5
  parameterisation; V/F is a fixed-effect-only parameter.
- **No covariate effect on V/F.** Del Frari 2018 Section 3.4
  paragraph 3 tested but did not retain a weight effect on V/F. The
  packaged model has no `(WT / 6.3)^e_wt_vc` scaling on the volume of
  distribution.
- **Cohort N for Figure 7 reproduction.** Del Frari 2018 used 1000
  replicates of the n = 22 cohort (~22,000 trajectories per regimen);
  the packaged vignette uses 50 replicates (1100 trajectories per
  regimen) to fit inside the pkgdown 5-minute per-vignette wall-time
  budget. The smaller N adds some Monte Carlo noise to the percentile
  bands but preserves the qualitative findings (Cmin / Cmax differ by
  less than 20% / 10% across regimens; the BID-regular and TID-regular
  envelopes overlap substantially).
- **Day 7 used as the steady-state day.** The source paper simulated
  one week of repeated dosing at the target 3 mg/kg/day dose to reach
  steady state before recording Cmin / Cmax; with a terminal-phase
  half-life of ~5-6 h (V/F = 122 L, CL/F = 19.3 L/h at 6.3 kg yields
  t1/2 = ln(2) * V / CL = 4.4 h; longer with allometric scaling at
  lower body weights), 7 days = 14 doses is well past 5 half-lives
  for any plausible WT.
- **Residual error reported on the variance scale.** Del Frari 2018
  Table 5 reports the proportional residual variance sigma^2 = 0.0953
  per the NONMEM $SIGMA convention; the packaged model uses
  `propSd = sqrt(0.0953) = 0.309` to obtain the equivalent
  linear-scale proportional SD (~30.9% CV).
- **Bioavailability not identifiable.** Del Frari 2018 Methods
  Section 2.3.1 paragraph 1 notes that absolute bioavailability could
  not be determined from oral-only data; the structural parameters
  reported in Table 5 are therefore apparent (CL/F, V/F) and the
  packaged model has no explicit F1 parameter. Downstream users
  performing dose-to-exposure simulations should treat the dose
  variable as an oral mg amount; the apparent-volume scaling absorbs
  the missing F term.
