---
title: "Sapropterin (Feillet 2008)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sapropterin (Feillet 2008)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

```{r model-header, echo = FALSE, results = "asis"}
mod <- rxode2::rxode(readModelDb("Feillet_2008_sapropterin"))
cat("* **Citation:** ", mod$reference, "\n", sep = "")
cat("* **Description:** ", mod$description, "\n", sep = "")
```

* **Article:** <https://doi.org/10.2165/0003088-200847120-00006>

# Population

Feillet 2008 characterised the population pharmacokinetics of oral
sapropterin dihydrochloride (Kuvan) in 78 adolescent and adult patients
with BH4-responsive phenylketonuria (PKU) enrolled in the 12-week
fixed-dose phase of the open-label PKU-004 extension study. The trial
was conducted at 26 centres across North America (Canada, USA) and
Europe (France, Germany, Ireland, Italy, Poland, UK). Patients were
randomised in the parent PKU-003 study to receive daily oral doses of
5 (n = 6), 10 (n = 37), or 20 mg/kg (n = 34); the dose was not recorded
for one patient (excluded from analysis). Baseline demographics per
Feillet 2008 Table II: age 9-50 years (mean 21.1, SD 9.64), body weight
28.2-144 kg (mean 67.2, SD 21.8), 45 (58%) males, 76 (97%) White, mean
baseline plasma phenylalanine 811 umol/L (SD 393, range 53-2190).
Median (range) creatinine clearance was 114 (48-231) mL/min; only
patients with normal or mildly impaired renal function were enrolled.

A D-optimal sparse sampling design was used, with four samples per
patient during weeks 6, 10, and 12. The final analysis dataset
contained 265 plasma BH4 observations from 76 patients (315
observations minus 38 below the limit of quantification, 4 unreported,
and 8 excluded as unreliable). BH4 concentrations were measured
indirectly via oxidation of the analyte to L-biopterin followed by
reversed-phase HPLC with tandem mass spectrometry; the nominal BH4 to
L-biopterin conversion ratio was 47.3%. Bodyweight was the only
covariate retained at p<0.001 out of a screen of sex, race, age,
height, body surface area, serum creatinine, albumin, ALT, AST, total
bilirubin, and phenylalanine.

The same information is available programmatically via
`readModelDb("Feillet_2008_sapropterin")()$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to
each `ini()` entry in
`inst/modeldb/specificDrugs/Feillet_2008_sapropterin.R`. The table
below collects them in one place for review.

| Equation / parameter | Value              | Source location (Feillet 2008)                                     |
|----------------------|--------------------|--------------------------------------------------------------------|
| `ltlag`              | `log(0.275)`       | Table III: tlag = 0.275 h                                          |
| `lka`                | `log(0.518)`       | Table III: ka = 0.518 h^-1                                         |
| `lcl`                | `log(2100)`        | Table III: CL/F = 2100 L/h/70 kg                                   |
| `lvc`                | `log(8350)`        | Table III: V1/F = 8350 L/70 kg                                     |
| `lvp`                | `log(4240)`        | Table III: V2/F = 4240 L (not weight-scaled)                       |
| `lq`                 | `log(862)`         | Table III: Q/F = 862 L/h (not weight-scaled)                       |
| `lrbase`             | `log(13.5)`        | Table III: BASE = 13.5 ng/mL                                       |
| `e_wt_cl`            | `0.586`            | Table III: power function on CL/F = 0.586                          |
| `e_wt_vc`            | `1.13`             | Table III: power function on V1/F = 1.13                           |
| IIV CL/F (omega^2)   | `0.539^2 = 0.290521` | Table III: SD of IIV for CL/F = 0.539                            |
| IIV V1/F (omega^2)   | `0.557^2 = 0.310249` | Table III: SD of IIV for V1/F = 0.557                            |
| cov(CL,V1) on log    | `0.336 * 0.539 * 0.557 = 0.100869` | Table III: correlation R = 0.336 between CL/F and V1/F |
| `propSd`             | `0.217`            | Table III: constant CV residual error = 21.7%                      |
| Equation: CL/F       | `theta * (WT/70)^0.586 * exp(eta)` | Equation 3 (power-form covariate model, reference 70 kg)   |
| Equation: V1/F       | `theta * (WT/70)^1.13  * exp(eta)` | Equation 3 (power-form covariate model, reference 70 kg)   |
| Equation: BASE       | `theta` (no IIV, no covariate)     | Section 2.4 and Table III (BASE has no IIV)                |
| ODE: 2-compartment   | first-order absorption + linear disposition | Fig. 1 (schematic) and Results (Section 3.1)      |
| Observation          | `Cc = 1000 * central / vc + rbase` | Fig. 1 + Table III BASE additive to plasma BH4             |
| Residual             | proportional (constant CV)         | Section 2.4 and Table III                                  |

The paper also reports mean (SD) individual half-lives from the final
model that are used below as external validation targets (Section
3.2): initial (distribution) half-life 1.45 (0.47) h and terminal
(elimination) half-life 6.69 (2.29) h.

# Virtual cohort

Individual observed data are not public. The cohort below approximates
the PKU-004 fixed-dose phase: adolescent and adult patients whose body
weight distribution matches Feillet 2008 Table II (67.2 +/- 21.8 kg,
range 28.2-144), evenly distributed across the three dose groups. Per
the vignette-cohort cap, 200 subjects are simulated per dose arm.

```{r cohort}
set.seed(20082008)

