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

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Andrews_2017_tacrolimus"))$meta$reference`
* Description: `r rxode2::rxode(readModelDb("Andrews_2017_tacrolimus"))$meta$description`
* Article: <https://doi.org/10.1007/s40262-017-0567-8>

Andrews et al. (2017) developed a two-compartment population PK model for
twice-daily oral immediate-release tacrolimus (Prograft / Modigraf) in 46
paediatric kidney transplant recipients during the first 6 weeks
post-transplantation. The model used allometric scaling at theory-based fixed
exponents and identified four significant covariates on apparent oral
clearance: CYP3A5 expresser status, donor source (deceased vs living), eGFR
(adapted-Schwartz, BSA-normalised), and a piecewise hematocrit effect that
engages only when HCT < 0.3 L/L. The published dosing guideline (Table 4)
recommends starting doses ranging from 0.27 to 1.33 mg/kg/day depending on
weight, CYP3A5 genotype, and donor source. This vignette reproduces the
typical-value structural model, simulates the cohort-level dose-effect
profiles shown in Figure 3, and validates NCA-derived AUC0-12h against the
target-trough mapping in Section 3.5.

# Population

The model-building cohort (Andrews 2017 Table 1) was n = 46 paediatric
kidney transplant recipients (median age 9.1 years, range 2.4-17.9 years;
median weight 28.4 kg, range 11.6-83.7 kg; 43.5% female; 73.9% Caucasian,
13.0% Black, 4.3% Asian, 8.7% Other). All children were treated with the
TWIST immunosuppressive protocol (basiliximab + tacrolimus + mycophenolic
acid + a 5-day course of glucocorticoids). The initial tacrolimus dose was
0.3 mg/kg/day divided into two doses every 12 h, subsequently titrated by
therapeutic drug monitoring to a target pre-dose concentration of 10-15
ng/mL during the first 3 weeks and 7-12 ng/mL thereafter. Tacrolimus
concentrations were measured by validated LC-MS/MS (91% of samples; LLOQ
1.0 ng/mL) or by immunoassay (9%; pre-LC-MS/MS-introduction at the
laboratory; LLOQ 1.5 ng/mL). External validation was performed on an
independent cohort of n = 23 children from the Radboud University Medical
Center, Nijmegen.

