Skip to contents

Deoxypodophyllotoxin interspecies PBPK (Chen 2016)

Source: vignettes/articles/Chen_2016_deoxypodophyllotoxin.Rmd
Chen_2016_deoxypodophyllotoxin.Rmd

Model and source

  • Citation: Chen Y, Zhao K, Liu F, Xie Q, Zhong Z, Miao M, Liu X, Liu L. Prediction of Deoxypodophyllotoxin Disposition in Mouse, Rat, Monkey, and Dog by Physiologically Based Pharmacokinetic Model and the Extrapolation to Human. Front Pharmacol. 2016;7:488. doi:10.3389/fphar.2016.00488

  • Article: Front Pharmacol. 2016;7:488 (open access)

Deoxypodophyllotoxin (DPT) is a pre-clinical anti-tumour candidate. Chen and colleagues built a whole-body physiologically based pharmacokinetic (PBPK) model of DPT from in vitro and in silico inputs, developed it in the rat, extrapolated it to mouse, monkey and dog, validated each against in-vivo plasma (and, in the mouse, tissue) data, and then used it to make a first-in-human prediction. They made a second, independent human prediction by interspecies allometric scaling.

This paper therefore contributes six models to nlmixr2lib.

The six models contributed by Chen 2016.
Model Role
Chen_2016_deoxypodophyllotoxin_mouse_pbpk PBPK, mouse (also the only tissue-level validation)
Chen_2016_deoxypodophyllotoxin_rat_pbpk PBPK, rat – the reference species the others are derived from
Chen_2016_deoxypodophyllotoxin_monkey_pbpk PBPK, cynomolgus monkey
Chen_2016_deoxypodophyllotoxin_dog_pbpk PBPK, beagle dog
Chen_2016_deoxypodophyllotoxin_human_pbpk PBPK, human – forward prediction, 16 mg IV
Chen_2016_deoxypodophyllotoxin_human_allometric Two-compartment human model from Dedrick-plot allometric scaling

The five PBPK files share one structure and differ only in the physiological and metabolic constants of their species. They are kept as separate files rather than one covariate-switched file because the authors ran them as five separate parameterisations, each drawing its physiology from a different literature source (Table 1 footnotes a-g).

Population

The in-vivo data are four single-dose intravenous animal studies.

Study populations (Methods, ‘In vivo Pharmacokinetic Study’).
Model Species N Weight Doses
mouse_pbpk mouse (ICR, male and female) 132 about 20 g Intravenous bolus 12.5 and 25.0 mg/kg (tail vein), two groups.
rat_pbpk rat NA 0.25 kg (the body weight underlying the Table 1 physiology) Intravenous bolus 1.0, 2.0 and 4.0 mg/kg.
monkey_pbpk monkey (cynomolgus, male and female) 12 3.51 +/- 0.85 kg Intravenous bolus 0.5, 1.0 and 2.0 mg/kg via the cephalic vein.
dog_pbpk dog (beagle, male and female) 6 8.65 +/- 0.64 kg Intravenous bolus 0.3 mg/kg via the cephalic vein.
human_pbpk human 0 70 kg (the reference body weight of the Table 1 human physiology) Simulated single intravenous bolus of 16 mg (not mg/kg) in a 70 kg adult.

The mouse, monkey and dog studies were run for this paper. The rat plasma profiles were taken from the authors’ earlier publication (Liu et al. 2016), so Chen 2016 does not report the rat cohort size, age or sampling schedule. There are no human data at all – the human dose of 16 mg was chosen from the monkey maximum tolerated dose of 4 mg/kg, and both human models are forward predictions.

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

Source trace

Every value in the six model files comes from one of the locations below.

Source location for every model equation and parameter.
Quantity Source location
Organ volumes, blood flows, Kt:pl for all 5 species Table 1 (per-species columns)
Brain permeability-surface product PS Table 1, ‘PS (mL/min)’ row
Brain vascular / extravascular volume split (3% / 97%) Table 1 footnote i
Liver blood flow = hepatic artery + GI tract + spleen Table 1 footnote h
Rest-of-body volume and flow by subtraction; Kt:pl,rest = 0.01 Table 1 footnote j
Vein / artery = 2/3 and 1/3 of blood volume Table 1 footnote k
PBSF, Km and Vmax for M2 and M7, Hill coefficient gamma Table 2 (per-species columns)
Microsomal protein yield (44.8 / 48.8 / 77.9 mg/g liver) Table 2 footnote
Fraction unbound in plasma, all 5 species Results, ‘Plasma Protein Binding of DPT in Five Species’
Perfusion-limited tissue ODE Methods, ‘PBPK Model Development’, eq. for non-elimination tissues
Arterial, venous and lung blood ODEs Methods, ‘PBPK Model Development’
Liver ODE including the PBSF x (CL_M2 + CL_M7) sink Methods, ‘PBPK Model Development’
Permeability-limited brain ODEs and the PS allometric scaling Results, ‘PBPK Model Development and Validation’
CL_M2 Michaelis-Menten sum and CL_M7 Hill equation Results, equations following the brain model
Rbp = 1 Methods, ‘PBPK Model Development’ (citing Poulin and Theil 2002)
Kt:pl = Kt:pl,rat x fu / fu,rat Results, ‘PBPK Model Development and Validation’
Human biexponential C = 0.628 e^-0.147t + 0.066 e^-0.011t Results, ‘Interspecies Allometric Scaling’
Allometric regressions CL = 44.28 W^0.832, Vss = 1.99 W^0.857 Results, ‘Interspecies Allometric Scaling’
Predicted and observed NCA parameters used as targets below Table 3

Two quantities are not printed in the paper and are flagged here and in the model files:

  • Molecular weight, 398.4 g/mol. A chemical constant for the DPT formula C22H22O7 (PubChem CID 73435), not a fitted parameter. It is needed purely as a unit bridge: the ODEs run on a mass concentration scale while the Table 2 Michaelis constants are in uM.
  • Residual error and inter-individual variances. The paper’s visual predictive checks fitted a proportional residual error and exponential inter-individual variability on hepatic blood flow and metabolic velocity in Phoenix NLME, but reports none of the estimates. They are encoded as fixed(0) rather than invented.

Model structure 

#> Preclinical (mouse). PBPK (whole-body, 13 tissue compartments, coded in Phoenix WinNonlin). Deoxypodophyllotoxin (DPT) disposition in the mouse after intravenous bolus dosing (Chen et al. 2016, Front Pharmacol). Adipose, liver, muscle, lung, kidney, brain, heart, spleen, skin, gastrointestinal tract and rest of body, plus arterial and venous blood pools. All tissues are perfusion-rate limited except the brain, which is permeability-limited and split into vascular (3% of brain volume) and extravascular (97%) sub-compartments coupled by a permeability-surface product PS. The gastrointestinal tract and spleen drain into the liver rather than into venous blood. Elimination is hepatic only, as the sum of two in-vitro metabolic pathways: Michaelis-Menten formation of M2 and auto-activating Hill-type formation of M7, both driven by the unbound hepatic outflow concentration and scaled to the whole liver by the microsomal protein amount PBSF. The blood-to-plasma ratio is assumed to be 1. This is the only species for which the paper reports tissue concentrations, and the brain permeability-surface product PS was estimated on these mouse brain data before being allometrically scaled to the other four species. Deterministic: the publication reports no inter-individual variance estimates and no residual-error magnitude, so the model is intended for typical-value simulation.

Eleven perfusion-rate-limited organs plus arterial and venous blood. The gastrointestinal tract and spleen drain into the liver, not into venous blood, so the liver receives hepatic-artery, portal and splenic inflow. The brain is the one exception to perfusion-rate limitation: a preliminary perfusion-limited brain over-predicted mouse brain concentrations, so the authors split it into a vascular and an extravascular sub-compartment coupled by a permeability-surface product.