make_cohort <- function(n,
                        dose_mg_per_kg,
                        n_doses = 5L,
                        dose_interval_hr = 24,
                        obs_times_post_dose_hr = c(0, 0.25, 0.5, 1, 1.5, 2,
                                                   3, 4, 6, 8, 12, 16, 20, 24),
                        weight_mean = 67.2,
                        weight_sd   = 21.8,
                        weight_min  = 28.2,
                        weight_max  = 144,
                        id_offset   = 0L) {

  WT <- pmax(weight_min,
             pmin(weight_max, rnorm(n, weight_mean, weight_sd)))

  pop <- data.frame(
    id        = id_offset + seq_len(n),
    WT        = WT,
    treatment = paste0(dose_mg_per_kg, " mg/kg")
  )

  dose_times <- seq(0, (n_doses - 1) * dose_interval_hr,
                    by = dose_interval_hr)

  d_dose <- pop[rep(seq_len(n), each = length(dose_times)), ] |>
    dplyr::mutate(
      time = rep(dose_times, times = n),
      amt  = dose_mg_per_kg * WT,
      evid = 1L,
      cmt  = "depot",
      dv   = NA_real_
    )

  last_dose_time <- (n_doses - 1) * dose_interval_hr
  obs_grid <- sort(unique(c(
    dose_times[1] + obs_times_post_dose_hr,
    last_dose_time + obs_times_post_dose_hr
  )))

  d_obs <- pop[rep(seq_len(n), each = length(obs_grid)), ] |>
    dplyr::mutate(
      time = rep(obs_grid, times = n),
      amt  = 0,
      evid = 0L,
      cmt  = "central",
      dv   = NA_real_
    )

  dplyr::bind_rows(d_dose, d_obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid)) |>
    dplyr::select(id, time, amt, evid, cmt, dv, WT, treatment)
}

events <- dplyr::bind_rows(
  make_cohort(n = 200L, dose_mg_per_kg =  5, id_offset =   0L),
  make_cohort(n = 200L, dose_mg_per_kg = 10, id_offset = 200L),
  make_cohort(n = 200L, dose_mg_per_kg = 20, id_offset = 400L)
)

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

# Simulation

```{r simulate}
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("treatment", "WT")) |>
  as.data.frame() |>
  dplyr::mutate(
    treatment      = factor(treatment,
                            levels = c("5 mg/kg", "10 mg/kg", "20 mg/kg")),
    time_post_last = time - 96   # last (5th) dose at t = 96 h
  )
```

# Replicate published figures

## Figure 3 analogue: visual predictive check of the final covariate model

Feillet 2008 Figure 3 shows a visual predictive check of the final
covariate model with 5000 simulated patients: the observed BH4
concentrations plotted against time after dose, overlaid with the
10th, 50th, and 90th simulated percentiles. Panels a-c below reproduce
the equivalent VPC bands stratified by dose group across the first
dosing interval.