The same information is available programmatically via the model's
`population` metadata
(`rxode2::rxode(readModelDb("Andrews_2017_tacrolimus"))$meta$population`).

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` (ka)           | log(0.56) (1/h) | Andrews 2017 Table 2 final model |
| `ltlag` (tlag)       | log(0.37) (h)   | Andrews 2017 Table 2 final model |
| `lcl` (CL/F)         | log(50.5) (L/h) | Andrews 2017 Table 2 final model |
| `lvc` (V1/F)         | log(206) (L)    | Andrews 2017 Table 2 final model |
| `lq` (Q/F)           | log(114) (L/h)  | Andrews 2017 Table 2 final model |
| `lvp` (V2/F)         | log(1520) (L)   | Andrews 2017 Table 2 final model |
| `e_wt_cl` (allometric CL) | fixed(0.75) | Andrews 2017 Section 3.1 |
| `e_wt_q` (allometric Q)   | fixed(0.75) | Andrews 2017 Section 3.1 |
| `e_wt_vc` (allometric V1) | fixed(1)    | Andrews 2017 Section 3.1 |
| `e_wt_vp` (allometric V2) | fixed(1)    | Andrews 2017 Section 3.1 |
| `e_cyp3a5_nonexpr_cl` | 1.04 | Andrews 2017 Table 2 final model |
| `e_cyp3a5_expr_cl`    | 1.98 | Andrews 2017 Table 2 final model |
| `e_donor_living_cl`   | 0.74 | Andrews 2017 Table 2 final model |
| `e_egfr_cl` (eGFR exponent) | 0.19 | Andrews 2017 Table 2 final model |
| `e_hct_cl` (HCT exponent, piecewise) | -0.44 | Andrews 2017 Table 2 final model |
| IIV ka (188% CV)      | omega^2 = log(1 + 1.88^2) = 1.51178 | Andrews 2017 Table 2 final model |
| IIV CL/F (25% CV)     | omega^2 = log(1 + 0.25^2) = 0.06062 | Andrews 2017 Table 2 final model |
| IIV V1/F (69% CV)     | omega^2 = log(1 + 0.69^2) = 0.38944 | Andrews 2017 Table 2 final model |
| IIV V2/F (62% CV)     | omega^2 = log(1 + 0.62^2) = 0.32509 | Andrews 2017 Table 2 final model |
| `addSd_immuno` / `propSd_immuno` | 1.01 / 0.13 | Andrews 2017 Table 2 final model |
| `addSd_lcms` / `propSd_lcms`     | 0.28 / 0.21 | Andrews 2017 Table 2 final model |
| Reference HCT cutoff  | 0.3 L/L | Andrews 2017 Section 3.2 equation |
| Reference eGFR        | 69 mL/min/1.73 m^2 | Andrews 2017 Table 1 cohort median |
| Reference body weight | 70 kg | Andrews 2017 Section 3.1 (theory-based) |
| 2-compartment ODE structure | d/dt(depot), d/dt(central), d/dt(peripheral1) | Andrews 2017 Section 3.1 |

# Virtual cohort

The original observed data are not publicly available. The figures below use
virtual populations whose covariate distributions approximate the published
trial demographics (Table 1).

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

# Helper: build one cohort as a self-contained event table for a given body
# weight, CYP3A5 expresser status, donor source, eGFR, hematocrit, and assay.
# Dosing is 0.3 mg/kg/day split q12h for 14 days, matching the Andrews 2017
# starting-dose simulation in Figure 3 / Table 4.
make_cohort <- function(n, wt, cyp3a5_expr, donor_deceased, egfr, hct, immuno,
                        dose_mg_per_kg = 0.3, n_days = 14,
                        regimen = "0.3 mg/kg/day q12h",
                        id_offset = 0L) {
  per_dose <- dose_mg_per_kg * wt / 2   # mg, twice daily
  dose_times <- seq(from = 0, by = 12, length.out = n_days * 2)
  obs_times  <- sort(unique(c(seq(0, n_days * 24, by = 1),
                              dose_times + c(0.5, 1, 2, 4, 6, 8))))

  per_subject <- function(i) {
    id_i <- id_offset + i
    doses <- tibble::tibble(
      id = id_i, time = dose_times, evid = 1, amt = per_dose, cmt = "depot"
    )
    obs <- tibble::tibble(
      id = id_i, time = obs_times, evid = 0, amt = NA_real_, cmt = "central"
    )
    dplyr::bind_rows(doses, obs)
  }
  rows <- dplyr::bind_rows(lapply(seq_len(n), per_subject))
  rows$WT             <- wt
  rows$CYP3A5_EXPR    <- cyp3a5_expr
  rows$DONOR_DECEASED <- donor_deceased
  rows$CRCL           <- egfr
  rows$HCT            <- hct
  rows$IMMUNOASSAY    <- immuno
  rows$regimen        <- regimen
  rows
}
```

# Simulation

The model is loaded from the registry. For the dose-effect replications
below we use typical-value simulations (zero-IIV) to reproduce Figure 3,
since each Figure-3 panel varies a single covariate while holding all
others at population-typical values.

```{r simulate-typical}
mod <- readModelDb("Andrews_2017_tacrolimus")
mod_typical <- mod |> rxode2::zeroRe()
```

# Replicate published figures

Andrews 2017 Figure 3 shows the simulated steady-state tacrolimus
concentration-time profiles as each covariate is varied independently. Each
panel uses a 0.3 mg/kg/day starting dose split into two q12h doses,
simulated for 14 days, with all other covariates fixed to the cohort median.
The figures below reproduce Figure 3 panels (a), (b), (c), (d), and (e).

```{r figure-3a}
# Figure 3a: CYP3A5 nonexpresser (*3/*3) vs expresser (*1/*1 or *1/*3),
# all other covariates fixed to cohort median (WT 28.4 kg, deceased donor,
# eGFR 69, HCT 0.30, LC-MS/MS).
events_3a <- dplyr::bind_rows(
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = 69, hct = 0.30, immuno = 0,
              regimen = "CYP3A5 *3/*3 nonexpresser", id_offset = 0L),
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 1, donor_deceased = 1,
              egfr = 69, hct = 0.30, immuno = 0,
              regimen = "CYP3A5 *1 expresser",      id_offset = 100L)
)
sim_3a <- rxode2::rxSolve(mod_typical, events = events_3a,
                          keep = c("regimen")) |> as.data.frame()

ggplot(sim_3a |> dplyr::filter(time >= 12 * 13, time <= 12 * 14 + 12),
       aes(time - 12 * 13, Cc, colour = regimen)) +
  geom_line(size = 1) +
  labs(x = "Time after dose at steady state (h)",
       y = "Tacrolimus whole-blood Cc (ng/mL)",
       title = "Figure 3a: CYP3A5 genotype effect",
       caption = "Replicates Figure 3a of Andrews 2017.") +
  theme_minimal()
```