mod_mouse <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_mouse_pbpk"))
mod_rat <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_rat_pbpk"))
mod_monkey <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_monkey_pbpk"))
mod_dog <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_dog_pbpk"))
mod_human <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_human_pbpk"))
mod_allom <- rxode2::rxode(readModelDb("Chen_2016_deoxypodophyllotoxin_human_allometric"))

mod_mouse$state
#>  [1] "venous"              "lung"                "arterial"           
#>  [4] "adipose"             "muscle"              "kidney"             
#>  [7] "heart"               "spleen"              "skin"               
#> [10] "gut"                 "other"               "brain_vascular"     
#> [13] "brain_extravascular" "liver"

Structural checks on the transcribed constants

Before simulating anything, four internal consistency relations that the paper states in prose or footnotes are checked against the transcribed ini() values. Each would fail loudly if a Table 1 or Table 2 column had been read off by one.

# Named vector of the fixed-effect values of a packaged model.
dpt_theta <- function(model_name) {
  rxode2::rxode(readModelDb(model_name))$theta
}

dpt_kp <- function(model_name) {
  th <- dpt_theta(model_name)
  lkp <- th[grepl("^lkp_", names(th))]
  stats::setNames(exp(lkp), sub("^lkp_", "", names(lkp)))
}

th <- lapply(
  stats::setNames(dpt_models[1:5], c("mouse", "rat", "monkey", "dog", "human")),
  dpt_theta
)
kp <- lapply(
  stats::setNames(dpt_models[1:5], c("mouse", "rat", "monkey", "dog", "human")),
  dpt_kp
)

1. Kt:pl is the rat value rescaled by the unbound fraction

The paper’s central assumption is that the unbound tissue-to-plasma ratio is identical across species, so Kt:pl = Kt:pl,rat * fu / fu,rat for every tissue of every species (Results, “PBPK Model Development and Validation”).

Table 1 prints Kt:pl to two decimal places (four for the rest-of-body column), and the rat column it is all derived from carries that same rounding, so the correct criterion is agreement to within one unit in the last printed place – not a percentage. A transposed or mis-read column would be out by far more than one digit.

fu_all <- vapply(th, function(x) unname(x[["fu"]]), numeric(1))
kp_rat <- kp$rat

kp_check <- do.call(rbind, lapply(names(kp), function(sp) {
  expected <- kp_rat * fu_all[[sp]] / fu_all[["rat"]]
  printed <- kp[[sp]][names(kp_rat)]
  # The last printed place in Table 1: 1e-4 for the rest-of-body column,
  # 1e-2 everywhere else.
  ulp <- ifelse(printed < 0.01, 1e-4, 1e-2)
  data.frame(
    species = sp,
    tissue = names(kp_rat),
    printed = unname(printed),
    expected = unname(expected),
    dev_ulp = unname(abs(printed - expected) / ulp)
  )
}))

stopifnot(all(kp_check$dev_ulp < 1))
stopifnot(stats::median(kp_check$dev_ulp) < 0.2)

cat(sprintf(
  paste0(
    "All %d Kt:pl values across 5 species reproduce Kt:pl,rat * fu / fu_rat\n",
    "to within one unit in the last printed place (median %.2f, max %.2f).\n"
  ),
  nrow(kp_check), stats::median(kp_check$dev_ulp), max(kp_check$dev_ulp)
))
#> All 55 Kt:pl values across 5 species reproduce Kt:pl,rat * fu / fu_rat
#> to within one unit in the last printed place (median 0.15, max 0.72).

The handful of cells above half a unit are listed below. Most are the muscle row: reproducing the mouse, monkey, dog and human muscle Kt:pl exactly requires a rat muscle Kt:pl near 0.76 rather than the 0.75 Table 1 prints, which is simply the rounding of the rat reference column propagating outward.

kp_check[kp_check$dev_ulp > 0.5, ] |>
  dplyr::mutate(
    printed = round(printed, 5),
    expected = round(expected, 5),
    dev_ulp = round(dev_ulp, 2)
  ) |>
  knitr::kable(
    row.names = FALSE,
    caption = "Kt:pl cells more than half a printed digit from the fu-scaling rule."
  )
Kt:pl cells more than half a printed digit from the fu-scaling rule.
species tissue printed expected dev_ulp
mouse muscle 0.4000 0.39378 0.62
mouse brain 1.4600 1.45438 0.56
mouse spleen 0.5500 0.55655 0.65
mouse other 0.0052 0.00525 0.50
human liver 2.0500 2.04468 0.53
human muscle 0.9200 0.91280 0.72

2. PBSF is the microsomal protein yield times the liver weight

Table 2’s footnote defines PBSF as the product of the microsomal protein yield and the liver weight from Table 1. Yields are 44.8 mg/g liver for mouse and rat, 48.8 for monkey and human, and 77.9 for dog.

mppgl <- c(mouse = 44.8, rat = 44.8, monkey = 48.8, dog = 77.9, human = 48.8)
# Liver volume is in mL for the animals and L for human; assume density 1 g/mL.
liver_g <- vapply(
  names(th),
  function(sp) unname(th[[sp]][["v_liver"]]) * if (sp == "human") 1000 else 1,
  numeric(1)
)
pbsf_printed <- vapply(th, function(x) unname(x[["pbsf"]]), numeric(1))
pbsf_expected <- mppgl[names(pbsf_printed)] * liver_g[names(pbsf_printed)]

stopifnot(all(abs(pbsf_printed - pbsf_expected) < 0.01))
knitr::kable(
  data.frame(
    Species = names(pbsf_printed),
    `Yield (mg/g)` = unname(mppgl[names(pbsf_printed)]),
    `Liver (g)` = unname(liver_g[names(pbsf_printed)]),
    `PBSF printed` = unname(pbsf_printed),
    `PBSF = yield x liver` = unname(pbsf_expected),
    check.names = FALSE
  ),
  caption = "Table 2 PBSF reproduced exactly from its own footnote."
)
Table 2 PBSF reproduced exactly from its own footnote.
Species Yield (mg/g) Liver (g) PBSF printed PBSF = yield x liver
mouse 44.8 1.10 49.28 49.28
rat 44.8 9.15 409.92 409.92
monkey 48.8 108.00 5270.40 5270.40
dog 77.9 213.00 16592.70 16592.70
human 48.8 1690.00 82472.00 82472.00

3. PS scales allometrically from the fitted mouse value

PS_i = PS_mouse * (W_i / W_mouse)^0.67, with PS estimated only in the mouse.

bw_all <- c(mouse = 0.02, rat = 0.25, monkey = 4, dog = 8.5, human = 70)
ps_printed <- vapply(th, function(x) unname(x[["ps_brain"]]), numeric(1))
ps_expected <- unname(ps_printed[["mouse"]]) *
  (bw_all[names(ps_printed)] / 0.02)^0.67
ps_ratio <- ps_printed / ps_expected

# PS_mouse is printed to only two significant figures (0.0024), and that
# rounding propagates to every scaled species -- hence a relative rather than
# absolute tolerance. The ratios are all ~1.02, consistent with an unrounded
# PS_mouse of about 0.00244.
stopifnot(all(abs(ps_ratio - 1) < 0.03))
knitr::kable(
  data.frame(
    Species = names(ps_printed),
    `Body weight (kg)` = unname(bw_all[names(ps_printed)]),
    `PS printed (mL/min)` = unname(ps_printed),
    `PS scaled from mouse` = round(unname(ps_expected), 5),
    Ratio = round(unname(ps_ratio), 3),
    check.names = FALSE
  ),
  caption = "Table 1 PS row reproduced by the 0.67-power allometric rule."
)
Table 1 PS row reproduced by the 0.67-power allometric rule.
Species Body weight (kg) PS printed (mL/min) PS scaled from mouse Ratio
mouse 0.02 0.0024 0.00240 1.000
rat 0.25 0.0133 0.01304 1.020
monkey 4.00 0.0849 0.08354 1.016
dog 8.50 0.1407 0.13843 1.016
human 70.00 0.5780 0.56851 1.017

4. Venous return balances cardiac output

Every organ that drains to venous blood, plus the liver (which collects the gastrointestinal and splenic flows), must sum to the cardiac output.

