---
title: "Tacrolimus (Prytula 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tacrolimus (Prytula 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Prytula_2016_tacrolimus"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Prytula_2016_tacrolimus"))$description`
* Article: <https://doi.org/10.1007/s40262-016-0390-7>

# Population

The model was developed from 120 abbreviated 2-h or 4-h pharmacokinetic
profiles obtained from 54 stable paediatric kidney transplant recipients at
the Erasmus MC-Sophia Children's Hospital in Rotterdam (n = 45) and Emma
Children's Hospital in Amsterdam (n = 9), both in the Netherlands (Prytula
2016 Table 1). Median age was 11.1 years (range 3.8-18.4); median body
weight 38.6 kg (range 15-86); 48% of subjects were female; and 70% were
recorded as Caucasian. All patients were beyond the first year after kidney
transplantation (median 16.2 months posttransplant) on twice-daily oral
tacrolimus (Prograft), with the dose unchanged for at least 14 days at
baseline and titrated by therapeutic drug monitoring to a posttransplant
year >=1 target trough (C0) of 4-8 ug/L. CYP3A5 6986A>G (rs776746) genotype
was determined in 49 of 54 patients (Table 1): *1/*1 n = 1 (2%), *1/*3 n =
12 (24.5%), *3/*3 n = 36 (73.5%); 5 patients had no DNA available and the
source paper estimated a separate CL/F shift for them (h(CYP3A5 missing) =
0.41 with 61% SE) that is not encoded here. Tacrolimus was measured in
whole blood by LC-MS/MS (83% of samples) or MEIA, with the assay-method
covariate tested but not retained.

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

# Source trace

Every parameter in the model file carries an inline source-location comment.
The table below collects the entries in one place.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` (ka at 70 kg) | 0.23 1/h | Table 2, Estimated typical value, ka row |
| `ltlag` (lag time) | 0.47 h | Table 2, Estimated typical value, tlag row |
| `lcl` (CL/F at WT 70 kg, *3/*3, gGT 13, HCT 0.34) | 35 L/h | Table 2, Estimated typical value, CL/F (CYP3A5 *3/*3) row |
| `lvc` (V1/F at 70 kg) | 12 L | Table 2, Estimated typical value, V1/F row |
| `lq` (Q/F at 70 kg) | 68 L/h | Table 2, Estimated typical value, Q/F row |
| `lvp` (V2/F at 70 kg; fixed) | 1090 L | Table 2, Estimated typical value, V2/F row (FIX) |
| `e_cyp3a5_expr_cl` (CYP3A5*1 carrier shift) | 0.45 | Table 2, h(CYP3A5 *1/*1, *1/*3) row |
| `e_ggt_cl` (gGT power exponent) | -0.21 | Table 2, h(gGT) row |
| `e_hct_cl` (HCT power exponent) | -0.59 | Table 2, h(Ht) row |
| `e_wt_cl` (allometric WT on CL/F; fixed) | 0.75 | Section 3.2.1, allometric scaling paragraph |
| `e_wt_q` (allometric WT on Q/F; fixed) | 0.75 | Section 3.2.1, allometric scaling paragraph |
| `e_wt_vc` (allometric WT on V1/F; fixed) | 1 | Section 3.2.1, allometric scaling paragraph |
| `e_wt_vp` (allometric WT on V2/F; fixed) | 1 | Section 3.2.1, allometric scaling paragraph |
| `e_wt_ka` (allometric WT on ka; fixed) | -0.25 | Section 3.2.1, allometric scaling paragraph |
| `q_eta_scale` (IIV-CL-Q scaling) | 2.0 | Table 2, IIV-CL-Q row |
| IIV ka (omega^2 = log(1 + 0.58^2)) | 58% CV | Table 2, Inter-individual variability ka row |
| IIV CL/F (omega^2 = log(1 + 0.45^2)) | 45% CV | Table 2, Inter-individual variability CL/F row |
| IIV V1/F (omega^2 = log(1 + 1.70^2)) | 170% CV | Table 2, Inter-individual variability V1/F row |
| Correlation eta_ka . eta_CL | 0.11 | Table 2 footnote (rho_{ka.CL}) |
| Correlation eta_ka . eta_V1 | 0.19 | Table 2 footnote (rho_{ka.V1}) |
| Correlation eta_CL . eta_V1 | 0.38 | Table 2 footnote (rho_{CL.V1}) |
| Proportional residual error | 21% | Table 2, Residual variability proportional row |
| CL/F covariate equation | -- | Equations on p. 1136 ("CL/F = 35 * (1 + 0.45)^CYP3A5_EXPR * (WT/70)^0.75 * (gGT/13)^-0.21 * (HCT/0.34)^-0.59 * exp(eta_CL + kappa_CL)") |
| Q/F covariate equation | -- | Equations on p. 1136 ("Q/F = 68 * (WT/70)^0.75 * exp(2.0 * eta_CL)") |
| 2-cmt structure, oral first-order absorption + lag | -- | Section 3.2.1, Structural model paragraph |

# Virtual cohort

The published dataset is not openly available, so the virtual cohorts below
mirror the demographics in Prytula 2016 Table 1 and are stratified by the
two pediatric-dosing scenarios the paper highlights in Figure 4 (20 kg and
40 kg children) and by CYP3A5 genotype.

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

n_per_stratum <- 200L

make_cohort <- function(n, wt_target, cyp3a5_expr, label, id_offset = 0L) {
  tibble(
    id          = id_offset + seq_len(n),
    # Body weight: log-normal noise around the target weight; SD chosen to
    # span roughly +/- 25% so each stratum has a narrow weight band around
    # the Figure 4 anchor (20 kg or 40 kg) rather than the full cohort
    # range (15-86 kg).
    WT          = exp(rnorm(n, mean = log(wt_target), sd = 0.12)),
    # gGT: study median 13 U/L (range 4-118); log-normal with moderate
    # spread, truncated to a reasonable clinical band.
    GGT         = pmin(pmax(exp(rnorm(n, mean = log(13), sd = 0.45)), 3), 60),
    # Haematocrit: study median 0.34 (range 0.21-0.44); centred truncated
    # normal.
    HCT         = pmin(pmax(rnorm(n, mean = 0.34, sd = 0.05), 0.20), 0.50),
    CYP3A5_EXPR = cyp3a5_expr,
    cohort      = label
  )
}

# Four CYP3A5 x weight strata for the Figure 4 reproduction. Disjoint IDs.
demo_fig4 <- bind_rows(
  make_cohort(n_per_stratum, wt_target = 20, cyp3a5_expr = 1L,
              label = "20 kg expresser",       id_offset = 0L),
  make_cohort(n_per_stratum, wt_target = 40, cyp3a5_expr = 1L,
              label = "40 kg expresser",       id_offset = 1L * n_per_stratum),
  make_cohort(n_per_stratum, wt_target = 20, cyp3a5_expr = 0L,
              label = "20 kg nonexpresser",    id_offset = 2L * n_per_stratum),
  make_cohort(n_per_stratum, wt_target = 40, cyp3a5_expr = 0L,
              label = "40 kg nonexpresser",    id_offset = 3L * n_per_stratum)
)
stopifnot(!anyDuplicated(demo_fig4$id))
```

# Simulation

Two scenarios are simulated:

1. **Figure 4 reproduction** -- 7 days of 0.1 mg/kg twice daily (0.2 mg/kg/day)
   to reach steady state, with sampling concentrated around the trough of
   the final dosing interval.
2. **Small-child target-trough scenario** -- 7 days of 0.15 mg/kg twice daily
   (0.3 mg/kg/day) for a 15 kg CYP3A5 nonexpresser, to test the paper's
   recommendation that this regimen produces a C0 of 4-8 ug/L in children
   below 20 kg.

```{r build-events}
build_events <- function(demo, mg_per_kg_per_dose, n_days = 7L) {
  # Twice-daily oral dosing -- one dose every 12 h for n_days days.
  doses <- demo |>
    mutate(amt  = round(mg_per_kg_per_dose * WT, 2),
           evid = 1L, cmt = "depot",
           ii   = 12, addl = (2L * n_days) - 1L,
           time = 0) |>
    select(id, time, amt, evid, cmt, ii, addl, cohort,
           WT, GGT, HCT, CYP3A5_EXPR)

  # Observation grid: dense over the final 24 h so we capture the trough
  # and the AUC0-12 interval.
  last_dose_time <- 12 * (2L * n_days - 1L)
  obs_times <- sort(unique(c(seq(0, 12, by = 1),
                              seq(last_dose_time,
                                  last_dose_time + 12, by = 0.25))))
  obs <- demo |>
    select(id, cohort, WT, GGT, HCT, CYP3A5_EXPR) |>
    tidyr::crossing(time = obs_times) |>
    mutate(amt = NA_real_, evid = 0L, cmt = NA_character_,
           ii = NA_real_, addl = NA_integer_)

  bind_rows(doses, obs) |> arrange(id, time, desc(evid))
}

events_fig4 <- build_events(demo_fig4, mg_per_kg_per_dose = 0.1, n_days = 7L)
last_dose_time <- 12 * (2L * 7L - 1L)  # 156 h = day 7 morning dose
```

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

sim_fig4 <- rxode2::rxSolve(
  mod, events = events_fig4,
  keep = c("cohort")
) |> as.data.frame()

mod_typical <- mod |> rxode2::zeroRe()
sim_typ_fig4 <- rxode2::rxSolve(
  mod_typical, events = events_fig4,
  keep = c("cohort")
) |> as.data.frame()
```

# Replicate published figures

## Figure 4 -- trough vs. CYP3A5 genotype and body weight

Prytula 2016 Figure 4 plots simulated tacrolimus trough concentrations
(medians and 25th-75th percentile bands) versus body weight under
twice-daily oral dosing, stratified by CYP3A5 genotype. The text on p. 1138
gives explicit typical-value targets for 0.2 mg/kg/day:

| Stratum | Reported C0 (ug/L) |
|---------|--------------------|
| 20 kg expresser | 5 |
| 40 kg expresser | 6 |
| 20 kg nonexpresser | 8 |
| 40 kg nonexpresser | 10 |

```{r figure-4, fig.cap = "Replicates Figure 4 of Prytula 2016: simulated tacrolimus concentration vs. time after the last day-7 dose, stratified by body weight and CYP3A5 genotype, at 0.1 mg/kg twice daily (0.2 mg/kg/day)."}
fig4_data <- sim_fig4 |>
  filter(time >= last_dose_time, time <= last_dose_time + 12) |>
  mutate(time_after_dose = time - last_dose_time)

fig4 <- fig4_data |>
  group_by(cohort, time_after_dose) |>
  summarise(Q25 = quantile(Cc, 0.25),
            Q50 = quantile(Cc, 0.50),
            Q75 = quantile(Cc, 0.75),
            .groups = "drop")

ggplot(fig4, aes(time_after_dose, Q50)) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.3) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~ cohort) +
  scale_x_continuous(breaks = seq(0, 12, by = 2)) +
  labs(x = "Time after dose (h)",
       y = "Tacrolimus concentration (ng/mL)",
       title = "Day-7 dosing-interval VPC at 0.1 mg/kg twice daily",
       caption = "Replicates Figure 4 of Prytula 2016 (lower-right panel: 0.1 mg/kg/dose).")
```

```{r figure-4-typical-trough}
# Typical-value (zeroRe) C0 at the day-7 trough -- compare to Section 3.2.4
# explicit-target values.
typical_c0 <- sim_typ_fig4 |>
  filter(time == last_dose_time + 12) |>
  group_by(cohort) |>
  summarise(C0_typical = mean(Cc), .groups = "drop")

reported_c0 <- tribble(
  ~cohort,                 ~C0_reported,
  "20 kg expresser",       5,
  "40 kg expresser",       6,
  "20 kg nonexpresser",    8,
  "40 kg nonexpresser",    10
)

typical_c0 |>
  left_join(reported_c0, by = "cohort") |>
  mutate(pct_diff = round(100 * (C0_typical - C0_reported) / C0_reported, 1)) |>
  knitr::kable(
    digits  = 2,
    caption = "Typical-value (zeroRe) trough at day 7 vs. Prytula 2016 Section 3.2.4 reported C0 at 0.2 mg/kg/day."
  )
```

## Small-child scenario: 0.3 mg/kg/day in a <20 kg child

Prytula 2016 Discussion: "In children 1 year or more post-transplantation
whose body weight is lower than 20 kg, the initial dose of 0.3 mg/kg/day
results in a tacrolimus C0 between 4 and 8 ug/L."

```{r small-child}
demo_small <- tibble(
  id          = 1L:200L,
  WT          = exp(rnorm(200L, mean = log(15), sd = 0.10)),
  GGT         = pmin(pmax(exp(rnorm(200L, mean = log(13), sd = 0.45)), 3), 60),
  HCT         = pmin(pmax(rnorm(200L, mean = 0.34, sd = 0.05), 0.20), 0.50),
  CYP3A5_EXPR = 0L,
  cohort      = "15 kg nonexpresser, 0.3 mg/kg/day"
)

events_small <- build_events(demo_small, mg_per_kg_per_dose = 0.15, n_days = 7L)
sim_small <- rxode2::rxSolve(mod, events = events_small,
                             keep = c("cohort")) |> as.data.frame()

trough_small <- sim_small |>
  filter(time == last_dose_time + 12) |>
  summarise(Q10 = quantile(Cc, 0.10),
            Q50 = quantile(Cc, 0.50),
            Q90 = quantile(Cc, 0.90))

knitr::kable(
  trough_small,
  digits  = 2,
  caption = "Simulated day-7 trough for a 15 kg CYP3A5 nonexpresser at 0.15 mg/kg twice daily (0.3 mg/kg/day). Prytula 2016 Discussion reports a 4-8 ug/L target for this scenario."
)
```

# PKNCA validation

NCA is performed over the day-7 dosing interval (treated as steady state).
The treatment grouping variable is the body-weight / CYP3A5 stratum.

```{r pknca}
# Concentration object for PKNCA -- restrict to the last dosing interval and
# re-zero time so PKNCA sees a clean 0-12 h window.
nca_window <- sim_fig4 |>
  filter(time >= last_dose_time, time <= last_dose_time + 12) |>
  mutate(time_after_dose = time - last_dose_time) |>
  select(id, time = time_after_dose, Cc, cohort)

dose_df <- demo_fig4 |>
  mutate(time = 0, amt = round(0.1 * WT, 2)) |>
  select(id, time, amt, cohort)

conc_obj <- PKNCA::PKNCAconc(nca_window, Cc ~ time | cohort + id)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | cohort + id)
intervals <- data.frame(start = 0, end = 12,
                        cmax = TRUE, tmax = TRUE,
                        auclast = TRUE, cmin = TRUE)
nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- suppressMessages(suppressWarnings(PKNCA::pk.nca(nca_data)))
knitr::kable(
  summary(nca_res),
  caption = "Day-7 NCA on the simulated paediatric cohort at 0.1 mg/kg twice daily (steady-state 12-h interval)."
)
```

## Comparison against published AUC12

Prytula 2016 Section 3.3 reports a Bayesian-estimated AUC12 distribution
from 120 observed profiles: median 97 h.ug/L (IQR 80-120, range 39-209)
across the full cohort. The 0.2 mg/kg/day x median-cohort scenario should
land near this central value; cohorts with lighter body weight and
expresser status are expected to fall lower (5 ug/L troughs map to lower
AUC12), and 40 kg nonexpressers higher (10 ug/L troughs map to higher
AUC12).

```{r auc-comparison}
auc_per_cohort <- nca_window |>
  group_by(cohort, id) |>
  summarise(auc12 = PKNCA::pk.calc.auc(conc = Cc, time = time,
                                       interval = c(0, 12),
                                       method = "linear"),
            .groups = "drop") |>
  group_by(cohort) |>
  summarise(median = median(auc12, na.rm = TRUE),
            Q25    = quantile(auc12, 0.25, na.rm = TRUE),
            Q75    = quantile(auc12, 0.75, na.rm = TRUE),
            .groups = "drop")

knitr::kable(
  auc_per_cohort,
  digits = 1,
  caption = "Simulated day-7 AUC0-12 (h.ng/mL) by CYP3A5/body-weight stratum at 0.1 mg/kg twice daily. Prytula 2016 reports a cohort-wide median AUC12 of 97 h.ug/L (IQR 80-120) for trough levels in the 4-8 ug/L target range."
)
```

# Assumptions and deviations

- **CYP3A5 missing-genotype arm omitted.** Prytula 2016 Table 2 carries a
  separate CL/F shift `h(CYP3A5 missing) = 0.41` (SE 61%, bootstrap 95%
  CI -0.08 to 1.52) for the 5 of 54 model-building subjects with no DNA
  available. The estimate is statistically indistinguishable from both the
  carrier effect (0.45) and from zero given the wide confidence interval,
  and is essentially a noise-absorbing parameter for the ungenotyped
  subset; it has no predictive value for downstream simulation when
  CYP3A5 genotype is known. The model file encodes only the
  `CYP3A5_EXPR` binary covariate; users with subjects of unknown
  genotype should impute the most prevalent reference category
  (`CYP3A5_EXPR = 0`) and inspect the impact qualitatively.
- **Intra-individual (inter-occasion) variability omitted.** Prytula 2016
  Table 2 reports an intra-patient variability of 25% CV on CL/F
  representing the within-subject change in CL/F from one year to the
  next. nlmixr2 model files do not have a native IOV construct in
  `ini()` / `model()`, and the source paper does not specify an `OCC`
  indicator pattern. For typical-value simulation the IIV-only form is
  sufficient; for stochastic prediction across visits, users with
  multi-occasion data should implement per-occasion etas in their own
  extension of the model and reduce the IIV variance accordingly.
- **Haematocrit kept as a fraction, not a percent.** The canonical-register
  `HCT` entry uses percent units (e.g. 33%), but Prytula 2016 reports HCT
  as a fraction (e.g. 0.34), and the published power exponent (-0.59
  centred at 0.34) is inseparable from that unit choice. The model
  file's `covariateData[[HCT]]$units` field documents this override; user
  datasets that record HCT in percent must be rescaled
  (`HCT_fraction = HCT_percent / 100`) before passing to this model.
- **Perfect eta_CL / eta_Q correlation encoded via deterministic scaling.**
  Prytula 2016 Section 3.2.1 reports the correlation between
  inter-patient variability on CL/F and Q/F as fixed to 1 with a single
  scaling parameter (IIV-CL-Q = 2.0 in Table 2). The model encodes this
  as `q = exp(lq + q_eta_scale * etalcl) * (WT/70)^0.75` rather than
  carrying a separate `etalq` -- the two formulations are mathematically
  equivalent for a correlation of exactly 1, and the deterministic form
  avoids the singular-covariance numerical issues that a 1-correlation
  block can trigger.
- **ABCB1 polymorphism not encoded.** Prytula 2016 tested ABCB1 3435 C>T
  genotype as a covariate but did not retain it in the final model. The
  cohort's ABCB1 distribution is recorded in `population$abcb1_distribution`
  for completeness but the covariate does not enter `model()`.
- **Race/ethnicity not encoded.** The cohort was 70% Caucasian (Table 1)
  but race was not retained as a covariate; the model has no race effect.
- **Vignette uses 200 subjects per CYP3A5 x weight stratum.** This is
  small enough to render under the 5-minute pkgdown gate but large
  enough to give stable percentile bands for the Figure 4 reproduction.
  The Prytula 2016 simulations used 1000 subjects per scenario.