```{r fig3-vpc, fig.width = 8, fig.height = 4}
sim |>
  dplyr::filter(time >= 0, time <= 24) |>
  dplyr::group_by(treatment, time) |>
  dplyr::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(aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90),
              fill = "#4682b4", alpha = 0.25) +
  geom_line(colour = "#4682b4", linewidth = 0.7) +
  facet_wrap(~ treatment) +
  labs(
    x       = "Time after dose (h)",
    y       = "BH4 concentration (ng/mL)",
    title   = "Figure 3 analogue: VPC of plasma BH4 by dose group",
    caption = "Simulated 10/50/90 percentiles from 200 subjects per arm.\nReplicates Figure 3 of Feillet 2008."
  ) +
  theme_bw()
```

## Figure 4 analogue: no accumulation of BH4 across five daily doses

Feillet 2008 Figure 4 shows stochastic simulations of the BH4
concentration-time profile across five daily doses at 5, 10, and
20 mg/kg to illustrate the absence of accumulation at any dose. We
reproduce the same window (0 to 120 h) using the packaged model.

```{r fig4-noaccum, fig.width = 8, fig.height = 4}
sim |>
  dplyr::group_by(treatment, time) |>
  dplyr::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(aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90),
              fill = "#b00000", alpha = 0.25) +
  geom_line(colour = "#b00000", linewidth = 0.7) +
  facet_wrap(~ treatment) +
  labs(
    x       = "Time (h)",
    y       = "BH4 concentration (ng/mL)",
    title   = "Figure 4 analogue: five daily doses of oral sapropterin",
    caption = "Simulated 10/50/90 percentiles across 5 x 24 h daily dosing.\nTrough at t=120 h should be close to the pre-dose baseline.\nReplicates Figure 4 of Feillet 2008."
  ) +
  theme_bw()
```

Feillet 2008 reports "little evidence of accumulation, even at the
highest dose". The median trough at the end of the fifth dosing
interval should therefore return close to the endogenous baseline.

```{r fig4-troughs}
troughs <- sim |>
  dplyr::group_by(treatment, id) |>
  dplyr::summarise(
    trough_dose1 = Cc[which.min(abs(time -  24))],
    trough_dose5 = Cc[which.min(abs(time - 120))],
    .groups = "drop"
  ) |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    median_trough_dose1 = median(trough_dose1),
    median_trough_dose5 = median(trough_dose5),
    ratio               = median_trough_dose5 / median_trough_dose1,
    .groups = "drop"
  )

troughs |>
  dplyr::rename(
    "Dose group"                = treatment,
    "Median trough after dose 1 (ng/mL)" = median_trough_dose1,
    "Median trough after dose 5 (ng/mL)" = median_trough_dose5,
    "Ratio (dose5 / dose1)"     = ratio
  ) |>
  knitr::kable(
    digits  = 2,
    caption = paste0(
      "Median trough plasma BH4 after the first vs. fifth daily dose. ",
      "A ratio close to 1 confirms no meaningful accumulation."
    )
  )
```

## Body-weight relationship on CL/F and V1/F

Feillet 2008 uses a power-form covariate model (Equation 3, reference
70 kg) with exponents 0.586 on CL/F and 1.13 on V1/F. The analytical
relationships are:

```{r wt-relationships, fig.width = 7, fig.height = 4}
wt_grid  <- seq(28, 145, length.out = 201)
cl_grid  <- 2100 * (wt_grid / 70)^0.586
vc_grid  <- 8350 * (wt_grid / 70)^1.13

p_cl <- ggplot(data.frame(wt = wt_grid, cl = cl_grid),
               aes(x = wt, y = cl)) +
  geom_line(colour = "#4682b4", linewidth = 0.8) +
  geom_vline(xintercept = 70, linetype = "dotted") +
  labs(x = "Body weight (kg)", y = "CL/F (L/h)",
       title = "CL/F vs body weight (Feillet 2008 Eq. 3)") +
  theme_bw()

p_vc <- ggplot(data.frame(wt = wt_grid, vc = vc_grid),
               aes(x = wt, y = vc)) +
  geom_line(colour = "#4682b4", linewidth = 0.8) +
  geom_vline(xintercept = 70, linetype = "dotted") +
  labs(x = "Body weight (kg)", y = "V1/F (L)",
       title = "V1/F vs body weight (Feillet 2008 Eq. 3)") +
  theme_bw()

if (requireNamespace("patchwork", quietly = TRUE)) {
  patchwork::wrap_plots(p_cl, p_vc, ncol = 2)
} else {
  print(p_cl)
  print(p_vc)
}
```