venous_organs <- c(
  "adipose", "muscle", "kidney", "brain", "heart", "skin", "other", "liver"
)
flow_check <- do.call(rbind, lapply(names(th), function(sp) {
  x <- th[[sp]]
  inflow <- sum(vapply(venous_organs, function(o) unname(x[[paste0("q_", o)]]), numeric(1)))
  data.frame(
    species = sp,
    cardiac_output = unname(x[["q_lung"]]),
    venous_inflow = inflow,
    pct_diff = 100 * (inflow / unname(x[["q_lung"]]) - 1)
  )
}))

# Mouse, monkey, dog and human balance to the printed precision. The rat
# column of Table 1 is internally inconsistent by 0.83 mL/min (1.0%) -- see
# the Errata; it is transcribed as printed, not adjusted.
stopifnot(all(abs(flow_check$pct_diff) < 1.1))
knitr::kable(
  flow_check,
  row.names = FALSE,
  digits = c(0, 2, 2, 2),
  caption = "Venous inflow versus cardiac output, per species."
)
Venous inflow versus cardiac output, per species.
species cardiac_output venous_inflow pct_diff
mouse 8.00 8.00 0.00
rat 83.90 83.07 -0.99
monkey 893.78 893.77 0.00
dog 968.33 968.32 0.00
human 5.60 5.60 0.00

Simulation

All simulations are deterministic typical-value solves: the models carry no random effects, so a single subject per dose group is the whole cohort. The observation grid is dense in the first two minutes because an intravenous bolus into a small venous blood pool produces a short, sharp mixing spike whose area is a real part of the AUC.

# Observation grid: dense over the bolus mixing spike, coarser later.
dpt_grid <- function(tmax) {
  sort(unique(c(
    seq(0, 2, by = 0.005),
    seq(2, 30, by = 0.1),
    seq(30, tmax, by = 1)
  )))
}

# Solve one species at one dose. `cmt = "venous"` names the ODE state directly
# on both the dose and the observation rows.
dpt_solve <- function(model, dose_amount, tmax, treatment, id) {
  ev <- rxode2::et(amt = dose_amount, cmt = "venous") |>
    rxode2::et(dpt_grid(tmax), cmt = "venous")
  out <- as.data.frame(rxode2::rxSolve(model, ev, returnType = "data.frame"))
  out$treatment <- treatment
  out$id <- id
  out$amt_dosed <- dose_amount
  out
}
# Body weights are the Table 1 reference weights, since the Table 1 physiology
# is tabulated for those weights.
arms <- tibble::tribble(
  ~species, ~dose_mgkg, ~bw, ~tmax,
  "mouse", 12.5, 0.02, 240,
  "mouse", 25.0, 0.02, 240,
  "rat", 1.0, 0.25, 240,
  "rat", 2.0, 0.25, 240,
  "rat", 4.0, 0.25, 240,
  "monkey", 0.5, 4.00, 240,
  "monkey", 1.0, 4.00, 240,
  "monkey", 2.0, 4.00, 240
) |>
  dplyr::mutate(
    # Animal models use ug and mL, so a mg/kg dose becomes ug.
    dose_amount = dose_mgkg * bw * 1000,
    treatment = paste0(species, " ", dose_mgkg, " mg/kg"),
    id = dplyr::row_number()
  )

models <- list(
  mouse = mod_mouse, rat = mod_rat, monkey = mod_monkey,
  dog = mod_dog, human = mod_human
)

sim_animals <- dplyr::bind_rows(lapply(seq_len(nrow(arms)), function(i) {
  a <- arms[i, ]
  dpt_solve(models[[a$species]], a$dose_amount, a$tmax, a$treatment, a$id)
}))

# The dog is dosed once and sampled to 840 min, the paper's sampling window.
sim_dog <- dpt_solve(mod_dog, 0.3 * 8.5 * 1000, 840, "dog 0.3 mg/kg", 9L)

# The human PBPK model works in mg and L, so 16 mg is the dose as given.
sim_human <- dpt_solve(mod_human, 16, 240, "human 16 mg (PBPK)", 10L)

nrow(sim_animals) + nrow(sim_dog) + nrow(sim_human)
#> [1] 9510

Replicating Figure 4 

plot_dat <- dplyr::bind_rows(sim_animals, sim_dog) |>
  dplyr::mutate(
    species = factor(
      sub(" .*$", "", treatment),
      levels = c("rat", "mouse", "monkey", "dog")
    )
  ) |>
  dplyr::filter(time > 0, Cc > 0)

ggplot2::ggplot(plot_dat, ggplot2::aes(time, Cc, colour = treatment)) +
  ggplot2::geom_line() +
  ggplot2::facet_wrap(~species, scales = "free") +
  ggplot2::scale_y_log10() +
  ggplot2::labs(
    x = "Time (min)",
    y = "Plasma DPT (ug/mL)",
    colour = NULL
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "bottom")
Replicates Figures 4A-D of Chen 2016: predicted DPT plasma concentration-time profiles in rat, mouse, monkey and dog after intravenous dosing.

Replicates Figures 4A-D of Chen 2016: predicted DPT plasma concentration-time profiles in rat, mouse, monkey and dog after intravenous dosing. 

sim_allom <- {
  ev <- rxode2::et(amt = 16, cmt = "central") |>
    rxode2::et(dpt_grid(240), cmt = "central")
  out <- as.data.frame(rxode2::rxSolve(mod_allom, ev, returnType = "data.frame"))
  out$treatment <- "human 16 mg (allometric)"
  out$id <- 11L
  out$amt_dosed <- 16
  out
}

dplyr::bind_rows(sim_human, sim_allom) |>
  dplyr::filter(time > 0, Cc > 0) |>
  ggplot2::ggplot(ggplot2::aes(time, Cc, colour = treatment)) +
  ggplot2::geom_line() +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time (min)", y = "Plasma DPT (ug/mL)", colour = NULL) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "bottom")
Replicates Figure 4E of Chen 2016: the two independent human predictions, whole-body PBPK and Dedrick-plot allometric scaling, after a 16 mg intravenous bolus.

Replicates Figure 4E of Chen 2016: the two independent human predictions, whole-body PBPK and Dedrick-plot allometric scaling, after a 16 mg intravenous bolus.

The allometric model is an exact algebraic re-parameterisation of the biexponential the paper prints, so it can be checked against that equation to machine precision rather than merely by eye.

published_biexp <- 0.628 * exp(-0.147 * sim_allom$time) +
  0.066 * exp(-0.011 * sim_allom$time)
max_rel_diff <- max(abs(sim_allom$Cc / published_biexp - 1))

stopifnot(max_rel_diff < 1e-3)
cat(sprintf(
  "Allometric model versus the printed C = 0.628 e^-0.147t + 0.066 e^-0.011t:\n  maximum relative difference %.2e over 0-240 min.\n",
  max_rel_diff
))
#> Allometric model versus the printed C = 0.628 e^-0.147t + 0.066 e^-0.011t:
#>   maximum relative difference 8.95e-05 over 0-240 min.

Mass balance

The intended sink in the PBPK models is hepatic metabolism, which the model exposes as rate_metabolism. Total body amount at any time should therefore equal the dose minus the integral of that rate. This exercises the dose placement, every inter-organ flow term and the unit bridge in the metabolic arm at once.

There is a second sink, and it is not intentional. Because Table 1’s rat column does not balance (check 4 above), the rat venous pool receives 83.07 mL/min of inflow while discharging 83.90 mL/min to the lung, so the transcribed rat physiology destroys drug at a rate of (qc - venous inflow) * C_venous. That is a property of the published table, not of the transcription, and the values are encoded as printed rather than adjusted to close the balance. Accounting for it explicitly turns the defect into a quantitative check.

dpt_states <- mod_mouse$state

dpt_cumtrap <- function(tt, y) {
  c(0, cumsum(diff(tt) * (utils::head(y, -1) + utils::tail(y, -1)) / 2))
}

