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

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Kim_2018_tacrolimus"))$reference`
* Description: `r rxode2::rxode(readModelDb("Kim_2018_tacrolimus"))$description`
* Article: <https://doi.org/10.1038/s41598-018-20071-3>

This is the **final integrated model** of Kim et al. (2018), which jointly
describes tacrolimus (TAC) and the mycophenolate mofetil (MMF) metabolite
cascade -- mycophenolic acid (MPA), MPA 7-O-glucuronide (MPAG), and MPA acyl
glucuronide (AcMPAG) -- and quantifies the TAC-MMF drug-drug interaction. The
independent (TAC-alone / MMF-alone) models in Table 3 are intermediate
development steps and are not packaged separately.

# Population

The model was developed from a three-period, fixed-sequence, open-label,
single-dose drug-drug-interaction study in 17 healthy Korean male volunteers
(NCT02743247; conducted 2015-2016). Subjects received 1,000 mg MMF alone, then
5 mg TAC alone, then the 1,000 mg MMF + 5 mg TAC combination, with one-week
washout periods. Median age was 25 years (range 20-42), median weight 69.7 kg
(57.4-88.3 kg), and all subjects were of Korean ethnicity with normal renal and
hepatic function (Table 1). Four of 17 subjects (23.5%) were CYP3A5 expressers
(CYP3A5*1/*1 or CYP3A5*1/*3). Because baseline demographic variation was small,
the CYP3A5 genotype was the only covariate retained in the final model.

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

# Source trace

Per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Kim_2018_tacrolimus.R`. All structural
estimates are the **Integrated model** column of Kim 2018 Table 3.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` (Ka TAC) | 1.78 1/h | Table 3, Ka TAC (integrated) |
| `lcl` (CL/F TAC, typical) | 13.8 L/h | Table 3, CL/F TAC (integrated) |
| `lvc` (V2/F TAC) | 93 L | Table 3, V2/F (integrated) |
| `lk23`, `lk32` (TAC distribution) | 0.313, 0.0719 1/h | Table 3, k23 / k32 (integrated) |
| `ltlag` (TAC lag time) | 0.59 h | Table 3, Lag time (integrated) |
| `e_cyp3a5_cl` (CYP3A5 on CL/F) | 1.48 | Table 3 + equation (2) |
| `e_cmpa_cl` (MPA-TAC interaction slope) | 0.0294 mL/ug | Table 3, Interaction + equation (2) |
| `lka_mpa` (Ka MPA) | 2.29 1/h | Table 3, Ka MPA (integrated) |
| `lcl_mpa` (CL/F MPA) | 16.3 L/h | Table 3, CL/F MPA (integrated) |
| `lvc_mpa` (V5/F MPA) | 19.7 L | Table 3, V5/F (integrated) |
| `lk56`, `lk65` (MPA distribution) | 1.12, 0.131 1/h | Table 3, k56 / k65 (integrated) |
| `f_mpag` (fraction MPA -> MPAG, fixed) | 0.85 | Table 3, fMPA (fixed) |
| `lvc_mpag` (V7/F MPAG) | 5.83 L | Table 3, V7/F (integrated) |
| `lk70` (MPAG elimination) | 0.251 1/h | Table 3, k70 (integrated) |
| `lehc` (fraction MPAG EHC) | 0.367 | Table 3, EHC + equation (1) |
| `lk84` (gallbladder emptying) | 18.4 1/h | Table 3, k84 (integrated) |
| `lmtime1` (meal time) | 7.96 h | Table 3, MTIME1 (integrated) |
| `lmtime2` (emptying duration, fixed) | 1.0 h | Table 3, MTIME2 (fixed) |
| `lvc_acmpag` (V9/F AcMPAG, fixed) | 23 L | Table 3, V9/F (fixed) |
| `lk90` (AcMPAG elimination, fixed) | 2.15 1/h | Table 3, k90 (fixed) |
| Residual error (TAC / MPA / MPAG / AcMPAG) | see `ini()` | Table 3, sigma rows (integrated) |
| IIV (Ka/CL/V TAC; Ka/CL/V MPA; EHC) | see `ini()` | Table 3, IIV CV% (integrated) |
| `d/dt(...)` 9-compartment structure | n/a | Figure 1 (model schematic) |
| `k57 = fMPA*CL_MPA/V5`; `k59 = (1-fMPA)*CL_MPA/V5` | derived | Figure 1 (MPA exits only via k57/k59) |
| `k78 = EHC*k70/(1-EHC)` | derived | equation (1) |

The metabolite-formation rate constants `k57` (MPA -> MPAG) and `k59`
(MPA -> AcMPAG) and the biliary-transfer constant `k78` (MPAG -> gallbladder)
are not tabulated in Table 3. They are recovered algebraically from the
Figure 1 structure: MPA central is eliminated solely through metabolism, so
`k57 + k59 = CL/F_MPA / V5/F`, split by the fixed fraction `fMPA = 0.85`; and
equation (1), `EHC = k78 / (k70 + k78)`, inverts to `k78 = EHC*k70/(1-EHC)`.

# Virtual cohort

Original observed data are not publicly available. The figures below use a
virtual population matching the published CYP3A5 expresser frequency (23.5%);
the model carries no continuous demographic covariate, so weight/age are not
needed for simulation. Two treatment arms reproduce the study's interaction
contrast: TAC alone (5 mg) and TAC + MMF (5 mg TAC + 1,000 mg MMF).

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

n_per_arm <- 250L
obs_times <- sort(unique(c(seq(0, 12, by = 0.25), seq(12, 72, by = 1))))

make_cohort <- function(n, arm, with_mmf, id_offset = 0L) {
  cyp <- stats::rbinom(n, 1L, 4 / 17)   # 23.5% CYP3A5 expressers (Table 1)
  dplyr::bind_rows(lapply(seq_len(n), function(i) {
    ev <- rxode2::et(amt = 5, cmt = "depot", time = 0)
    if (with_mmf) ev <- rxode2::et(ev, amt = 1000, cmt = "depot_mpa", time = 0)
    ev <- rxode2::et(ev, obs_times, cmt = "Cc")
    df <- as.data.frame(ev)
    df$id <- id_offset + i
    df$arm <- arm
    df$CYP3A5_EXPR <- cyp[i]
    df
  }))
}

events <- dplyr::bind_rows(
  make_cohort(n_per_arm, "TAC alone",  with_mmf = FALSE, id_offset = 0L),
  make_cohort(n_per_arm, "TAC + MMF",  with_mmf = TRUE,  id_offset = n_per_arm)
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Kim_2018_tacrolimus")

# Stochastic simulation (between-subject variability) for the VPC bands.
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("arm", "CYP3A5_EXPR"),
                       returnType = "data.frame", addDosing = FALSE,
                       atol = 1e-8, rtol = 1e-6)
```

For the typical-value interaction figure, the random effects are zeroed:

```{r simulate-typical}
mod_typical <- rxode2::zeroRe(mod)

typ_events <- dplyr::bind_rows(
  make_cohort(1L, "TAC alone", with_mmf = FALSE, id_offset = 0L),
  make_cohort(1L, "TAC + MMF", with_mmf = TRUE,  id_offset = 1L)
)
typ_events$CYP3A5_EXPR <- 0   # nonexpresser (majority genotype)

sim_typ <- rxode2::rxSolve(mod_typical, events = typ_events,
                           keep = c("arm"), returnType = "data.frame",
                           addDosing = FALSE)
```

# Replicate published figures

## Figure 2 -- VPC of the four analytes (combination arm)

Kim 2018 Figure 2 shows prediction-corrected VPCs of TAC, MPA, MPAG and AcMPAG
in the combination treatment. Without the observed data we plot the simulated
median and 5th-95th percentile prediction interval for the `TAC + MMF` arm.

```{r figure-2}
combo <- dplyr::filter(sim, arm == "TAC + MMF")

vpc_long <- combo |>
  dplyr::select(time, TAC = Cc, MPA = Cc_mpa, MPAG = Cc_mpag, AcMPAG = Cc_acmpag) |>
  tidyr::pivot_longer(-time, names_to = "analyte", values_to = "conc") |>
  dplyr::mutate(analyte = factor(analyte, levels = c("TAC", "MPA", "MPAG", "AcMPAG")))

vpc_bands <- vpc_long |>
  dplyr::group_by(analyte, time) |>
  dplyr::summarise(
    Q05 = quantile(conc, 0.05, na.rm = TRUE),
    Q50 = quantile(conc, 0.50, na.rm = TRUE),
    Q95 = quantile(conc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_bands, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~analyte, scales = "free_y") +
  labs(x = "Time (h)",
       y = "Concentration (TAC: ng/mL; MPA/MPAG/AcMPAG: ug/mL)",
       title = "Figure 2 -- integrated-model VPC (combination arm)",
       caption = "Replicates Figure 2 of Kim 2018 (predicted median + 5th-95th percentile; observed data not public).") +
  theme_bw()
```

The MPA panel shows the characteristic secondary peak around 8 h driven by
the meal-triggered gallbladder emptying (MTIME1 = 7.96 h) that completes the
enterohepatic-recirculation loop -- matching the "secondary peak at 6-12 h"
described in the paper.

## Tacrolimus interaction (typical value)

The integrated model encodes the MMF interaction as an inverse-exponential
suppression of TAC CL/F by MPA concentration (equation (2)). The typical-value
profiles below show the higher TAC exposure when co-administered with MMF.

```{r figure-interaction}
ggplot(sim_typ, aes(time, Cc, colour = arm)) +
  geom_line(linewidth = 0.9) +
  labs(x = "Time (h)", y = "Tacrolimus (ng/mL)",
       colour = NULL,
       title = "Tacrolimus typical-value profile: alone vs. with MMF",
       caption = "MMF co-administration lowers TAC CL/F (equation (2)), raising exposure.") +
  theme_bw()
```

# PKNCA validation

NCA parameters for tacrolimus are computed with PKNCA per treatment arm and
compared against Kim 2018 Table 2 (TAC alone vs. TAC + MMF).

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, arm)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | arm + id,
                             concu = "ng/mL", timeu = "hr")

dose_df <- events |>
  dplyr::filter(evid == 1, cmt == "depot") |>   # the tacrolimus dose
  dplyr::select(id, time, amt, arm)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id, doseu = "mg")

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))
knitr::kable(summary(nca_res),
             caption = "Simulated tacrolimus NCA by treatment arm.")
```

## Comparison against published NCA (Table 2)

```{r pknca-compare}
published <- tibble::tribble(
  ~arm,         ~cmax, ~tmax, ~aucinf.obs,
  "TAC alone",   33.5,   1,    322.4,
  "TAC + MMF",   35.0,   2,    393.6
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "arm",
  params        = c("cmax", "tmax", "aucinf.obs"),
  units         = c(cmax = "ng/mL", tmax = "h", aucinf.obs = "ng*h/mL"),
  tolerance_pct = 20
)
knitr::kable(
  cmp,
  caption = paste(
    "Tacrolimus NCA: simulated (median across n =",
    length(unique(sim_nca$id)),
    ") vs. published (Kim 2018 Table 2). * differs from reference by",
    ">20%."
  )
)
```

The Kim 2018 Table 2 also reports CL/F (TAC alone 20.1 L/h; TAC + MMF
14.2 L/h). The integrated model's typical CL/F is 13.8 L/h; the
TAC + MMF arm matches closely, while the TAC-alone arm runs lower than
the independent-model NCA-derived value (see Assumptions and deviations
below).

The simulated combination-arm exposure (`TAC + MMF`) closely matches the
published combination NCA, and the direction of the interaction (higher TAC
exposure / lower CL/F with MMF) is reproduced. The `TAC alone` simulated
exposure is modestly higher than the published TAC-alone NCA because the
**integrated** model's typical CL/F (13.8 L/h) was estimated jointly with the
interaction term and is lower than the independent TAC-alone estimate
(17.8 L/h) or the NCA-derived value (20.1 L/h); see the deviations below. No
parameters were tuned to match the published values.

# Assumptions and deviations

- **Integrated (final) model only.** The independent TAC-alone and MMF-alone
  models (Table 3, "Independent model" columns) are development steps and are
  not packaged. Simulating "TAC alone" with this model sets the MPA
  concentration to zero (no interaction); the resulting CL/F (13.8 L/h) is the
  integrated-model base value, which is lower than the independent-model
  TAC-alone CL/F (17.8 L/h). This is an intrinsic feature of the integrated
  parameterisation, not a transcription error, and explains the modest
  over-prediction of TAC-alone exposure relative to Table 2.

- **Derived metabolite / EHC rate constants.** `k57`, `k59` and `k78` are not
  tabulated; they are derived from the reported `CL/F_MPA`, `V5/F`, `fMPA`,
  `EHC` and `k70` together with the Figure 1 structure and equation (1) (see
  Source trace). MPA central is assumed to be eliminated solely via metabolism
  to MPAG and AcMPAG, consistent with Figure 1 (no separate elimination arrow
  from the MPA central compartment).

- **Dosing in apparent units.** TAC (5 mg) and MMF (1,000 mg) doses are entered
  as the administered amounts; the apparent clearance / volume parameters
  (CL/F, V/F) absorb both bioavailability and the MMF -> MPA conversion. TAC
  concentrations are reported in ng/mL (a factor of 1000 over the model's
  internal mg/L); MPA, MPAG and AcMPAG are in ug/mL.

- **Gallbladder emptying window.** The study is single-dose, so the meal-
  triggered emptying window is `[MTIME1, MTIME1 + MTIME2]` in absolute time.
  For multi-dose simulation a recurring (per-interval) window would be needed.

- **EHC fraction stochastic guard.** The source modelled IIV on the EHC
  fraction as exponential (log-normal, 35.5% CV), which can exceed 1 in the
  rare upper stochastic tail and would make `k78 = EHC*k70/(1-EHC)` diverge.
  The packaged model clips the per-subject EHC at 0.99 (`min(...)`) purely for
  numerical stability; the typical value (0.367) and the bulk of the
  distribution are unchanged, and roughly 0.5% of simulated subjects are
  affected.

- **Non-canonical parameter / compartment names.** The micro-rate-constant
  parameterisation (`lk23`, `lk32`, `lk56`, `lk65`, `lk70`, `lk84`, `lmtime1`,
  `lmtime2`, `lehc`, `f_mpag`) and the mechanistic compartment names
  (`gallbladder`, the `_mpa` / `_mpag` / `_acmpag` analyte suffixes) follow the
  paper's notation rather than the canonical CL/Q/V scheme, so
  `checkModelConventions()` reports naming warnings. These are intentional and
  preserve a direct line-by-line correspondence with Figure 1 and Table 3. The
  interaction parameter `e_cmpa_cl` is a structural drug-interaction slope
  driven by the model-internal MPA concentration, not a data covariate column.

- **Virtual cohort.** CYP3A5 expresser status was sampled at the published
  frequency (4/17 = 23.5%); no other covariate is used by the model. Observed
  concentration data are not public, so the Figure 2 reproduction shows the
  predicted median and 5th-95th percentile interval rather than a true
  prediction-corrected VPC.

- **Errata.** No published erratum or corrigendum to Kim 2018 was located. The
  supplementary information (Supplementary Tables S1-S3, Figures S1-S4) was not
  on disk; it contains genotype frequencies and MPA/MPAG/AcMPAG NCA values
  (Table S3) but no additional structural parameters, so its absence does not
  affect the packaged model.
