---
title: "Tacrolimus (Lu 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tacrolimus (Lu 2015)}
  %\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::rxode2(readModelDb("Lu_2015_tacrolimus"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Lu_2015_tacrolimus"))$description`
* Article: <https://doi.org/10.1038/aps.2014.110>

# Population

Lu 2015 pooled two adult Chinese cohorts in a single popPK analysis (Table 1):

* **Liver-transplant recipients** -- 112 adults (86 male, 26 female; mean age
  58.4 +/- 11.6 years, range 28 - 78; 36 subjects >= 65 years) enrolled at the
  Organ Transplant Center of the General Hospital of the Armed Police Forces
  from 2004 - 2012. Underlying liver disease: primary hepatic carcinoma (50.9%),
  liver cirrhosis or chronic severe hepatitis (30.6%), hepatic cancer recurrence
  (6.5%), alcoholic cirrhosis (4.6%), other. All received the triple
  immunosuppressive regimen of tacrolimus + mycophenolate mofetil 750 mg bid +
  corticosteroid taper. Tacrolimus was initiated at 0.05 mg/kg/day in two divided
  oral doses and titrated by therapeutic drug monitoring (TDM); 1100 trough
  whole-blood samples were collected over POD 2 - 137 (mean 19.4 +/- 17.8 days)
  and measured by MEIA (linear range 1.5 - 30 ug/L).

* **Healthy Chinese volunteers** -- 40 adult males (mean age 28.7 +/- 3.47
  years, range 24 - 37) recruited at the PLA Second Artillery General Hospital
  for a randomised two-formulation crossover bioequivalence study. Each received
  a single 2 mg oral dose of tacrolimus; 851 whole-blood samples were drawn at
  0, 0.5, 1, 1.5, 2, 3, 4, 6, 12, 24, 48, and 72 h after dose and measured by
  HPLC-MS (linear range 0.1 - 25 ug/L).

The pooled analysis dataset comprises 152 subjects and 1951 observations.
Baseline laboratory values differ markedly between cohorts: alanine
aminotransferase (ALT) in transplant patients averaged 146 +/- 290 IU/L vs.
33 +/- 21 IU/L in healthy volunteers, motivating the ALT covariate effect on
CL/F (Table 2). Hematocrit, hemoglobin, weight, and bilirubin were also lower
in patients but were not retained in the final model (Table 3).

The population metadata is available programmatically:

```{r population}
str(rxode2::rxode2(readModelDb("Lu_2015_tacrolimus"))$population)
```

# Source trace

The per-parameter origin is also recorded as an in-file comment in
`inst/modeldb/specificDrugs/Lu_2015_tacrolimus.R`.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| Two-compartment oral PK with first-order absorption and lag (`d/dt(depot)`, `d/dt(central)`, `d/dt(peripheral1)`, `alag(depot)`) | n/a | Methods "Population pharmacokinetic model development" (paper p.283; NONMEM ADVAN4) |
| `CL/F` typical (patient at ALT = 0) `= theta1 * theta7 = 18.4 L/h` | `lcl <- log(32.8 * 0.562)` | Table 2 final-model column: CL/F = 32.8 L/h x CL_SubPop = 0.562 |
| `V2/F = 22.7 L` | `lvc <- log(22.7)` | Table 2 final-model column |
| `Q/F = 76.3 L/h` | `lq <- log(76.3)` | Table 2 final-model column |
| `V3/F = 916 L` (fixed) | `lvp <- fixed(log(916))` | Table 2 final-model column (asterisk = fixed; Methods text: "V3/F was fixed as the estimate from the healthy Chinese volunteer group") |
| `Ka = 0.419 1/h` | `lka <- log(0.419)` | Table 2 final-model column |
| `ALAG1 = 0.404 h` | `ltlag <- log(0.404)` | Table 2 final-model column |
| `CL_SubPop = theta7 = 0.562` (multiplicative ratio patient/healthy) | `e_dis_healthy_cl <- 1 / 0.562 = 1.78` | Table 2 final-model column; re-expressed as healthy/patient ratio for the canonical DIS_HEALTHY = 1 orientation |
| `CL_ALT = theta8 = -0.0237` (ALT exponent in patient CL/F) | `e_alt_cl <- -0.0237` | Table 2 final-model column |
| `Eq. 10`: `CL/F = theta1 * theta7^SubPop * Exp(ALT/40 * theta8)^SubPop * eta` | re-expressed with SubPop = 1 - DIS_HEALTHY | Paper p.284 Eq. 10 (final-model regression) |
| `Eq. 11 - 15`: `V2/F`, `Q/F`, `V3/F`, `Ka`, `ALAG1` with multiplicative eta | identity-mapped to `exp(l<param> + etal<param>)` | Paper p.284 Eqs. 11 - 15 |
| `CL-IIV = 46.6% CV` | `etalcl ~ 0.1965` (= log(0.466^2 + 1)) | Table 2 final-model column (Methods: IIV expressed exponentially) |
| `V2-IIV = 57.3% CV` | `etalvc ~ 0.2841` | Table 2 final-model column |
| `Q-IIV = 46.0% CV` | `etalq  ~ 0.1919` | Table 2 final-model column |
| `V3-IIV = 93.5% CV` | `etalvp ~ 0.6280` | Table 2 final-model column |
| `Ka-IIV = 0% (fixed)` | no eta on ka | Table 2 final-model column (asterisk = fixed) |
| `Proportional residual = 39.8%` | `propSd <- 0.398` | Table 2 final-model column (Methods: residual expressed proportionally) |

# Virtual cohort

Original observed concentrations are not publicly available. The chunks below
build two virtual cohorts that mirror the published study design:

* Cohort A -- 40 healthy male volunteers receiving a single 2 mg oral dose
  with dense sampling out to 72 h (Lu 2015 Methods "Dosage regimen" and
  "Blood sample collecting and analysis").
* Cohort B -- 50 liver-transplant patients receiving 1.5 mg orally q12h
  (a typical TDM-adjusted regimen for a 60 kg patient at 0.05 mg/kg/day) for
  7 days, with trough samples just before each morning dose.

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

make_cohort <- function(n, dis_healthy, alt_mean, alt_sd, dose_amt, dose_times,
                        obs_times, id_offset = 0L) {
  alt_vals <- pmax(1, rnorm(n, mean = alt_mean, sd = alt_sd))
  cohort_label <- if (dis_healthy == 1L) "Healthy volunteer (single 2 mg PO)" else "Liver-transplant patient (1.5 mg PO q12h)"

  dose_rows <- expand.grid(id = id_offset + seq_len(n), time = dose_times) |>
    dplyr::mutate(evid = 1L, amt = dose_amt, cmt = "depot")
  obs_rows <- expand.grid(id = id_offset + seq_len(n), time = obs_times) |>
    dplyr::mutate(evid = 0L, amt = 0,        cmt = "central")
  rows <- dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid)) |>
    dplyr::mutate(
      DIS_HEALTHY = dis_healthy,
      ALT         = alt_vals[match(id, id_offset + seq_len(n))],
      cohort      = cohort_label
    )
  rows
}