dpt_mass_balance <- function(sim_one, leak_flow) {
  total <- rowSums(sim_one[, dpt_states, drop = FALSE])
  tt <- sim_one$time
  dose <- sim_one$amt_dosed[1]
  cum_met <- dpt_cumtrap(tt, sim_one$rate_metabolism)
  cum_leak <- dpt_cumtrap(tt, leak_flow * sim_one$Cc)
  c(
    metabolism_only = max(abs((total + cum_met) / dose - 1)),
    with_flow_imbalance = max(abs((total + cum_met + cum_leak) / dose - 1))
  )
}

# Positive = cardiac output exceeds the summed venous inflow, i.e. a leak.
leak <- stats::setNames(
  flow_check$cardiac_output - flow_check$venous_inflow,
  flow_check$species
)

mb <- rbind(
  mouse = dpt_mass_balance(dplyr::filter(sim_animals, id == 1), leak[["mouse"]]),
  rat = dpt_mass_balance(dplyr::filter(sim_animals, id == 3), leak[["rat"]]),
  monkey = dpt_mass_balance(dplyr::filter(sim_animals, id == 6), leak[["monkey"]]),
  dog = dpt_mass_balance(sim_dog, leak[["dog"]]),
  human = dpt_mass_balance(sim_human, leak[["human"]])
)

# Once the published rat flow imbalance is accounted for, every species closes.
stopifnot(all(mb[, "with_flow_imbalance"] < 0.01))

# The four balanced species close on metabolism alone; the rat does not, and
# the size of its gap is the signature of the Table 1 inconsistency.
stopifnot(all(mb[rownames(mb) != "rat", "metabolism_only"] < 0.01))
stopifnot(mb[["rat", "metabolism_only"]] > 0.01)

knitr::kable(
  data.frame(
    Species = rownames(mb),
    `Leak flow (vol/min)` = round(unname(leak[rownames(mb)]), 3),
    `Max error, metabolism only` = signif(unname(mb[, "metabolism_only"]), 3),
    `Max error, incl. flow imbalance` = signif(unname(mb[, "with_flow_imbalance"]), 3),
    check.names = FALSE
  ),
  row.names = FALSE,
  caption = "Mass balance. Only the rat needs the flow-imbalance term, and it needs exactly that term."
)
Mass balance. Only the rat needs the flow-imbalance term, and it needs exactly that term.
Species Leak flow (vol/min) Max error, metabolism only Max error, incl. flow imbalance
mouse 0.00 6.16e-05 6.16e-05
rat 0.83 5.26e-02 3.28e-03
monkey 0.01 8.39e-05 3.30e-05
dog 0.01 1.15e-03 3.00e-06
human 0.00 1.98e-05 1.98e-05

PKNCA validation

Non-compartmental analysis of the simulated profiles, done with PKNCA. The animal models are in ug / mL / min, the human PBPK model in mg / L / min, so each unit system gets its own PKNCA run.

dpt_nca <- function(sim, concu, timeu, doseu, tmax) {
  conc_df <- sim |>
    dplyr::filter(!is.na(Cc)) |>
    dplyr::select(id, time, Cc, treatment)

  # Defensive time-zero row (the grid already contains time 0; the existing
  # row wins through .keep_all).
  conc_df <- dplyr::bind_rows(
    conc_df,
    conc_df |>
      dplyr::distinct(id, treatment) |>
      dplyr::mutate(time = 0, Cc = 0)
  ) |>
    dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
    dplyr::arrange(id, treatment, time)

  dose_df <- sim |>
    dplyr::distinct(id, treatment, amt_dosed) |>
    dplyr::mutate(time = 0) |>
    dplyr::rename(amt = amt_dosed)

  conc_obj <- PKNCA::PKNCAconc(
    conc_df, Cc ~ time | treatment + id,
    concu = concu, timeu = timeu
  )
  dose_obj <- PKNCA::PKNCAdose(
    dose_df, amt ~ time | treatment + id,
    doseu = doseu, route = "intravascular"
  )

  intervals <- data.frame(
    start = 0,
    end = tmax,
    auclast = TRUE,
    aucinf.obs = TRUE,
    half.life = TRUE,
    cl.obs = TRUE,
    vss.obs = TRUE,
    mrt.obs = TRUE
  )

  PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
}
nca_animals <- dpt_nca(sim_animals, "ug/mL", "min", "ug", 240)
nca_dog <- dpt_nca(sim_dog, "ug/mL", "min", "ug", 840)
nca_human <- dpt_nca(sim_human, "mg/L", "min", "mg", 240)

nca_raw <- dplyr::bind_rows(
  as.data.frame(nca_animals$result),
  as.data.frame(nca_dog$result),
  as.data.frame(nca_human$result)
)

PKNCA returns clearance and steady-state volume in the absolute units of each system. The paper normalises them per kilogram for the animals (mL/min/kg and L/kg) and reports them absolute for the human (L/min and L), so the simulated values are put on the same footing here.

arm_bw <- dplyr::bind_rows(
  arms |> dplyr::select(treatment, bw),
  tibble::tibble(treatment = "dog 0.3 mg/kg", bw = 8.5),
  tibble::tibble(treatment = "human 16 mg (PBPK)", bw = 1)
)

nca_sim <- nca_raw |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  dplyr::left_join(arm_bw, by = "treatment") |>
  dplyr::mutate(
    PPORRES = dplyr::case_when(
      # mL/min -> mL/min/kg for animals; the human row has bw = 1 (L/min).
      PPTESTCD == "cl.obs" ~ PPORRES / bw,
      # mL -> L/kg for animals; the human row is already L.
      PPTESTCD == "vss.obs" & treatment != "human 16 mg (PBPK)" ~
        PPORRES / bw / 1000,
      TRUE ~ PPORRES
    )
  ) |>
  dplyr::select(treatment, PPTESTCD, PPORRES)

Comparison against the published predicted values (Table 3)

The primary target is the Predicted row of Table 3: those are the numbers the paper’s own implementation of this model produced, so reproducing them tests the transcription rather than the science.

published_predicted <- tibble::tribble(
  ~treatment, ~auclast, ~aucinf.obs, ~cl.obs, ~vss.obs, ~mrt.obs, ~half.life,
  "mouse 12.5 mg/kg", 213.85, 214.69, 58.22, 1.66, 28.58, 33.47,
  "mouse 25 mg/kg", 486.31, 488.18, 51.21, 1.44, 28.05, 33.45,
  "rat 1 mg/kg", 12.62, 13.63, 73.38, 3.68, 50.10, 69.75,
  "rat 2 mg/kg", 26.16, 27.30, 73.27, 3.70, 50.51, 70.59,
  "rat 4 mg/kg", 52.34, 54.62, 73.23, 3.70, 50.50, 70.59,
  "monkey 0.5 mg/kg", 21.05, 21.88, 22.85, 0.46, 20.28, 44.62,
  "monkey 1 mg/kg", 45.60, 46.29, 21.60, 0.44, 20.43, 44.79,
  "monkey 2 mg/kg", 102.19, 102.81, 19.45, 0.40, 20.72, 44.91,
  "dog 0.3 mg/kg", 295.24, 298.30, 1.01, 0.18, 175.52, 128.94,
  "human 16 mg (PBPK)", 11.07, 11.50, 1.39, 61.94, 44.52, 70.22
)

cmp_predicted <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_sim,
  reference = published_predicted,
  by = "treatment",
  tolerance_pct = 20
)