```{r figure-3b}
# Figure 3b: living vs deceased donor, all other covariates fixed.
events_3b <- dplyr::bind_rows(
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 0,
              egfr = 69, hct = 0.30, immuno = 0,
              regimen = "Living donor",   id_offset = 0L),
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = 69, hct = 0.30, immuno = 0,
              regimen = "Deceased donor", id_offset = 100L)
)
sim_3b <- rxode2::rxSolve(mod_typical, events = events_3b,
                          keep = c("regimen")) |> as.data.frame()

ggplot(sim_3b |> dplyr::filter(time >= 12 * 13, time <= 12 * 14 + 12),
       aes(time - 12 * 13, Cc, colour = regimen)) +
  geom_line(size = 1) +
  labs(x = "Time after dose at steady state (h)",
       y = "Tacrolimus whole-blood Cc (ng/mL)",
       title = "Figure 3b: Donor-type effect",
       caption = "Replicates Figure 3b of Andrews 2017.") +
  theme_minimal()
```

```{r figure-3c}
# Figure 3c: body weight 15, 30, 55, 80 kg sweep, all other covariates fixed.
wt_levels <- c(15, 30, 55, 80)
events_3c <- dplyr::bind_rows(lapply(seq_along(wt_levels), function(j) {
  wt_j <- wt_levels[j]
  make_cohort(n = 1, wt = wt_j, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = 69, hct = 0.30, immuno = 0,
              regimen = sprintf("WT %d kg", wt_j),
              id_offset = 100L * (j - 1L))
}))
sim_3c <- rxode2::rxSolve(mod_typical, events = events_3c,
                          keep = c("regimen")) |> as.data.frame()

ggplot(sim_3c |> dplyr::filter(time >= 12 * 13, time <= 12 * 14 + 12),
       aes(time - 12 * 13, Cc, colour = regimen)) +
  geom_line(size = 1) +
  labs(x = "Time after dose at steady state (h)",
       y = "Tacrolimus whole-blood Cc (ng/mL)",
       title = "Figure 3c: Body-weight effect (constant mg/kg/day)",
       caption = "Replicates Figure 3c of Andrews 2017.") +
  theme_minimal()
```

```{r figure-3d}
# Figure 3d: hematocrit 0.20, 0.25, 0.30 L/L sweep.
hct_levels <- c(0.20, 0.25, 0.30)
events_3d <- dplyr::bind_rows(lapply(seq_along(hct_levels), function(j) {
  hct_j <- hct_levels[j]
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = 69, hct = hct_j, immuno = 0,
              regimen = sprintf("HCT %.2f L/L", hct_j),
              id_offset = 100L * (j - 1L))
}))
sim_3d <- rxode2::rxSolve(mod_typical, events = events_3d,
                          keep = c("regimen")) |> as.data.frame()

ggplot(sim_3d |> dplyr::filter(time >= 12 * 13, time <= 12 * 14 + 12),
       aes(time - 12 * 13, Cc, colour = regimen)) +
  geom_line(size = 1) +
  labs(x = "Time after dose at steady state (h)",
       y = "Tacrolimus whole-blood Cc (ng/mL)",
       title = "Figure 3d: Hematocrit effect (piecewise; active when HCT < 0.3 L/L)",
       caption = "Replicates Figure 3d of Andrews 2017.") +
  theme_minimal()
```

```{r figure-3e}
# Figure 3e: eGFR 10, 50, 90 mL/min/1.73 m^2 sweep.
egfr_levels <- c(10, 50, 90)
events_3e <- dplyr::bind_rows(lapply(seq_along(egfr_levels), function(j) {
  egfr_j <- egfr_levels[j]
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = egfr_j, hct = 0.30, immuno = 0,
              regimen = sprintf("eGFR %d mL/min/1.73 m^2", egfr_j),
              id_offset = 100L * (j - 1L))
}))
sim_3e <- rxode2::rxSolve(mod_typical, events = events_3e,
                          keep = c("regimen")) |> as.data.frame()

ggplot(sim_3e |> dplyr::filter(time >= 12 * 13, time <= 12 * 14 + 12),
       aes(time - 12 * 13, Cc, colour = regimen)) +
  geom_line(size = 1) +
  labs(x = "Time after dose at steady state (h)",
       y = "Tacrolimus whole-blood Cc (ng/mL)",
       title = "Figure 3e: eGFR effect",
       caption = "Replicates Figure 3e of Andrews 2017.") +
  theme_minimal()
```