# Cohort A: healthy volunteers, single 2 mg PO dose
obs_grid_hv <- c(0.1, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 6, 8, 10, 12, 18, 24, 36, 48, 60, 72)
events_hv <- make_cohort(
  n = 40, dis_healthy = 1L,
  alt_mean = 33.0, alt_sd = 20.9,
  dose_amt = 2, dose_times = 0,
  obs_times = obs_grid_hv,
  id_offset = 0L
)

# Cohort B: liver-transplant patients, 1.5 mg PO q12h for 7 days; trough samples
# just before each morning dose.
dose_times_tx <- seq(0, 24 * 7 - 12, by = 12)
obs_grid_tx   <- seq(24, 24 * 7, by = 24) - 0.1  # ~10 minutes before the morning dose
events_tx <- make_cohort(
  n = 50, dis_healthy = 0L,
  alt_mean = 146.4, alt_sd = 50,        # heavy-tailed in the paper; truncate at sd = 50 to avoid extreme outliers
  dose_amt = 1.5, dose_times = dose_times_tx,
  obs_times = obs_grid_tx,
  id_offset = 200L
)

events <- dplyr::bind_rows(events_hv, events_tx)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode2(readModelDb("Lu_2015_tacrolimus"))
sim <- rxode2::rxSolve(mod, events = events, keep = c("cohort", "DIS_HEALTHY", "ALT"))
sim_df <- as.data.frame(sim)
```

A deterministic (zero-random-effects) typical-value run is useful for
comparing the structural model against published mean curves:

```{r simulate-typical}
mod_typ <- rxode2::zeroRe(mod)
typ_hv  <- rxode2::rxSolve(mod_typ, events = events_hv[events_hv$id == 1L, ])
```

# Replicate published figures

Lu 2015 Figure 2A shows the visual predictive check for the healthy-volunteer
cohort (concentration vs. time after a single 2 mg oral dose, 5th / 50th / 95th
percentiles of the simulation). The chunk below reproduces the same plot from
the packaged model.

```{r figure-2A}
sim_hv <- sim_df |> dplyr::filter(DIS_HEALTHY == 1L, Cc > 0)