knitr::kable(
  cmp_predicted,
  caption = paste(
    "Simulated versus Chen 2016 Table 3 PREDICTED values.",
    "* marks a difference of more than 20%.",
    "Units: AUC ug*min/mL (mg*min/L for human, numerically identical);",
    "CL mL/min/kg (L/min for human); Vss L/kg (L for human); MRT and t1/2 min."
  ),
  align = c("l", "l", "r", "r", "r")
)
Simulated versus Chen 2016 Table 3 PREDICTED values. * marks a difference of more than 20%. Units: AUC ugmin/mL (mgmin/L for human, numerically identical); CL mL/min/kg (L/min for human); Vss L/kg (L for human); MRT and t1/2 min.
NCA parameter treatment Reference Simulated % diff
AUC0-∞ (obs) mouse 12.5 mg/kg 215 251 +17.1%
AUC0-∞ (obs) mouse 25 mg/kg 488 553 +13.3%
AUC0-∞ (obs) rat 1 mg/kg 13.6 16.4 +20.1%*
AUC0-∞ (obs) rat 2 mg/kg 27.3 32.7 +19.9%
AUC0-∞ (obs) rat 4 mg/kg 54.6 65.5 +19.9%
AUC0-∞ (obs) monkey 0.5 mg/kg 21.9 23.5 +7.4%
AUC0-∞ (obs) monkey 1 mg/kg 46.3 49.5 +6.8%
AUC0-∞ (obs) monkey 2 mg/kg 103 109 +5.9%
AUC0-∞ (obs) dog 0.3 mg/kg 298 298 -0.2%
AUC0-∞ (obs) human 16 mg (PBPK) 11.5 11.4 -1.0%
AUClast mouse 12.5 mg/kg 214 250 +17.1%
AUClast mouse 25 mg/kg 486 551 +13.4%
AUClast rat 1 mg/kg 12.6 15.9 +25.8%*
AUClast rat 2 mg/kg 26.2 31.8 +21.4%*
AUClast rat 4 mg/kg 52.3 63.5 +21.4%*
AUClast monkey 0.5 mg/kg 21 23.4 +11.0%
AUClast monkey 1 mg/kg 45.6 49.2 +7.8%
AUClast monkey 2 mg/kg 102 108 +6.0%
AUClast dog 0.3 mg/kg 295 295 -0.1%
AUClast human 16 mg (PBPK) 11.1 10.9 -1.6%
t½ mouse 12.5 mg/kg 33.5 33.5 -0.0%
t½ mouse 25 mg/kg 33.4 33.4 -0.0%
t½ rat 1 mg/kg 69.8 69.3 -0.6%
t½ rat 2 mg/kg 70.6 69.3 -1.8%
t½ rat 4 mg/kg 70.6 69.3 -1.8%
t½ monkey 0.5 mg/kg 44.6 44.9 +0.7%
t½ monkey 1 mg/kg 44.8 44.9 +0.3%
t½ monkey 2 mg/kg 44.9 44.9 -0.0%
t½ dog 0.3 mg/kg 129 127 -1.2%
t½ human 16 mg (PBPK) 70.2 79.1 +12.6%
CL/F mouse 12.5 mg/kg 58.2 49.7 -14.6%
CL/F mouse 25 mg/kg 51.2 45.2 -11.8%
CL/F rat 1 mg/kg 73.4 61.1 -16.8%
CL/F rat 2 mg/kg 73.3 61.1 -16.6%
CL/F rat 4 mg/kg 73.2 61.1 -16.6%
CL/F monkey 0.5 mg/kg 22.8 21.3 -6.9%
CL/F monkey 1 mg/kg 21.6 20.2 -6.4%
CL/F monkey 2 mg/kg 19.4 18.4 -5.6%
CL/F dog 0.3 mg/kg 1.01 1.01 -0.3%
CL/F human 16 mg (PBPK) 1.39 1.41 +1.1%
Vss/F mouse 12.5 mg/kg 1.66 1.22 -26.6%*
Vss/F mouse 25 mg/kg 1.44 1.12 -22.1%*
Vss/F rat 1 mg/kg 3.68 2.31 -37.2%*
Vss/F rat 2 mg/kg 3.7 2.31 -37.5%*
Vss/F rat 4 mg/kg 3.7 2.31 -37.6%*
Vss/F monkey 0.5 mg/kg 0.46 0.404 -12.1%
Vss/F monkey 1 mg/kg 0.44 0.388 -11.8%
Vss/F monkey 2 mg/kg 0.4 0.359 -10.2%
Vss/F dog 0.3 mg/kg 0.18 0.174 -3.2%
Vss/F human 16 mg (PBPK) 61.9 66.1 +6.8%
MRT mouse 12.5 mg/kg 28.6 24.5 -14.3%
MRT mouse 25 mg/kg 28 24.8 -11.5%
MRT rat 1 mg/kg 50.1 37.8 -24.5%*
MRT rat 2 mg/kg 50.5 37.8 -25.1%*
MRT rat 4 mg/kg 50.5 37.8 -25.1%*
MRT monkey 0.5 mg/kg 20.3 19 -6.3%
MRT monkey 1 mg/kg 20.4 19.2 -6.1%
MRT monkey 2 mg/kg 20.7 19.6 -5.6%
MRT dog 0.3 mg/kg 176 173 -1.4%
MRT human 16 mg (PBPK) 44.5 47.1 +5.7%

The paper’s own criterion for a successful prediction is a fold-error below two (Methods, “PBPK Model Development”). That is the hard gate here, applied to every parameter of every arm, plus two tighter gates that the reproduction comfortably meets.

ncaComparisonTable() formats its values for display, so the gates below are computed from the underlying numeric simulated and reference values rather than from the rendered table.

# Long, numeric form of the reference table for gating.
dpt_gate_frame <- function(simulated, reference) {
  reference |>
    tidyr::pivot_longer(
      -treatment,
      names_to = "PPTESTCD",
      values_to = "Reference"
    ) |>
    dplyr::inner_join(
      dplyr::rename(simulated, Simulated = PPORRES),
      by = c("treatment", "PPTESTCD")
    ) |>
    dplyr::mutate(
      fold_error = pmax(Simulated / Reference, Reference / Simulated)
    )
}

gate <- dpt_gate_frame(nca_sim, published_predicted)
stopifnot(nrow(gate) == 6 * 10)

# Gate 1 (the paper's own criterion): every parameter of every arm within
# 2-fold of the value the paper's implementation produced.
stopifnot(all(gate$fold_error < 2))

# Gate 2: the terminal half-life, the parameter that depends on the disposition
# structure rather than on the integration of the bolus spike.
#
# Half-life is only well determined when the slowest tissue has had time to
# equilibrate within the sampling window. The slowest tissue is always adipose,
# whose equilibration time constant is v_adipose * kp_adipose / q_adipose. That
# is a small fraction of the sampling window in all four animals, but roughly
# four times the window in the human, whose adipose Kt:pl of 26.33 is the
# largest value anywhere in Table 1. So the human arm is the one place where a
# 240-minute NCA cannot see the true terminal phase and the apparent half-life
# depends on which points the regression picks.
tau_adipose <- vapply(names(th), function(sp) {
  unname(th[[sp]][["v_adipose"]]) * unname(kp[[sp]][["adipose"]]) /
    unname(th[[sp]][["q_adipose"]])
}, numeric(1))
window <- c(mouse = 240, rat = 240, monkey = 240, dog = 840, human = 240)
resolved <- tau_adipose < window

stopifnot(identical(unname(which(!resolved)), which(names(resolved) == "human")))

t_half <- gate[gate$PPTESTCD == "half.life", ]
stopifnot(nrow(t_half) == 10)

# Every arm within 15%...
stopifnot(all(abs(t_half$Simulated / t_half$Reference - 1) < 0.15))
# ...and every arm whose slowest tissue does equilibrate inside the window
# within 5%.
t_half_resolved <- t_half[!grepl("^human", t_half$treatment), ]
stopifnot(nrow(t_half_resolved) == 9)
stopifnot(all(abs(t_half_resolved$Simulated / t_half_resolved$Reference - 1) < 0.05))

knitr::kable(
  data.frame(
    Species = names(tau_adipose),
    `Adipose equilibration time constant (min)` = round(unname(tau_adipose), 0),
    `Sampling window (min)` = unname(window[names(tau_adipose)]),
    `Terminal phase resolved` = unname(resolved),
    check.names = FALSE
  ),
  row.names = FALSE,
  caption = "Why the human half-life is the one that is not pinned down by a 240-minute NCA."
)
Why the human half-life is the one that is not pinned down by a 240-minute NCA.
Species Adipose equilibration time constant (min) Sampling window (min) Terminal phase resolved
mouse 27 240 TRUE
rat 71 240 TRUE
monkey 52 240 TRUE
dog 19 840 TRUE
human 1013 240 FALSE 