# PKNCA validation

Feillet 2008 does not report classical NCA parameters (Cmax, AUC), but
does report mean (SD) individual half-lives: initial (distribution)
1.45 (0.47) h and terminal (elimination) 6.69 (2.29) h (Section 3.2).
We simulate a typical-value single-dose profile (20 mg/kg into a 70 kg
adult) and compute the terminal half-life via PKNCA over a 48-hour
window, subtracting the endogenous BH4 baseline first so the NCA
reflects drug-derived exposure.

```{r pknca-single}
mod_typical <- rxode2::zeroRe(mod)

obs_times_single <- c(0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 5, 6,
                      8, 10, 12, 16, 20, 24, 30, 36, 48)

ev_single <- data.frame(
  id   = 1L,
  time = c(0, obs_times_single),
  amt  = c(20 * 70, rep(0, length(obs_times_single))),
  evid = c(1L,      rep(0L, length(obs_times_single))),
  cmt  = c("depot", rep("central", length(obs_times_single))),
  dv   = NA_real_,
  WT   = 70,
  treatment = "single_20mgkg_70kg"
)

sim_single <- rxode2::rxSolve(mod_typical, events = ev_single,
                              keep = c("treatment", "WT")) |>
  as.data.frame() |>
  dplyr::mutate(id = 1L)  # rxSolve drops id for single-subject sims

# Feillet 2008 reports drug-only half-lives; subtract the endogenous
# BH4 baseline (BASE = 13.5 ng/mL) before running PKNCA.
c0_typical <- 13.5
sim_single <- sim_single |>
  dplyr::mutate(Cc_drug = pmax(Cc - c0_typical, 0))

sim_nca <- sim_single |>
  dplyr::filter(!is.na(Cc_drug)) |>
  dplyr::select(id, time, Cc_drug, treatment) |>
  dplyr::rename(Cc = Cc_drug)

# Time-zero guarantee.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

dose_df <- ev_single |>
  dplyr::filter(evid == 1) |>
  dplyr::transmute(id = id, time = time, amt = amt,
                   treatment = "single_20mgkg_70kg")

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

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

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

## Comparison against published half-lives

```{r nca-comparison}
nca_tbl <- as.data.frame(nca_res$result)

get_param <- function(res_df, ppname) {
  val <- res_df$PPORRES[res_df$PPTESTCD == ppname]
  if (length(val) == 0) return(NA_real_)
  val[1]
}

hl_obs_sim <- get_param(nca_tbl, "half.life")
cmax_sim   <- get_param(nca_tbl, "cmax")
tmax_sim   <- get_param(nca_tbl, "tmax")
auc_sim    <- get_param(nca_tbl, "aucinf.obs")

# Analytical half-lives at WT = 70 kg from the packaged thetas.
ka_t   <- 0.518
cl_70  <- 2100
vc_70  <- 8350
vp_70  <- 4240
q_70   <- 862
kel_70 <- cl_70 / vc_70
k12_70 <- q_70  / vc_70
k21_70 <- q_70  / vp_70
sum_r  <- kel_70 + k12_70 + k21_70
prod_r <- kel_70 * k21_70
disc   <- sqrt(pmax(sum_r^2 - 4 * prod_r, 0))
alpha  <- (sum_r + disc) / 2   # fast (distribution)
beta   <- (sum_r - disc) / 2   # slow (elimination)
hl_alpha <- log(2) / alpha
hl_beta  <- log(2) / beta

comparison <- data.frame(
  Quantity = c(
    "Initial (alpha) half-life at WT = 70 kg (h)",
    "Terminal (beta) half-life at WT = 70 kg (h)",
    "PKNCA terminal half-life from simulated profile (h)",
    "PKNCA Cmax (ng/mL, drug-only)",
    "PKNCA Tmax (h)",
    "PKNCA AUCinf (ng*h/mL, drug-only)"
  ),
  Published = c(
    "1.45 (SD 0.47)",
    "6.69 (SD 2.29)",
    "-",
    "-",
    "~2 (peak occurs approximately 2 h post-dose per Discussion)",
    "-"
  ),
  Analytical = c(
    sprintf("%.2f", hl_alpha),
    sprintf("%.2f", hl_beta),
    "-",
    "-",
    "-",
    "-"
  ),
  Simulated_typical = c(
    "-",
    "-",
    sprintf("%.2f", hl_obs_sim),
    sprintf("%.1f", cmax_sim),
    sprintf("%.2f", tmax_sim),
    sprintf("%.0f", auc_sim)
  ),
  check.names = FALSE
)