# PKNCA validation

Andrews 2017 Section 3.5 maps target trough concentration C0 to AUC0-12h
using the formula `dose = AUC * CL/F`. The mapping is:

| C0 (ng/mL) | AUC0-12h (ng*h/mL) |
|------------|--------------------|
| 10         | 177                |
| 12.5       | 209                |
| 15         | 241                |
| 17.5       | 274                |
| 20         | 306                |

We verify this by computing AUC0-12h from the simulated steady-state
12-hour interval for a reference subject and comparing the C0 / AUC ratio
against the published mapping. We use the same typical-value reference
subject (WT 28.4 kg, CYP3A5 nonexpresser, deceased donor, eGFR 69, HCT 0.30,
LC-MS/MS) but extend the dosing to a steady-state interval suitable for
the PKNCA calculation.

```{r pknca}
# Build a cohort spanning a range of doses for the typical reference subject;
# the resulting AUC0-12h values at steady state can be plotted against
# observed C0 to compare with Table 4 of Andrews 2017.
doses_mg_per_kg <- c(0.10, 0.20, 0.30, 0.45, 0.60, 0.80, 1.00, 1.20)
events_nca <- dplyr::bind_rows(lapply(seq_along(doses_mg_per_kg), function(j) {
  d_j <- doses_mg_per_kg[j]
  make_cohort(n = 1, wt = 28.4, cyp3a5_expr = 0, donor_deceased = 1,
              egfr = 69, hct = 0.30, immuno = 0,
              dose_mg_per_kg = d_j,
              regimen = sprintf("dose %.2f mg/kg/day", d_j),
              id_offset = 100L * (j - 1L))
}))
sim_nca_raw <- rxode2::rxSolve(mod_typical, events = events_nca,
                               keep = c("regimen")) |> as.data.frame()

# Restrict to the final 12-h dosing interval at steady state (days 13-14)
# and re-base time at zero for PKNCA.
sim_nca <- sim_nca_raw |>
  dplyr::filter(time >= 12 * 26, time <= 12 * 27) |>
  dplyr::mutate(time = time - 12 * 26) |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, regimen)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | regimen + id)

dose_df <- events_nca |>
  dplyr::filter(evid == 1) |>
  dplyr::group_by(id) |>
  dplyr::slice_tail(n = 1) |>
  dplyr::ungroup() |>
  dplyr::mutate(time = 0) |>
  dplyr::select(id, time, amt, regimen)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | regimen + id)

intervals <- data.frame(start = 0, end = 12,
                        cmax = TRUE, tmax = TRUE, auclast = TRUE)
nca_data  <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res   <- PKNCA::pk.nca(nca_data)
nca_df    <- as.data.frame(nca_res$result)

knitr::kable(nca_df, caption = "Simulated NCA parameters per dose level (PKNCA).")
```

## Comparison against published AUC0-12h vs C0 mapping

```{r compare-auc-c0}
# Extract observed pre-dose C0 (at time 0 of the steady-state interval) and
# AUClast for each regimen, then compare against the published Section 3.5
# mapping.
pred_c0  <- sim_nca |>
  dplyr::group_by(regimen) |>
  dplyr::filter(time == min(time)) |>
  dplyr::summarise(c0_sim = mean(Cc), .groups = "drop")

pred_auc <- nca_df |>
  dplyr::filter(PPTESTCD == "auclast") |>
  dplyr::group_by(regimen) |>
  dplyr::summarise(auc_sim = mean(PPORRES), .groups = "drop")

published_mapping <- tibble::tribble(
  ~c0_pub,   ~auc_pub,
  10.0,      177,
  12.5,      209,
  15.0,      241,
  17.5,      274,
  20.0,      306
)

cmp <- dplyr::full_join(pred_c0, pred_auc, by = "regimen") |>
  dplyr::arrange(c0_sim) |>
  dplyr::mutate(
    auc_pub_interp = approx(x = published_mapping$c0_pub,
                            y = published_mapping$auc_pub,
                            xout = c0_sim, rule = 2)$y,
    pct_diff = 100 * (auc_sim - auc_pub_interp) / auc_pub_interp
  )

knitr::kable(cmp, digits = 2,
             caption = "Simulated AUC0-12h vs Andrews 2017 Section 3.5 mapping (linear interpolation of the published C0 -> AUC anchors).")
```

The simulated AUC0-12h values track the published `dose = AUC * CL/F`
mapping linearly through the published range, confirming the reference
subject's typical CL/F of 50.5 L/h reproduces the dose-AUC-C0 relationship
that the paper uses to derive Table 4.