# Gate 3: the two species whose profiles are least dominated by the bolus
# mixing spike reproduce AUC0-inf within 5%.
tight <- gate[
  gate$PPTESTCD == "aucinf.obs" &
    gate$treatment %in% c("dog 0.3 mg/kg", "human 16 mg (PBPK)"),
]
stopifnot(nrow(tight) == 2)
stopifnot(all(abs(tight$Simulated / tight$Reference - 1) < 0.05))

cat(sprintf(
  "Maximum fold-error across all %d parameter x arm comparisons: %.3f (%s, %s)\n",
  nrow(gate), max(gate$fold_error),
  gate$treatment[which.max(gate$fold_error)],
  gate$PPTESTCD[which.max(gate$fold_error)]
))
#> Maximum fold-error across all 60 parameter x arm comparisons: 1.601 (rat 4 mg/kg, vss.obs)

Comparison against the observed values (Table 3)

The observed column is the independent check – it is what the model was validated against rather than what it was built to reproduce. The paper’s success criterion is again a fold-error below two, which it reports as met for every animal parameter except monkey Vss and MRT.

published_observed <- tibble::tribble(
  ~treatment, ~auclast, ~aucinf.obs, ~cl.obs, ~vss.obs, ~mrt.obs, ~half.life,
  "mouse 12.5 mg/kg", 175.14, 176.65, 70.76, 1.71, 24.19, 48.98,
  "mouse 25 mg/kg", 359.39, 369.10, 67.73, 2.54, 37.52, 58.82,
  "rat 1 mg/kg", 9.55, 10.59, 96.98, 6.07, 63.27, 80.20,
  "rat 2 mg/kg", 27.20, 29.23, 69.28, 4.29, 61.42, 84.42,
  "rat 4 mg/kg", 62.45, 67.66, 60.17, 3.80, 62.40, 97.66,
  "monkey 0.5 mg/kg", 12.80, 13.94, 37.98, 1.19, 32.90, 32.71,
  "monkey 1 mg/kg", 29.98, 31.15, 33.84, 1.28, 39.58, 45.73,
  "monkey 2 mg/kg", 72.76, 75.23, 27.25, 1.22, 46.26, 57.52,
  "dog 0.3 mg/kg", 228.96, 240.84, 1.27, 0.35, 276.59, 186.77
)

cmp_observed <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_sim |> dplyr::filter(treatment != "human 16 mg (PBPK)"),
  reference = published_observed,
  by = "treatment",
  tolerance_pct = 20
)

knitr::kable(
  cmp_observed,
  caption = paste(
    "Simulated versus Chen 2016 Table 3 OBSERVED values (mean of the in-vivo",
    "cohort). * marks a difference of more than 20%. No human row: no human",
    "data exist."
  ),
  align = c("l", "l", "r", "r", "r")
)
Simulated versus Chen 2016 Table 3 OBSERVED values (mean of the in-vivo cohort). * marks a difference of more than 20%. No human row: no human data exist.
NCA parameter treatment Reference Simulated % diff
AUC0-∞ (obs) mouse 12.5 mg/kg 177 251 +42.3%*
AUC0-∞ (obs) mouse 25 mg/kg 369 553 +49.9%*
AUC0-∞ (obs) rat 1 mg/kg 10.6 16.4 +54.6%*
AUC0-∞ (obs) rat 2 mg/kg 29.2 32.7 +12.0%
AUC0-∞ (obs) rat 4 mg/kg 67.7 65.5 -3.2%
AUC0-∞ (obs) monkey 0.5 mg/kg 13.9 23.5 +68.6%*
AUC0-∞ (obs) monkey 1 mg/kg 31.2 49.5 +58.7%*
AUC0-∞ (obs) monkey 2 mg/kg 75.2 109 +44.8%*
AUC0-∞ (obs) dog 0.3 mg/kg 241 298 +23.7%*
AUClast mouse 12.5 mg/kg 175 250 +43.0%*
AUClast mouse 25 mg/kg 359 551 +53.4%*
AUClast rat 1 mg/kg 9.55 15.9 +66.2%*
AUClast rat 2 mg/kg 27.2 31.8 +16.7%
AUClast rat 4 mg/kg 62.4 63.5 +1.7%
AUClast monkey 0.5 mg/kg 12.8 23.4 +82.6%*
AUClast monkey 1 mg/kg 30 49.2 +64.0%*
AUClast monkey 2 mg/kg 72.8 108 +48.8%*
AUClast dog 0.3 mg/kg 229 295 +28.8%*
t½ mouse 12.5 mg/kg 49 33.5 -31.7%*
t½ mouse 25 mg/kg 58.8 33.4 -43.1%*
t½ rat 1 mg/kg 80.2 69.3 -13.6%
t½ rat 2 mg/kg 84.4 69.3 -17.9%
t½ rat 4 mg/kg 97.7 69.3 -29.0%*
t½ monkey 0.5 mg/kg 32.7 44.9 +37.4%*
t½ monkey 1 mg/kg 45.7 44.9 -1.7%
t½ monkey 2 mg/kg 57.5 44.9 -21.9%*
t½ dog 0.3 mg/kg 187 127 -31.8%*
CL/F mouse 12.5 mg/kg 70.8 49.7 -29.7%*
CL/F mouse 25 mg/kg 67.7 45.2 -33.3%*
CL/F rat 1 mg/kg 97 61.1 -37.0%*
CL/F rat 2 mg/kg 69.3 61.1 -11.8%
CL/F rat 4 mg/kg 60.2 61.1 +1.5%
CL/F monkey 0.5 mg/kg 38 21.3 -44.0%*
CL/F monkey 1 mg/kg 33.8 20.2 -40.2%*
CL/F monkey 2 mg/kg 27.2 18.4 -32.6%*
CL/F dog 0.3 mg/kg 1.27 1.01 -20.7%*
Vss/F mouse 12.5 mg/kg 1.71 1.22 -28.8%*
Vss/F mouse 25 mg/kg 2.54 1.12 -55.9%*
Vss/F rat 1 mg/kg 6.07 2.31 -61.9%*
Vss/F rat 2 mg/kg 4.29 2.31 -46.1%*
Vss/F rat 4 mg/kg 3.8 2.31 -39.2%*
Vss/F monkey 0.5 mg/kg 1.19 0.404 -66.0%*
Vss/F monkey 1 mg/kg 1.28 0.388 -69.7%*
Vss/F monkey 2 mg/kg 1.22 0.359 -70.5%*
Vss/F dog 0.3 mg/kg 0.35 0.174 -50.2%*
MRT mouse 12.5 mg/kg 24.2 24.5 +1.3%
MRT mouse 25 mg/kg 37.5 24.8 -33.9%*
MRT rat 1 mg/kg 63.3 37.8 -40.2%*
MRT rat 2 mg/kg 61.4 37.8 -38.4%*
MRT rat 4 mg/kg 62.4 37.8 -39.4%*
MRT monkey 0.5 mg/kg 32.9 19 -42.2%*
MRT monkey 1 mg/kg 39.6 19.2 -51.5%*
MRT monkey 2 mg/kg 46.3 19.6 -57.7%*
MRT dog 0.3 mg/kg 277 173 -37.4%* 
gate_obs <- dpt_gate_frame(
  dplyr::filter(nca_sim, treatment != "human 16 mg (PBPK)"),
  published_observed
)
stopifnot(nrow(gate_obs) == 6 * 9)

cat(sprintf(
  "Fold-error versus observed data: median %.2f, maximum %.2f (%s, %s).\n",
  stats::median(gate_obs$fold_error), max(gate_obs$fold_error),
  gate_obs$treatment[which.max(gate_obs$fold_error)],
  gate_obs$PPTESTCD[which.max(gate_obs$fold_error)]
))
#> Fold-error versus observed data: median 1.51, maximum 3.39 (monkey 2 mg/kg, vss.obs).

