Skip to contents

Tacrolimus (Andrews 2017)

Source: vignettes/articles/Andrews_2017_tacrolimus.Rmd
Andrews_2017_tacrolimus.Rmd

Model and source

  • Citation: Andrews LM, Hesselink DA, van Gelder T, Koch BCP, Cornelissen EAM, Bruggemann RJM, van Schaik RHN, de Wildt SN, Cransberg K, de Winter BCM. A Population Pharmacokinetic Model to Predict the Individual Starting Dose of Tacrolimus Following Pediatric Renal Transplantation. Clin Pharmacokinet. 2018;57(4):475-489. doi:10.1007/s40262-017-0567-8 (published online 5 July 2017).
  • Description: Two-compartment population PK model with first-order absorption and an absorption lag time for twice-daily oral immediate-release tacrolimus (Prograft and Modigraf) in paediatric renal transplant recipients during the first 6 weeks post-transplantation (Andrews 2017). Apparent oral clearance CL/F and apparent inter-compartmental clearance Q/F scale allometrically with body weight at a fixed exponent of 0.75 referenced to a 70 kg adult; apparent central volume V1/F and apparent peripheral volume V2/F scale at a fixed exponent of 1.0; ka has no body-weight scaling. CL/F additionally varies with CYP3A5 expresser status (1.04 multiplier for 3/3 or unknown genotype, 1.98 multiplier for 1/1 or 1/3 carriers; pooled with unknown because Andrews 2017 explicitly groups 3/3 with unknown in the final equation), donor source (0.74 multiplier for living-donor recipients vs deceased-donor reference; equivalent to deceased-donor recipients having ~35% higher CL/F), eGFR (power exponent 0.19 centred at the cohort median 69 mL/min/1.73 m^2 of adapted-Schwartz eGFR), and a piecewise hematocrit effect (power exponent -0.44 centred at 0.3 L/L applied only when HCT < 0.3 L/L). Inter-individual variability is diagonal on ka, CL/F, V1/F, and V2/F. Residual error is a combined additive + proportional model with separate immunoassay and LC-MS/MS magnitudes selected by the per-sample IMMUNOASSAY indicator. Inter-occasion variability (IOV) on CL/F (18% CV) and V2/F (35% CV) reported by Andrews 2017 Table 2 is NOT encoded structurally here (per the Brooks 2021 tacrolimus precedent) – the source paper does not define an operational occasion column for the model-library use case; downstream users who want to simulate IOV can add an OCC indicator and a per-occasion eta in rxode2.
  • 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).

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.

mod <- readModelDb("Andrews_2017_tacrolimus")
mod_typical <- mod |> rxode2::zeroRe()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: No sigma parameters in the model

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).

# 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)
)
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_3a <- rxode2::rxSolve(mod_typical, events = events_3a,
                          keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

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()
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

# 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)
)
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_3b <- rxode2::rxSolve(mod_typical, events = events_3b,
                          keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

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()

# 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))
}))
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_3c <- rxode2::rxSolve(mod_typical, events = events_3c,
                          keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

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()

# 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))
}))
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_3d <- rxode2::rxSolve(mod_typical, events = events_3d,
                          keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

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()

# 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))
}))
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_3e <- rxode2::rxSolve(mod_typical, events = events_3e,
                          keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

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.

# 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))
}))
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
#> Warning in dose_times + c(0.5, 1, 2, 4, 6, 8): longer object length is not a
#> multiple of shorter object length
sim_nca_raw <- rxode2::rxSolve(mod_typical, events = events_nca,
                               keep = c("regimen")) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalka', 'etalcl', 'etalvc', 'etalvp'
#> Warning: multi-subject simulation without without 'omega'

# 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).")
Simulated NCA parameters per dose level (PKNCA).
regimen id start end PPTESTCD PPORRES exclude
dose 0.10 mg/kg/day 1 0 12 auclast 53.36178 NA
dose 0.10 mg/kg/day 1 0 12 cmax 7.28013 NA
dose 0.10 mg/kg/day 1 0 12 tmax 2.00000 NA
dose 0.20 mg/kg/day 101 0 12 auclast 106.72356 NA
dose 0.20 mg/kg/day 101 0 12 cmax 14.56027 NA
dose 0.20 mg/kg/day 101 0 12 tmax 2.00000 NA
dose 0.30 mg/kg/day 201 0 12 auclast 160.08533 NA
dose 0.30 mg/kg/day 201 0 12 cmax 21.84040 NA
dose 0.30 mg/kg/day 201 0 12 tmax 2.00000 NA
dose 0.45 mg/kg/day 301 0 12 auclast 240.12800 NA
dose 0.45 mg/kg/day 301 0 12 cmax 32.76060 NA
dose 0.45 mg/kg/day 301 0 12 tmax 2.00000 NA
dose 0.60 mg/kg/day 401 0 12 auclast 320.17067 NA
dose 0.60 mg/kg/day 401 0 12 cmax 43.68080 NA
dose 0.60 mg/kg/day 401 0 12 tmax 2.00000 NA
dose 0.80 mg/kg/day 501 0 12 auclast 426.89423 NA
dose 0.80 mg/kg/day 501 0 12 cmax 58.24106 NA
dose 0.80 mg/kg/day 501 0 12 tmax 2.00000 NA
dose 1.00 mg/kg/day 601 0 12 auclast 533.61778 NA
dose 1.00 mg/kg/day 601 0 12 cmax 72.80133 NA
dose 1.00 mg/kg/day 601 0 12 tmax 2.00000 NA
dose 1.20 mg/kg/day 701 0 12 auclast 640.34134 NA
dose 1.20 mg/kg/day 701 0 12 cmax 87.36159 NA
dose 1.20 mg/kg/day 701 0 12 tmax 2.00000 NA

Comparison against published AUC0-12h vs C0 mapping 

# 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).")
Simulated AUC0-12h vs Andrews 2017 Section 3.5 mapping (linear interpolation of the published C0 -> AUC anchors).
regimen c0_sim auc_sim auc_pub_interp pct_diff
dose 0.10 mg/kg/day 2.89 53.36 177.00 -69.85
dose 0.20 mg/kg/day 5.78 106.72 177.00 -39.70
dose 0.30 mg/kg/day 8.67 160.09 177.00 -9.56
dose 0.45 mg/kg/day 13.01 240.13 215.54 11.41
dose 0.60 mg/kg/day 17.35 320.17 272.00 17.71
dose 0.80 mg/kg/day 23.13 426.89 306.00 39.51
dose 1.00 mg/kg/day 28.91 533.62 306.00 74.38
dose 1.20 mg/kg/day 34.70 640.34 306.00 109.26

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.