sim_hv_summary <- sim_hv |>
  dplyr::group_by(time) |>
  dplyr::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"
  )

ggplot(sim_hv_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue4", linewidth = 1) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)",
    y = "Whole-blood tacrolimus (ng/mL)",
    title = "Figure 2A -- healthy-volunteer VPC after a single 2 mg oral dose",
    caption = "Replicates Lu 2015 Figure 2A. Shaded band = 5th - 95th percentile of the simulated cohort."
  )
```

Lu 2015 Figure 2B reports the patient-cohort VPC over the 0 - 1500 h
post-transplant window. The dosing reconstruction below is illustrative
rather than a faithful reproduction (per-patient TDM-adjusted doses are not
published) so the figure here uses the simpler 1.5 mg q12h regimen of cohort B.

```{r figure-2B}
sim_tx <- sim_df |> dplyr::filter(DIS_HEALTHY == 0L, Cc > 0)

sim_tx_summary <- sim_tx |>
  dplyr::group_by(time) |>
  dplyr::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"
  )

ggplot(sim_tx_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "firebrick") +
  geom_line(colour = "firebrick4", linewidth = 1) +
  scale_y_log10() +
  labs(
    x = "Time (h)",
    y = "Whole-blood tacrolimus trough (ng/mL)",
    title = "Figure 2B-like -- liver-transplant trough simulation (1.5 mg PO q12h)",
    caption = "Approximate reconstruction of Lu 2015 Figure 2B; the published figure uses per-patient TDM-adjusted doses that are not tabulated."
  )
```

# PKNCA validation

Single-dose NCA on the healthy-volunteer cohort gives a direct check that the
simulated Cmax, Tmax, AUC, and half-life are in the range expected for a 2 mg
oral tacrolimus dose in healthy adults.

```{r pknca-hv}
sim_nca_hv <- sim_hv |> dplyr::select(id, time, Cc, cohort)

dose_df_hv <- events_hv |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(cohort = "Healthy volunteer (single 2 mg PO)") |>
  dplyr::select(id, time, amt, cohort)

conc_obj_hv <- PKNCA::PKNCAconc(
  sim_nca_hv,
  Cc ~ time | cohort + id,
  concu = "ng/mL",
  timeu = "h"
)
dose_obj_hv <- PKNCA::PKNCAdose(
  dose_df_hv,
  amt ~ time | cohort + id,
  doseu = "mg"
)

intervals_hv <- data.frame(
  start       = 0,
  end         = 72,
  cmax        = TRUE,
  tmax        = TRUE,
  auclast     = TRUE,
  aucinf.obs  = TRUE,
  half.life   = TRUE,
  clast.obs   = TRUE
)

nca_hv <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_hv, dose_obj_hv, intervals = intervals_hv)
)
knitr::kable(
  summary(nca_hv),
  caption = "Healthy-volunteer simulated NCA after a single 2 mg PO tacrolimus dose."
)
```

For the steady-state patient cohort, AUC over the final 12 h dosing interval
gives Cmax,ss / Cmin,ss / Cavg,ss:

```{r pknca-tx-ss}
sim_nca_tx <- sim_df |>
  dplyr::filter(DIS_HEALTHY == 0L) |>
  dplyr::select(id, time, Cc, cohort)

dose_df_tx <- events_tx |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(cohort = "Liver-transplant patient (1.5 mg PO q12h)") |>
  dplyr::select(id, time, amt, cohort)