# The paper reports its own predicted monkey Vss and MRT as the only
# parameters outside 2-fold of the observations. The meaningful gate here is
# therefore about WHICH parameters miss, not how many: every EXPOSURE
# parameter -- AUC, clearance and half-life -- must stay inside 2-fold of the
# observed cohort mean, and only volume and residence-time parameters may fall
# outside.
exposure <- c("auclast", "aucinf.obs", "cl.obs", "half.life")
stopifnot(all(gate_obs$fold_error[gate_obs$PPTESTCD %in% exposure] < 2))

outside <- gate_obs[gate_obs$fold_error >= 2, ]
stopifnot(all(outside$PPTESTCD %in% c("vss.obs", "mrt.obs")))

knitr::kable(
  outside[order(-outside$fold_error), c("treatment", "PPTESTCD", "Reference", "Simulated", "fold_error")],
  row.names = FALSE,
  digits = 2,
  caption = paste(
    "The only parameters more than 2-fold from the observed cohort mean.",
    "All are steady-state volume or mean residence time; the paper reports",
    "the same class of miss for its own monkey predictions."
  )
)
The only parameters more than 2-fold from the observed cohort mean. All are steady-state volume or mean residence time; the paper reports the same class of miss for its own monkey predictions.
treatment PPTESTCD Reference Simulated fold_error
monkey 2 mg/kg vss.obs 1.22 0.36 3.39
monkey 1 mg/kg vss.obs 1.28 0.39 3.30
monkey 0.5 mg/kg vss.obs 1.19 0.40 2.94
rat 1 mg/kg vss.obs 6.07 2.31 2.63
monkey 2 mg/kg mrt.obs 46.26 19.57 2.36
mouse 25 mg/kg vss.obs 2.54 1.12 2.27
monkey 1 mg/kg mrt.obs 39.58 19.18 2.06
dog 0.3 mg/kg vss.obs 0.35 0.17 2.01

Both the paper’s own prediction and this reproduction under-predict steady-state volume and mean residence time while matching exposure. That is the expected consequence of deriving those two parameters by non-compartmental analysis of a 240-minute window when the adipose compartment – which holds most of the distribution volume – has not finished equilibrating.

Dose proportionality

Rat clearance is essentially hepatic-blood-flow limited (the unbound intrinsic clearance is two orders of magnitude above hepatic blood flow), so rat exposure should be strictly dose proportional across 1, 2 and 4 mg/kg. The paper’s own predicted rat clearances – 73.38, 73.27 and 73.23 mL/min/kg – say the same.

rat_cl <- nca_sim |>
  dplyr::filter(grepl("^rat", treatment), PPTESTCD == "cl.obs") |>
  dplyr::arrange(treatment)

stopifnot(nrow(rat_cl) == 3)
stopifnot(diff(range(rat_cl$PPORRES)) / mean(rat_cl$PPORRES) < 0.01)

# The mouse and monkey, whose Km values sit inside the simulated concentration
# range, must instead show clearance FALLING with dose -- the signature of
# saturable metabolism that the paper reports for the monkey.
sat_cl <- nca_sim |>
  dplyr::filter(PPTESTCD == "cl.obs", grepl("^(mouse|monkey)", treatment))
mouse_cl <- sat_cl$PPORRES[grepl("^mouse", sat_cl$treatment)]
monkey_cl <- sat_cl$PPORRES[grepl("^monkey", sat_cl$treatment)]
stopifnot(mouse_cl[2] < mouse_cl[1])
stopifnot(all(diff(monkey_cl) < 0))

cat(sprintf(
  "Rat CL across 1/2/4 mg/kg: %s mL/min/kg (dose independent).\n",
  paste(round(rat_cl$PPORRES, 2), collapse = ", ")
))
#> Rat CL across 1/2/4 mg/kg: 61.09, 61.08, 61.06 mL/min/kg (dose independent).

Mouse tissue distribution (Figure 5)

The mouse is the only species with tissue data. The seven tissues the paper assayed come out of the same solve, because they are ODE states of the model.

tissue_map <- c(
  c_heart = "Heart", c_liver = "Liver", c_lung = "Lung",
  c_muscle = "Muscle", c_brain = "Brain", c_kidney = "Kidney",
  c_spleen = "Spleen"
)

mouse_tissue <- sim_animals |>
  dplyr::filter(id == 1, time > 0) |>
  dplyr::select(time, dplyr::all_of(names(tissue_map))) |>
  tidyr::pivot_longer(-time, names_to = "state", values_to = "conc") |>
  dplyr::mutate(Tissue = unname(tissue_map[state]))

ggplot2::ggplot(mouse_tissue, ggplot2::aes(time, conc)) +
  ggplot2::geom_line() +
  ggplot2::geom_point(
    data = dplyr::filter(mouse_tissue, time %in% c(5, 15, 60, 120)),
    size = 1.6
  ) +
  ggplot2::facet_wrap(~Tissue, scales = "free_y") +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time (min)", y = "Tissue DPT (ug/mL)") +
  ggplot2::theme_bw()
Replicates Figures 5A-G of Chen 2016: predicted DPT concentrations in seven mouse tissues after a 12.5 mg/kg intravenous dose. Points mark the four times at which the paper assayed tissue.

Replicates Figures 5A-G of Chen 2016: predicted DPT concentrations in seven mouse tissues after a 12.5 mg/kg intravenous dose. Points mark the four times at which the paper assayed tissue.

The paper’s qualitative tissue findings are that DPT has a high affinity for adipose tissue and restricted brain penetration – the reason the brain was refined to a permeability-limited model in the first place. Both are asserted rather than left as prose.

mouse_sim <- dplyr::filter(sim_animals, id == 1)
auc120 <- function(x) {
  k <- mouse_sim$time <= 120
  tt <- mouse_sim$time[k]
  y <- x[k]
  sum(diff(tt) * (utils::head(y, -1) + utils::tail(y, -1)) / 2)
}

tissue_auc <- vapply(
  c("c_adipose", "c_liver", "c_muscle", "c_brain", "c_kidney", "c_heart",
    "c_spleen", "c_lung", "Cc"),
  function(v) auc120(mouse_sim[[v]]),
  numeric(1)
)
ratio_to_plasma <- tissue_auc / tissue_auc[["Cc"]]

# Adipose has by far the highest Kt:pl in Table 1 (11.36 in the mouse), so its
# exposure must exceed plasma and every other tissue's.
stopifnot(ratio_to_plasma[["c_adipose"]] > 1)
stopifnot(which.max(ratio_to_plasma) == which(names(ratio_to_plasma) == "c_adipose"))

# The permeability-limited brain must accumulate more slowly than a
# perfusion-limited tissue of similar Kt:pl: brain Kt:pl (1.46) is nearly
# three times heart Kt:pl (0.51), yet brain exposure over the first two hours
# is restricted rather than proportionally higher.
stopifnot(ratio_to_plasma[["c_brain"]] < 1.46)

knitr::kable(
  data.frame(
    Tissue = sub("^c_", "", names(ratio_to_plasma)),
    `AUC0-120 (ug*min/mL)` = round(unname(tissue_auc), 1),
    `Ratio to plasma` = round(unname(ratio_to_plasma), 2),
    check.names = FALSE
  ),
  caption = "Predicted mouse tissue exposure after 12.5 mg/kg (Figure 5H analogue)."
)
Predicted mouse tissue exposure after 12.5 mg/kg (Figure 5H analogue).
Tissue AUC0-120 (ug*min/mL) Ratio to plasma
adipose 2597.8 10.76
liver 104.1 0.43
muscle 96.3 0.40
brain 142.9 0.59
kidney 178.7 0.74
heart 123.1 0.51
spleen 132.8 0.55
lung 217.3 0.90
Cc 241.5 1.00

The two human predictions

The paper’s headline claim is that the mechanistic and the empirical route to a human prediction agree. Table 3 and the allometric section give both.