comparison |>
  knitr::kable(
    caption = paste0(
      "Half-lives, Cmax, Tmax, and AUCinf from the packaged model at a ",
      "typical 70 kg adult given 20 mg/kg oral sapropterin, vs the mean ",
      "individual values in Feillet 2008 Section 3.2. Cmax / AUCinf are ",
      "drug-only (endogenous BASE = 13.5 ng/mL subtracted before NCA)."
    )
  )
```

The analytical terminal half-life reproduces Feillet 2008's mean
individual terminal half-life of 6.69 h to within ~10%; the initial
half-life sits close to the reported 1.45 h; and the simulated tmax of
~2 h matches the paper's description of "peak concentrations occurring
approximately 2 hours after dosing" (Discussion). Small residual
differences arise because the paper's reported half-lives are means
across individual weights (28.2-144 kg) while the analytical values
here are evaluated at the 70 kg reference used in the population
covariate model.

# Assumptions and deviations

- **Endogenous BH4 baseline added at the observation step.** Feillet
  2008 models endogenous BH4 as a constant additive offset (BASE)
  on the predicted plasma concentration with no IIV and no
  covariate structure (Table III lists BASE as a fixed-effect only
  parameter). The packaged model encodes this as
  `Cc = (central / vc) * 1000 + rbase` where `rbase = exp(lrbase)`
  and `lrbase = log(13.5)`. For the PKNCA half-life comparison in
  this vignette, BASE is subtracted from the simulated Cc so the
  NCA reflects drug-derived exposure only (the paper's reported
  half-lives are drug-derived).
- **No IIV on tlag, ka, V2/F, or Q/F.** Feillet 2008 explicitly
  states in the Discussion that "the model did not include any
  variability in ka or lag time (tlag) which might explain the
  slight underprediction at approximately 2.5 hours post-dose". The
  packaged encoding follows the paper: only CL/F and V1/F carry
  ETAs, with a correlation R = 0.336 between them (Table III).
- **Sex, race, and other covariates screened but not retained.**
  Feillet 2008 evaluated sex, race, age, height, body surface area,
  serum creatinine, albumin, ALT, AST, total bilirubin, and
  phenylalanine as covariates but found only bodyweight
  statistically significant at p<0.001. The packaged model therefore
  exposes only `WT` as a structural covariate.
- **Reference weight = 70 kg.** The paper reports CL/F and V1/F with
  the units "L/h/70 kg" and "L/70 kg" in Table III, which we
  interpret as a normalisation to a 70 kg adult reference in the
  power covariate model of Equation 3.
- **V2/F and Q/F are not weight-scaled.** Feillet 2008 does not
  attach the "/70 kg" normalisation to V2/F or Q/F in Table III and
  states in the Results that bodyweight was retained on CL/F and
  V1/F only. The packaged model consequently treats V2/F and Q/F as
  population-wide constants that do not vary with weight.
- **Virtual cohort weight distribution.** Individual patient weights
  are not published. The vignette draws a truncated normal weight
  distribution whose mean (67.2 kg) and range (28.2-144 kg) match
  Feillet 2008 Table II, evenly distributing 200 subjects across
  each of the three dose groups.
- **Apparent parameterisation (no separate bioavailability).**
  Sapropterin is administered orally as sapropterin dihydrochloride
  but the assay measures the pharmacologically active BH4 in
  plasma; the paper reports apparent parameters (CL/F, V1/F, V2/F,
  Q/F) that absorb both bioavailability and the nominal 47.3% drug
  to L-biopterin conversion factor of the assay. The packaged
  model reproduces the same apparent parameterisation; no
  `lfdepot` bioavailability parameter is exposed.

# Reference

```{r ref-footer, echo = FALSE, results = "asis"}
cat("* ", mod$reference, "\n", sep = "")
```