conc_obj_tx <- PKNCA::PKNCAconc(
  sim_nca_tx,
  Cc ~ time | cohort + id,
  concu = "ng/mL",
  timeu = "h"
)
dose_obj_tx <- PKNCA::PKNCAdose(
  dose_df_tx,
  amt ~ time | cohort + id,
  doseu = "mg"
)

intervals_tx <- data.frame(
  start    = max(dose_df_tx$time),
  end      = max(dose_df_tx$time) + 12,
  cmin     = TRUE,
  auclast  = TRUE,
  cav      = TRUE
)

nca_tx <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_obj_tx, dose_obj_tx, intervals = intervals_tx)
)
knitr::kable(
  summary(nca_tx),
  caption = "Liver-transplant simulated steady-state NCA over the final 12 h dosing interval (1.5 mg PO q12h)."
)
```

## Comparison against published NCA

Lu 2015 does not tabulate per-cohort NCA values (Cmax, AUC, half-life are not
reported in the Results text or any Table). The healthy-volunteer trough
concentrations in Lu 2015 Figure 2A range from ~ 0.5 - 30 ng/mL across the
72 h sampling window, with peak around 1 - 2 h after dose. The simulated
percentile band reproduces that envelope. For the transplant cohort, the
typical TDM target trough is 5 - 15 ng/mL in the early post-transplant period
and the simulated steady-state troughs at 1.5 mg q12h fall within that range.

# Assumptions and deviations

* **Re-parameterised covariate orientation.** Lu 2015 Eq. 10 writes the
  multiplicative cohort effect on CL/F using the source's `SubPop` indicator
  (SubPop = 0 for healthy volunteer, SubPop = 1 for liver-transplant
  patient). The canonical covariate-columns register orients `DIS_HEALTHY = 1`
  for healthy and `DIS_HEALTHY = 0` for patient (the reference category for
  the canonical column), with the identity `SubPop = 1 - DIS_HEALTHY`. The
  model file re-parameterises the typical CL/F as the patient typical
  (`lcl = log(32.8 * 0.562) = log(18.4)`, i.e. the paper's reported
  "patient CL/F = 18.4 L/h" at ALT = 0) and introduces the multiplicative
  factor `e_dis_healthy_cl = 1 / 0.562 = 1.78` for the healthy direction.
  The ALT exponent `e_alt_cl = -0.0237` is gated on `(1 - DIS_HEALTHY)` so it
  contributes only in the transplant cohort. This is mathematically identical
  to Lu 2015 Eq. 10; only the orientation of the typical-vs-deviation parameters
  changes.
* **ALT normalisation factor.** Lu 2015 Eq. 10 normalises ALT by 40 (IU/L)
  without stating the source of 40. We interpret it as the clinical upper
  limit of normal for serum ALT used by the central laboratory (commonly
  40 IU/L for adults); the paper's reported "typical patient CL/F = 18.4 L/h"
  in the abstract is consistent with the value of `theta1 * theta7` at ALT = 0
  (where the exponential ALT term equals 1).
* **Patient dosing reconstruction.** Per-patient TDM-adjusted dose histories
  are not published. The transplant-cohort simulation uses a uniform 1.5 mg
  q12h regimen as a representative case for a 60 kg patient at the protocol
  starting dose of 0.05 mg/kg/day, not a one-to-one VPC reproduction of
  Lu 2015 Figure 2B.
* **Patient ALT distribution.** Lu 2015 Table 1 reports patient ALT as
  146.4 +/- 290.0 IU/L (range 5 - 6300). The empirical distribution is heavily
  right-skewed by a small number of acute-rejection episodes. The simulation
  uses a truncated normal (mean 146.4, SD 50, lower bound 1 IU/L) to keep
  ALT in a clinically plausible range; this slightly under-disperses the CL/F
  spread relative to the actual cohort.
* **Healthy-volunteer ALT.** ALT enters the model only via the
  `(1 - DIS_HEALTHY)` gate, so the ALT value supplied for healthy volunteers
  has no effect on CL/F. We populate it with the cohort empirical mean
  (33 IU/L, SD 20.9) for the sake of dataset completeness.
* **No erratum.** A PubMed search for "Lu 2015 tacrolimus erratum" and a
  publisher search of the Acta Pharmacologica Sinica article landing page
  returned no corrigendum or correction notice.