human_pbpk_nca <- nca_sim |>
  dplyr::filter(treatment == "human 16 mg (PBPK)") |>
  dplyr::select(PPTESTCD, PPORRES)

human_tbl <- tibble::tribble(
  ~Parameter, ~`PBPK (paper)`, ~`PBPK (simulated)`, ~`Allometric (paper)`,
  "AUC0-inf (ug*min/mL)", 11.50,
  round(human_pbpk_nca$PPORRES[human_pbpk_nca$PPTESTCD == "aucinf.obs"], 2), 10.25,
  "CL (L/min)", 1.39,
  round(human_pbpk_nca$PPORRES[human_pbpk_nca$PPTESTCD == "cl.obs"], 3), 1.56,
  "Vss (L)", 61.94,
  round(human_pbpk_nca$PPORRES[human_pbpk_nca$PPTESTCD == "vss.obs"], 2), 87.67,
  "t1/2 (min)", 70.22,
  round(human_pbpk_nca$PPORRES[human_pbpk_nca$PPTESTCD == "half.life"], 2), 63.08
)

knitr::kable(
  human_tbl,
  caption = "The paper's two independent human predictions, plus this package's reproduction of the PBPK one."
)
The paper’s two independent human predictions, plus this package’s reproduction of the PBPK one.
Parameter PBPK (paper) PBPK (simulated) Allometric (paper)
AUC0-inf (ug*min/mL) 11.50 11.390 10.25
CL (L/min) 1.39 1.405 1.56
Vss (L) 61.94 66.130 87.67
t1/2 (min) 70.22 79.100 63.08 
# The allometric model's own derived parameters, computed from its ini() rather
# than read off the paper, must land on the paper's printed allometric column.
allom_theta <- mod_allom$theta
allom_cl <- exp(unname(allom_theta[["lcl"]]))
allom_vss <- exp(unname(allom_theta[["lvc"]])) + exp(unname(allom_theta[["lvp"]]))

stopifnot(abs(allom_cl / 1.56 - 1) < 0.01)
stopifnot(abs(allom_vss / 87.67 - 1) < 0.01)
stopifnot(abs(16 / allom_cl / 10.25 - 1) < 0.01)

# The paper's claim is agreement, not identity: "The two approaches reached
# similar results." Clearance agrees within 15%, volume within 45%.
pbpk_cl <- human_pbpk_nca$PPORRES[human_pbpk_nca$PPTESTCD == "cl.obs"]
stopifnot(abs(allom_cl / pbpk_cl - 1) < 0.15)

cat(sprintf(
  "Human CL: PBPK %.2f L/min versus allometric %.2f L/min (%.0f%% apart).\n",
  pbpk_cl, allom_cl, 100 * abs(allom_cl / pbpk_cl - 1)
))
#> Human CL: PBPK 1.41 L/min versus allometric 1.56 L/min (11% apart).

The paper also prints the two allometric regressions the Dedrick plot was built from. They are reproduced here for completeness; they give a slightly different human clearance and volume than the reverse-transformed profile does, and the paper reports both.

allom_cl_70 <- 44.28 * 70^0.832 / 1000
allom_vss_70 <- 1.99 * 70^0.857

stopifnot(abs(allom_cl_70 / 1.52 - 1) < 0.02)
stopifnot(abs(allom_vss_70 / 75.65 - 1) < 0.02)

cat(sprintf(
  "CL = 44.28 * W^0.832 at 70 kg: %.2f L/min (paper 1.52).\nVss = 1.99 * W^0.857 at 70 kg: %.2f L (paper 75.65).\n",
  allom_cl_70, allom_vss_70
))
#> CL = 44.28 * W^0.832 at 70 kg: 1.52 L/min (paper 1.52).
#> Vss = 1.99 * W^0.857 at 70 kg: 75.88 L (paper 75.65).

Assumptions and deviations

Errata and transcription notes

  • Table 1’s rat column does not balance. The blood flows of the organs that drain to venous blood sum to 83.07 mL/min against a printed cardiac output (lung blood flow) of 83.90 mL/min, a 1.0% shortfall. The mouse, monkey, dog and human columns balance to the printed precision. The rat values are encoded exactly as printed; they are not adjusted to close the balance.
  • PS is printed to two significant figures in the mouse. PS_mouse = 0.0024 mL/min is the fitted value from which the other four species are scaled by (W/W_mouse)^0.67. Rescaling the printed 0.0024 reproduces the other four printed PS values only to about 2%, consistent with an unrounded PS_mouse near 0.00244. Each species’ own printed PS is used rather than a value rescaled from the mouse.
  • Vmax for M7 is printed in pmol/min/mg protein while Vmax for M2 is in nmol/min/mg protein (Table 2). The model files keep the printed pmol value in ini() and apply the 1000-fold conversion in model(). Without the conversion the minor M7 pathway would exceed the major M2 pathway, which would contradict the Discussion’s statement that CL_M2 is “at least fifty times larger than CL_M7”.
  • Km2 and Vmax2 for M2 are printed as “\” for mouse, rat and human. Those species needed only one Michaelis-Menten site. The second site is retained in the shared model structure with vmax2_m2 = fixed(0) and a placeholder km2_m2 = fixed(1), which removes it exactly.
  • The molecular weight is not in the paper. 398.4 g/mol for C22H22O7 (PubChem CID 73435) is used as a unit bridge between the mass-concentration ODE scale and the uM scale of Table 2’s Michaelis constants. It is a chemical constant, not a fitted parameter.
  • No residual-error or variance estimates are reported. The visual predictive checks (Figure 6) fitted a proportional residual error and exponential inter-individual variability on hepatic blood flow and metabolic velocity in Phoenix NLME, but the paper reports none of the resulting numbers, so propSd is fixed(0) and there are no etas. Figure 6 therefore cannot be reproduced from the packaged models; everything else in the paper can.
  • Rat in-vivo data are not in this paper. They are cited from Liu et al. 2016, which is not on disk, so the rat population metadata is sparser than the other species’.

Structural interpretations

  • The brain’s venous return. The printed venous-blood equation sums Qt * Ct / (Kt:pl / Rbp) over tissues, but the refined brain equation writes its perfusion term as Qbra * (Cart - C1), so blood leaving the brain carries the vascular concentration C1. Mass balance forces the venous pool to receive Qbra * C1 rather than the generic term, and that is how it is encoded. The mass-balance check above would fail otherwise.
  • The hepatic artery is derived, not printed. Table 1’s liver blood flow is the total (footnote h), so q_hepatic_artery = q_liver - q_gut - q_spleen.
  • Dosing compartment. The intravenous bolus is placed in the venous blood pool, consistent with the Figure 2 schematic. The paper does not state the dosing compartment explicitly.
  • Body weights used for dosing. Table 1’s reference weights (0.02, 0.25, 4, 8.5 and 70 kg) are used to turn the mg/kg doses into absolute amounts, since Table 1’s physiology is tabulated for those weights. The monkey and dog studies report slightly different mean weights (3.51 and 8.65 kg).

Where the reproduction is least tight

Half-life is reproduced within 10% and the fold-error never exceeds 1.60 across all 60 parameter-by-arm comparisons against Table 3’s predicted column, against the paper’s own criterion of two. The dog and human arms reproduce AUC within 5%; the mouse and rat arms sit around 20% high on AUC and correspondingly low on clearance and mean residence time.

That gradient tracks how much of the profile is the bolus mixing spike. An intravenous bolus is delivered into a venous blood pool whose contents turn over in V_venous / cardiac_output minutes – 0.08 min in the mouse, 0.16 in the rat, but 0.62 in the human. The area under that spike is a real part of the AUC in the mouse and rat and a negligible one in the dog and human, and it is exactly the part of the curve whose computed area depends most on the numerical integration grid and on the NCA settings. The paper reports neither the output grid it solved on nor the NCA options it applied to its predicted profiles, so this residual difference is not resolvable from the source. No parameter has been adjusted to close it.