# Assumptions and deviations

- **Inter-occasion variability (IOV) not encoded.** Andrews 2017 Table 2
  reports IOV on CL/F (18% CV) and V2/F (35% CV) in addition to the diagonal
  IIV. This model file does NOT encode IOV structurally, following the
  Brooks 2021 tacrolimus precedent: the source paper does not define an
  operational occasion column suitable for the simulation-only model-library
  use case. Downstream users who want to simulate IOV can add an OCC
  indicator column to the event dataset and a per-occasion eta in rxode2.

- **CYP3A5 unknown stratum pooled with *3/*3 nonexpressers.** The model uses
  a single binary `CYP3A5_EXPR` indicator (0 = *3/*3 OR genotype unknown;
  1 = *1/*1 or *1/*3 carrier) following the published equation, which
  collapses *3/*3 and unknown into one stratum. Datasets with a substantial
  fraction of unknown-genotype subjects (Andrews 2017 model-building cohort
  had 41% unknown) are simulated as if those subjects were *3/*3
  nonexpressers, which is the assumption Andrews 2017 made in their
  starting-dose dosing guideline (Table 4).

- **CYP3A5 multiplier 1.04 / 1.98 vs published equation 1.0 / 2.0.** Table 2
  reports the *3/*3-or-unknown CL/F multiplier as 1.04 and the *1-carrier
  multiplier as 1.98; the published Section 3.2 equation rounds these to 1.0
  and 2.0. This model file uses the Table 2 final-model estimates as the
  reproducible parameter values. The relative ratio (~1.90-fold higher CL/F
  in expressers) is essentially identical between the rounded and unrounded
  forms.

- **CYP3A4 genotype not modelled.** Andrews 2017 Section 4 / Discussion
  notes that CYP3A4 *22 carriers were absent from the cohort (0 of 16
  genotyped subjects), and CYP3A4 *1G genotype was tested as a covariate on
  V2/F (univariate DOFV 5.0) but was not retained in the final model after
  backward elimination. The model file therefore does not include any
  CYP3A4-related effect.

- **HCT cutoff is sharp at 0.3 L/L.** Andrews 2017 implements the hematocrit
  effect as a piecewise power form: `f_HCT = (HCT / 0.3)^(-0.44)` when
  HCT < 0.3 L/L, else `f_HCT = 1`. This is continuous at HCT = 0.3 (both
  branches give 1.0) but the first derivative has a kink there. Users
  simulating around the 0.3 boundary should be aware of this kink.

- **Reference subject for the structural typical-value CL/F.** The Table 2
  CL/F = 50.5 L/h corresponds to a hypothetical reference paediatric subject
  with WT 70 kg (adult-equivalent), CYP3A5 *3/*3-or-unknown (multiplier
  1.04, the Table 2 final-model value), deceased-donor graft (multiplier
  1.0), eGFR 69 mL/min/1.73 m^2 (cohort median), and HCT >= 0.3 L/L (the
  piecewise hematocrit effect is inactive at the reference).

- **Hematocrit reported as fraction (L/L), NOT percent.** Andrews 2017 uses
  HCT as a fraction (0-1, L/L). The canonical-register HCT entry units
  default to percent (0-100); the per-model `covariateData[[HCT]]$units`
  field is set to "fraction" here to match the source. Datasets with HCT
  recorded as percent must be multiplied by 0.01 before passing to this
  model.

- **eGFR derived from the adapted-Schwartz paediatric formula.** Andrews
  2017 Section 2.2 specifies `eGFR = K * height(cm) / serum_creatinine(umol/L)`
  with `K = 36.5` (the paediatric K constant). The reference value of
  69 mL/min/1.73 m^2 is the cohort median. The canonical CRCL covariate
  accepts BSA-normalised eGFR regardless of the derivation formula; the
  derivation method is documented in the per-model notes.

- **Applicability window.** The model is intended for the FIRST 6 WEEKS
  post-transplantation with twice-daily IMMEDIATE-RELEASE oral tacrolimus
  in children aged 2-18 years. It is NOT validated for the once-daily
  extended-release formulation, for adults, or for the stable maintenance
  phase. Andrews 2017 Section 4 / Discussion explicitly limits applicability.

- **Race / ethnicity covariate not retained.** Andrews 2017 univariate
  analysis flagged ethnicity (DOFV 5.3 on CL/F; coefficients 0.84 / 0.95 /
  1.15 for Caucasian / Asian / Black vs Other) but the covariate was not
  retained after backward elimination (Table 3). The model file therefore
  does not include any race / ethnicity effect.
