Skip to contents

Oseltamivir in the ferret (Reddy 2015)

Source: vignettes/articles/Reddy_2015_oseltamivir_ferret.Rmd
Reddy_2015_oseltamivir_ferret.Rmd

Model and source 

mod <- readModelDb("Reddy_2015_oseltamivir_ferret")
ui  <- rxode2::rxode(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
  • Citation: Reddy MB, Yang K-H, Rao G, Rayner CR, Nie J, Pamulapati C, Marathe BM, Forrest A, Govorkova EA (2015). Oseltamivir Population Pharmacokinetics in the Ferret: Model Application for Pharmacokinetic/Pharmacodynamic Study Design. PLoS ONE 10(10):e0138069. doi:10.1371/journal.pone.0138069.
  • Description: Preclinical (ferret). Population PK model for oseltamivir carboxylate (OC), the active metabolite of oseltamivir, after oral dosing of oseltamivir phosphate or oseltamivir free base: an absorptive compartment feeding two transit compartments that model the delayed appearance of OC, followed by two-compartment OC disposition with first-order elimination. Clearance and volume terms are apparent values conditioned on oral bioavailability (F) and the fraction of parent converted to metabolite (Fm), and are normalised per kg of ferret body weight (dose is entered as ug of oseltamivir free base per kg). Parameter values are the pooled influenza A + B (ketamine-anaesthetised) Monte Carlo parameter set the authors used for all reported simulations.
  • Article (open access): https://doi.org/10.1371/journal.pone.0138069
  • Supporting information (PLOS ONE, open access): S1 Fig (individual model fits) and S1 Table (the modelling dataset, 451 oseltamivir-carboxylate concentrations in 65 ferrets).

Reddy 2015 is a meta-analysis of four ferret oseltamivir PK studies run at three sites. Its purpose is to make ferret influenza PK/PD studies designable: it produces a population PK model for the active metabolite oseltamivir carboxylate (OC), tests whether anaesthesia or influenza inoculation perturb that PK, and then uses the model to pick a ferret dose that matches human steady-state exposure.

Population

65 ferrets from 3 studies contributed 430 OC concentrations above the 10 ng/mL limit of quantitation (Reddy 2015 Methods, “OC PK Model”). The animals were young adults, 3-5 months (young adult), weighing 0.6-0.92 kg. Studies 1 (Beijing) and 2 (Memphis) used male ferrets and ketamine anaesthesia before blood sampling; Study 3 (London) used female ferrets and no anaesthesia, so the female fraction of the modelled cohort is 8/65. Study 4 (Cardiff), in which animals were maintained under Saffan anaesthesia for 12 h, was excluded from model development because its NCA profile differed materially from the others (Tmax 7 h for both prodrug and metabolite versus 1 h and 3-4 h elsewhere; Reddy 2015 Table 1 and Results, “Population PK Model”).

Seventeen animals were uninfected; the remainder were inoculated with influenza A/Shenzheng/406H/2006 (H5N1, n = 18), A/Hong Kong/433581/2009 (H3N2, n = 12) or B/Yamagata/16/1988 (n = 18). All three inoculations produced only mild, essentially subclinical illness, which the authors flag as the principal limitation on their “influenza does not change OC PK” conclusion.

str(ui$population)
#> List of 11
#>  $ species       : chr "ferret (Mustela putorius furo)"
#>  $ n_subjects    : num 65
#>  $ n_studies     : num 3
#>  $ n_observations: num 430
#>  $ age_range     : chr "3-5 months (young adult)"
#>  $ weight_range  : chr "0.6-0.92 kg"
#>  $ sex_female_pct: num 12.3
#>  $ disease_state : chr "uninfected (n = 17) or inoculated with influenza A/Shenzheng/406H/2006 (H5N1, n = 18), influenza A/Hong Kong/43"| __truncated__
#>  $ dose_range    : chr "0.76-25 mg/kg oseltamivir free base (equivalently 1.0-32.9 mg/kg oseltamivir phosphate) orally, as single doses"| __truncated__
#>  $ regions       : chr "Beijing (China), Memphis (USA), London (UK)"
#>  $ notes         : chr "Reddy 2015 Methods 'PK Studies in a Ferret Model' and 'OC PK Model'. Studies 1 and 2 used male ferrets, Study 3"| __truncated__

Unit system

The modelling dataset (S1 Table) carries dose as micrograms of oseltamivir free base per kg of body weight and OC concentration as ng/mL, and contains no body-weight column. The per-kg dose numbers are therefore the amounts the model was fitted against, so Vc, Vp, CLd and CLt – printed as L and L/h in Reddy 2015 Tables 4 and 6 – are per-kg values, and Cc = central / vc returns ug/L, which is numerically ng/mL. This vignette doses in ug of free base throughout.

The dataset also settles which salt form the dose column refers to. Study 2 administered oseltamivir phosphate at 1.0, 5.0 and 25.0 mg/kg, which the paper equates to 0.76, 3.8 and 19 mg/kg free base; the dataset’s dose values for that study are 760, 3800 and 19000 ug/kg, i.e. free base. Study 3’s 5.0 and 25.0 mg/kg phosphate doses likewise appear as 3800 and 19000 ug/kg. The conversions used below follow the paper: 1 mg/kg oseltamivir phosphate = 0.7616 mg/kg free base.

# Free-base equivalents exactly as Reddy 2015 states them: Study 2 Methods
# ("OP doses of 1.0, 5.0, or 25.0 mg/kg ... i.e., a single 0.76, 3.8, and
# 19 mg/kg dose of OFB") and Methods, "Simulations" ("3.87 mg/kg of OFB
# (5.08 mg/kg of OP)").
ofb_ug <- c(`1` = 760, `5` = 3800, `25` = 19000, `5.08` = 3870)
dose_ug <- function(op_mg_per_kg) {
  key <- as.character(op_mg_per_kg)
  stopifnot(key %in% names(ofb_ug))
  unname(ofb_ug[key])
}
ofb_ug
#>     1     5    25  5.08 
#>   760  3800 19000  3870

Source trace

Per-parameter origins are recorded as in-file comments beside each ini() entry in inst/modeldb/specificDrugs/Reddy_2015_oseltamivir_ferret.R. They are collected here for review. Reddy 2015 Table 6 is the parameter set the authors carried into every simulation they report; its means are identical to the “Studies 1 and 2, ketamine” column of Table 4.

Equation / parameter Value Source location
d/dt(depot) = -ktr * depot n/a Eq (1), p. 6
d/dt(transit1) = ktr * depot - ktr * transit1 n/a Eq (2), p. 6
d/dt(transit2) = ktr * transit1 - ka * transit2 n/a Eq (3), p. 6
d/dt(central) = ka * transit2 - q/vc * central + q/vp * peripheral1 - cl/vc * central n/a Eq (4), p. 6
d/dt(peripheral1) = q/vc * central - q/vp * peripheral1 n/a Eq (5), p. 6
Compartment roles (depot / 2 transit / central / peripheral) n/a Fig 1 schematic
Cc <- central / vc n/a Methods, “The concentration of OC in the plasma is calculated as X4/Vc”
lktr (Kt) 1.27 1/h Table 6, Kt mean
lka (Ka) 0.463 1/h Table 6, Ka mean
lq (CLd) 0.585 L/h Table 6, CLd mean
lcl (CLt) 1.52 L/h Table 6, CLt mean
lvc (Vc) 0.157 L Table 6, Vc mean
lvp (Vp) 5.59 L Table 6, Vp mean
etalktr var 0.402 Table 6 Kt variance; CV% imposed at 50 per Methods, “Simulations”
etalq var 0.0856 Table 6 CLd variance; CV% imposed at 50
etalvc var 0.00615 Table 6 Vc variance; CV% imposed at 50
etalka, etalcl, etalvp block var 0.0937, 0.704, 12.6 Table 6 diagonals
block covariances 0.104, -0.636, -0.863 Table 6 off-diagonals CLt, Ka; Vp, Ka; Vp, CLt
addSd 5 ng/mL Results, “The intercept for the additive error model was 5 ng/mL”
propSd 0.15 Results, “The slope was 0.15”

Reddy 2015 reports the variances and covariances on the natural parameter scale while stating that the parameters are log-normally distributed. The model file stores the moment-matched log-scale block; the chunk below confirms that it round-trips to the published natural-scale CV% and correlations.

omega <- as.matrix(ui$omega)
attr(omega, "lotriFix") <- NULL   # drop the fixed-flag attribute so prints stay readable
lvar  <- diag(omega)
nat_cv <- sqrt(exp(lvar) - 1)
published_cv <- c(etalktr = 0.499, etalq = 0.500, etalvc = 0.499,
                  etalka = 0.661, etalcl = 0.552, etalvp = 0.635)
round_trip <- data.frame(
  eta = names(lvar),
  `Published CV` = round(100 * published_cv[names(lvar)], 1),
  `Model CV` = round(100 * nat_cv, 1),
  check.names = FALSE
)
knitr::kable(round_trip, row.names = FALSE, caption = "Natural-scale CV% implied by the stored log-scale omega block versus Reddy 2015 Tables 4 and 6.")
Natural-scale CV% implied by the stored log-scale omega block versus Reddy 2015 Tables 4 and 6.
eta Published CV Model CV
etalktr 49.9 49.9
etalq 50.0 50.0
etalvc 49.9 50.0
etalka 66.1 66.1
etalcl 55.2 55.2
etalvp 63.5 63.5 

# Deterministic: a transcription error in any variance moves these by >1 point.
stopifnot(max(abs(100 * (nat_cv - published_cv[names(lvar)]))) < 0.15)

# Natural-scale correlations, Table 6 off-diagonals divided by the SDs.
m  <- c(ka = 0.463, cl = 1.52, vp = 5.59)
nv <- c(ka = 0.0937, cl = 0.704, vp = 12.6)
published_corr <- c(
  `cl~ka` = 0.104 / sqrt(nv[["cl"]] * nv[["ka"]]),
  `vp~ka` = -0.636 / sqrt(nv[["vp"]] * nv[["ka"]]),
  `vp~cl` = -0.863 / sqrt(nv[["vp"]] * nv[["cl"]])
)
lognormal_corr <- function(i, j) {
  cv_i <- sqrt(exp(log(1 + nv[[i]] / m[[i]]^2)) - 1)
  cv_j <- sqrt(exp(log(1 + nv[[j]] / m[[j]]^2)) - 1)
  (exp(omega[paste0("etal", i), paste0("etal", j)]) - 1) / (cv_i * cv_j)
}
recovered <- c(`cl~ka` = lognormal_corr("cl", "ka"),
               `vp~ka` = lognormal_corr("vp", "ka"),
               `vp~cl` = lognormal_corr("vp", "cl"))
stopifnot(max(abs(recovered - published_corr)) < 1e-6)
knitr::kable(
  data.frame(Pair = names(published_corr),
             `Table 6 correlation` = round(published_corr, 3),
             `Recovered from omega` = round(recovered, 3),
             check.names = FALSE, row.names = NULL),
  caption = "Natural-scale correlations recovered from the stored log-scale covariances."
)
Natural-scale correlations recovered from the stored log-scale covariances.
Pair Table 6 correlation Recovered from omega
cl~ka 0.405 0.405
vp~ka -0.585 -0.585
vp~cl -0.290 -0.290

The omega block must also be positive definite, or rxSolve cannot draw a cohort from it.

stopifnot(all(eigen(omega, symmetric = TRUE, only.values = TRUE)$values > 0))
round(eigen(omega, symmetric = TRUE, only.values = TRUE)$values, 4)
#> [1] 0.7021 0.2232 0.2227 0.2225 0.1990 0.0662

Typical-value structure checks

Steady-state exposure of a linear model is dose over clearance, and Reddy 2015 computes AUCss24h exactly that way (“the 24-h steady-state AUC was calculated by taking the 24-h dose divided by the fitted CLt”). That makes it a closed-form gate on the packaged ODEs: if a transfer rate were mis-signed or a distribution term dropped, the solved AUC would not match.

mod_typ <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'

solve_typical <- function(amt, ii = 12, addl = 9, obs) {
  ev <- rxode2::et(amt = amt, cmt = "depot", ii = ii, addl = addl) |>
    rxode2::et(obs)
  out <- as.data.frame(rxode2::rxSolve(mod_typ, ev))
  if (is.null(out$id)) out$id <- 1L
  out
}
trapz <- function(t, y) sum(diff(t) * (head(y, -1) + tail(y, -1)) / 2)

ss <- solve_typical(dose_ug(5), obs = seq(108, 120, by = 0.02))
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'
auc_ss12 <- trapz(ss$time, ss$Cc)
cl_typ   <- exp(ui$theta[["lcl"]])

c(`AUCss12h solved (ng*h/mL)` = round(auc_ss12, 1),
  `Dose / CLt (ng*h/mL)`      = round(dose_ug(5) / cl_typ, 1),
  `% difference`              = round(100 * (auc_ss12 / (dose_ug(5) / cl_typ) - 1), 3))
#> AUCss12h solved (ng*h/mL)      Dose / CLt (ng*h/mL)              % difference 
#>                  2499.900                  2500.000                    -0.004

# Deterministic solve against its own closed form: tight bound is correct here.
stopifnot(abs(auc_ss12 / (dose_ug(5) / cl_typ) - 1) < 0.002)

# Linearity: the paper assumes linear PK, so 25 times the dose gives 25 times
# the AUC (the free-base doses are 760 and 19000 ug, an exact 25-fold ratio).
ss_lo  <- solve_typical(dose_ug(1),  obs = seq(108, 120, by = 0.02))
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'
ss_hi  <- solve_typical(dose_ug(25), obs = seq(108, 120, by = 0.02))
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'
auc_lo <- trapz(ss_lo$time, ss_lo$Cc)
auc_hi <- trapz(ss_hi$time, ss_hi$Cc)
stopifnot(abs(auc_hi / auc_lo / 25 - 1) < 0.002)
c(`AUC ratio, OP 25 vs 1 mg/kg` = round(auc_hi / auc_lo, 3))
#> AUC ratio, OP 25 vs 1 mg/kg 
#>                          25

Replicating Table 1: observed single-dose NCA

Reddy 2015 Table 1 reports non-compartmental parameters for uninfected ferrets after a single oral dose. The Study 1 row (ketamine anaesthesia, 5.0 mg/kg free base, n = 3, rich sampling) is the row this model should reproduce: it is the same anaesthesia stratum the Table 6 parameters come from, and NCA values were never used to fit the model, so the comparison is not circular.

sd_obs <- seq(0, 24, by = 0.02)
sd_sim <- solve_typical(amt = 5000, ii = 12, addl = 0, obs = sd_obs) |>
  mutate(treatment = "OFB 5.0 mg/kg, single dose")
#> Warning: 'ii' requires non zero additional doses ('addl') or steady state
#> dosing ('ii': 12.000000, 'ss': 0; 'addl': 0), reset 'ii' to zero
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'

stopifnot(all(is.finite(sd_sim$Cc)), all(sd_sim$Cc >= 0))

sd_conc <- sd_sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)
sd_dose <- data.frame(id = 1L, time = 0, amt = 5000,
                      treatment = "OFB 5.0 mg/kg, single dose")

sd_nca <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(sd_conc, Cc ~ time | treatment + id),
  PKNCA::PKNCAdose(sd_dose, amt ~ time | treatment + id),
  intervals = data.frame(start = 0, end = 12,
                         cmax = TRUE, tmax = TRUE, auclast = TRUE)
))
published_t1 <- tibble::tribble(
  ~treatment,                     ~cmax, ~tmax, ~auclast,
  "OFB 5.0 mg/kg, single dose",   450,   3,     2610
)

cmp_t1 <- nlmixr2lib::ncaComparisonTable(
  simulated = sd_nca, reference = published_t1, by = "treatment",
  units = c(cmax = "ng/mL", tmax = "h", auclast = "ng*h/mL"),
  tolerance_pct = 20
)
knitr::kable(
  cmp_t1,
  caption = paste(
    "Typical-value prediction versus Reddy 2015 Table 1, Study 1 (uninfected,",
    "ketamine, 5.0 mg/kg oseltamivir free base, n = 3). Published values are",
    "arithmetic means; the reported SDs are 0.47 ug/mL on Cmax, 2 h on Tmax",
    "and 1.99 ug*h/mL on AUC0-12h, so all three simulated values sit well",
    "inside one observed SD. * differs from reference by >20%."
  ),
  align = c("l", "l", "r", "r", "r")
)
Typical-value prediction versus Reddy 2015 Table 1, Study 1 (uninfected, ketamine, 5.0 mg/kg oseltamivir free base, n = 3). Published values are arithmetic means; the reported SDs are 0.47 ug/mL on Cmax, 2 h on Tmax and 1.99 ugh/mL on AUC0-12h, so all three simulated values sit well inside one observed SD. differs from reference by >20%.
NCA parameter treatment Reference Simulated % diff
Cmax (ng/mL) OFB 5.0 mg/kg, single dose 450 539 +19.8%
Tmax (h) OFB 5.0 mg/kg, single dose 3 2.34 -22.0%*
AUClast (ng*h/mL) OFB 5.0 mg/kg, single dose 2610 2760 +5.9% 
attr(cmp_t1, "footnote")
#> [1] "* differs from reference by more than ±20%."

t1_pct <- as.numeric(gsub("[^0-9.-]", "", cmp_t1[["% diff"]]))
names(t1_pct) <- cmp_t1[[1]]
# Deterministic typical-value solve vs a published mean of n = 3 animals.
# Realised: Cmax +19.8%, Tmax -21.7%, AUC0-12h +5.9%. A mis-transcribed
# clearance, volume or dose moves these by tens of percent.
stopifnot(max(abs(t1_pct)) < 35)

Cmax runs high and Tmax early against a three-animal mean whose own SD on Cmax is larger than its mean; AUC0-12h, the parameter least sensitive to the shape of the absorption phase, agrees to 6%.

Replicating Table 7: Monte Carlo steady-state exposure

Table 7 is the paper’s own Monte Carlo output using exactly the Table 6 parameter set: 1000 ferrets dosed every 12 h for 5 days at oseltamivir phosphate doses of 1.0, 5.0 and 25.0 mg/kg. Reproducing it exercises the between-animal variability block as well as the structural model.

# rxSetSeed() fixes rxode2's stream per solver thread, not across thread
# counts, so this cohort differs between a workstation and a CI runner.
# Every assertion below is written to hold for any cohort the model can draw.
set.seed(20151013)
rxode2::rxSetSeed(20151013)

n_per_arm <- 200L   # skill cap: never more than 200 participants per arm

make_arm <- function(n, op_mg_per_kg, label, id_offset) {
  ids <- id_offset + seq_len(n)
  amt <- dose_ug(op_mg_per_kg)
  dose <- tidyr::crossing(id = ids, time = seq(0, 108, by = 12)) |>
    mutate(amt = amt, evid = 1L, cmt = "depot")
  obs <- tidyr::crossing(id = ids, time = seq(108, 120, by = 0.05)) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "central")
  bind_rows(dose, obs) |>
    mutate(treatment = label) |>
    arrange(id, time, desc(evid))
}

arms <- tibble::tribble(
  ~op,   ~label,
  1.0,   "OP 1.0 mg/kg q12h",
  5.0,   "OP 5.0 mg/kg q12h",
  25.0,  "OP 25.0 mg/kg q12h"
)

events <- bind_rows(lapply(seq_len(nrow(arms)), function(i) {
  make_arm(n_per_arm, arms$op[i], arms$label[i], id_offset = (i - 1L) * n_per_arm)
}))
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
sim <- as.data.frame(rxode2::rxSolve(mod, events = events, keep = "treatment"))
#> ℹ parameter labels from comments will be replaced by 'label()'
stopifnot(all(is.finite(sim$Cc)), all(sim$Cc >= 0), nrow(sim) > 0)

# Recentre on the final dosing interval so PKNCA sees a clean 0-12 h window
# whose t = 0 record is the steady-state trough.
ss_conc <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  mutate(time = time - 108) |>
  dplyr::select(id, time, Cc, treatment)

ss_dose <- events |>
  dplyr::filter(evid == 1, time == 108) |>
  mutate(time = 0) |>
  dplyr::select(id, time, amt, treatment)

ss_nca <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(ss_conc, Cc ~ time | treatment + id),
  PKNCA::PKNCAdose(ss_dose, amt ~ time | treatment + id),
  intervals = data.frame(start = 0, end = 12,
                         cmax = TRUE, tmax = TRUE, cmin = TRUE, auclast = TRUE)
))

Table 7 prints AUCss24h; because the paper defines it as the 24-h dose over clearance, it is exactly twice the 12-h interval AUC, and the reference column below is Table 7’s value halved. Table 7’s Cmax and Cmin columns are headed ug/L, but the values are ug/mL – see Errata – so they are multiplied by 1000 to reach ng/mL.

published_t7 <- tibble::tribble(
  ~treatment,              ~tmax, ~cmax,   ~cmin, ~auclast,
  "OP 1.0 mg/kg q12h",     2.6,   98,      18,    1060 / 2,
  "OP 5.0 mg/kg q12h",     2.6,   494,     90,    5310 / 2,
  "OP 25.0 mg/kg q12h",    2.6,   2470,    460,   26600 / 2
)

cmp_t7 <- nlmixr2lib::ncaComparisonTable(
  simulated = ss_nca, reference = published_t7, by = "treatment",
  units = c(cmax = "ng/mL", cmin = "ng/mL", tmax = "h", auclast = "ng*h/mL"),
  tolerance_pct = 20
)
knitr::kable(
  cmp_t7,
  caption = paste(
    "Cohort medians (n = 200 per arm) versus the medians of Reddy 2015",
    "Table 7 (1000 ferrets). * differs from reference by >20%."
  ),
  align = c("l", "l", "r", "r", "r")
)
Cohort medians (n = 200 per arm) versus the medians of Reddy 2015 Table 7 (1000 ferrets). * differs from reference by >20%.
NCA parameter treatment Reference Simulated % diff
Cmax (ng/mL) OP 1.0 mg/kg q12h 98 91.1 -7.0%
Cmax (ng/mL) OP 5.0 mg/kg q12h 494 424 -14.2%
Cmax (ng/mL) OP 25.0 mg/kg q12h 2470 2070 -16.1%
Cmin (ng/mL) OP 1.0 mg/kg q12h 18 11.6 -35.6%*
Cmin (ng/mL) OP 5.0 mg/kg q12h 90 63.1 -29.9%*
Cmin (ng/mL) OP 25.0 mg/kg q12h 460 286 -37.8%*
Tmax (h) OP 1.0 mg/kg q12h 2.6 2.2 -15.4%
Tmax (h) OP 5.0 mg/kg q12h 2.6 2.35 -9.6%
Tmax (h) OP 25.0 mg/kg q12h 2.6 2.4 -7.7%
AUClast (ng*h/mL) OP 1.0 mg/kg q12h 530 493 -7.1%
AUClast (ng*h/mL) OP 5.0 mg/kg q12h 2660 2500 -5.9%
AUClast (ng*h/mL) OP 25.0 mg/kg q12h 13300 11900 -10.8% 
attr(cmp_t7, "footnote")
#> [1] "* differs from reference by more than ±20%."
t7 <- cmp_t7
t7$pct  <- as.numeric(gsub("[^0-9.-]", "", t7[["% diff"]]))
t7$code <- ifelse(grepl("^Cmax", t7[[1]]), "cmax",
           ifelse(grepl("^Cmin", t7[[1]]), "cmin",
           ifelse(grepl("^Tmax", t7[[1]]), "tmax", "auclast")))
# Guard against a label-mismatch silently emptying the gate (pattern 10).
stopifnot(nrow(t7) == 12L, sum(t7$code == "cmin") == 3L, !anyNA(t7$pct))

# Cmin is excluded from the gate as a documented deviation (see Errata): it has
# a published CV of 100%, so its median is the least well determined quantity
# in Table 7 and the model reproducibly sits about a third below it.
gated <- t7[t7$code != "cmin", ]
# Realised across arms: Cmax -7%, Tmax -13%, AUC -4%. Cohort-derived, so the
# bound must admit the between-run spread; 30 still breaks on a mis-transcribed
# clearance, volume or dose, all of which move these by tens of percent.
stopifnot(max(abs(gated$pct)) < 30)
knitr::kable(t7[, c(1, 2, ncol(t7) - 1, ncol(t7))],
             row.names = FALSE,
             caption = "Percent differences by parameter; cmin is a documented deviation.")
Percent differences by parameter; cmin is a documented deviation.
NCA parameter treatment pct code
Cmax (ng/mL) OP 1.0 mg/kg q12h -7.0 cmax
Cmax (ng/mL) OP 5.0 mg/kg q12h -14.2 cmax
Cmax (ng/mL) OP 25.0 mg/kg q12h -16.1 cmax
Cmin (ng/mL) OP 1.0 mg/kg q12h -35.6 cmin
Cmin (ng/mL) OP 5.0 mg/kg q12h -29.9 cmin
Cmin (ng/mL) OP 25.0 mg/kg q12h -37.8 cmin
Tmax (h) OP 1.0 mg/kg q12h -15.4 tmax
Tmax (h) OP 5.0 mg/kg q12h -9.6 tmax
Tmax (h) OP 25.0 mg/kg q12h -7.7 tmax
AUClast (ng*h/mL) OP 1.0 mg/kg q12h -7.1 auclast
AUClast (ng*h/mL) OP 5.0 mg/kg q12h -5.9 auclast
AUClast (ng*h/mL) OP 25.0 mg/kg q12h -10.8 auclast

Replicating Figure 6

Figure 6 shows the Monte Carlo steady-state OC profile at the 5.08 mg/kg oseltamivir phosphate dose, with the median and the 10th-90th percentile band.

rxode2::rxSetSeed(20151013)
fig6_events <- make_arm(n_per_arm, 5.08, "OP 5.08 mg/kg q12h", id_offset = 0L)
fig6 <- as.data.frame(rxode2::rxSolve(mod, events = fig6_events)) |>
  mutate(time = time - 108)

fig6 |>
  group_by(time) |>
  summarise(Q10 = quantile(Cc, 0.10), Q50 = median(Cc), Q90 = quantile(Cc, 0.90),
            .groups = "drop") |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.25) +
  geom_line(linewidth = 1) +
  labs(x = "Time after the last dose (h)", y = "OC concentration (ng/mL)",
       title = "Steady-state OC PK at 5.08 mg/kg oseltamivir phosphate q12h",
       caption = "Median with 10th-90th percentile band; replicates Figure 6 of Reddy 2015.")
Replicates Figure 6 of Reddy 2015: simulated steady-state OC concentrations at 5.08 mg/kg oseltamivir phosphate every 12 h.

Replicates Figure 6 of Reddy 2015: simulated steady-state OC concentrations at 5.08 mg/kg oseltamivir phosphate every 12 h. 

fig6_auc <- fig6 |>
  group_by(id) |>
  summarise(auc12 = trapz(time, Cc), .groups = "drop")

human_target <- 3220   # ng*h/mL, human AUCss12h at 75 mg BID (Reddy 2015 Discussion)
c(`Simulated median AUCss12h (ng*h/mL)` = round(median(fig6_auc$auc12)),
  `Human target (ng*h/mL)`              = human_target,
  `% difference`                        = round(100 * (median(fig6_auc$auc12) / human_target - 1), 1),
  `Dose reproducing the target (mg/kg OP)` =
    round(5.08 * human_target / median(fig6_auc$auc12), 2))
#>    Simulated median AUCss12h (ng*h/mL)                 Human target (ng*h/mL) 
#>                                 2518.0                                 3220.0 
#>                           % difference Dose reproducing the target (mg/kg OP) 
#>                                  -21.8                                    6.5

The paper’s headline dose recommendation does not reproduce from the paper’s own parameters: at 5.08 mg/kg the model returns a median AUCss12h roughly 20% below the 3220 ng*h/mL human target it is meant to match. This is recorded as a deviation rather than gated – see Errata.

Anaesthesia: reproducing Table 4

The paper’s second finding is that ketamine anaesthesia perturbs OC PK. It was established by a Kruskal-Wallis test on the individual post-hoc parameter estimates (Table 4), not by fitting an anaesthesia covariate, so the packaged model carries no covariate term; the no-anaesthesia parameter set is a re-parameterisation of the same structure. Applying it to Study 3’s dosing (3.8 mg/kg free base, no anaesthesia) tests the finding against Table 1’s observed NCA for that study.

mod_noanes <- suppressMessages(rxode2::ini(
  mod_typ,
  lktr = log(4.18), lka = log(0.335), lq = log(0.878),
  lcl = log(0.919), lvc = log(1.00),  lvp = log(2.08)
))

profile_for <- function(model, amt) {
  ev <- rxode2::et(amt = amt, cmt = "depot") |> rxode2::et(seq(0, 24, by = 0.02))
  out <- as.data.frame(rxode2::rxSolve(model, ev))
  out[out$time <= 12, ]
}

s3_ket <- profile_for(mod_typ,    3800)
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'
s3_non <- profile_for(mod_noanes, 3800)
#> ℹ omega/sigma items treated as zero: 'etalktr', 'etalq', 'etalvc', 'etalka', 'etalcl', 'etalvp'

study3 <- data.frame(
  Parameter = c("Cmax (ng/mL)", "Tmax (h)", "AUC0-12h (ng*h/mL)"),
  `Observed, Study 3` = c(580, 3, 4030),
  `Ketamine parameters` = c(round(max(s3_ket$Cc)),
                            s3_ket$time[which.max(s3_ket$Cc)],
                            round(trapz(s3_ket$time, s3_ket$Cc))),
  `No-anaesthesia parameters` = c(round(max(s3_non$Cc)),
                                  s3_non$time[which.max(s3_non$Cc)],
                                  round(trapz(s3_non$time, s3_non$Cc))),
  check.names = FALSE
)
knitr::kable(study3, digits = 0, row.names = FALSE, caption = paste(
  "Reddy 2015 Table 1 Study 3 (no anaesthesia, 3.8 mg/kg free base, n = 4)",
  "against typical-value predictions from the two Table 4 parameter sets."
))
Reddy 2015 Table 1 Study 3 (no anaesthesia, 3.8 mg/kg free base, n = 4) against typical-value predictions from the two Table 4 parameter sets.
Parameter Observed, Study 3 Ketamine parameters No-anaesthesia parameters
Cmax (ng/mL) 580 410 518
Tmax (h) 3 2 2
AUC0-12h (ng*h/mL) 4030 2100 3520 

auc_ket <- trapz(s3_ket$time, s3_ket$Cc)
auc_non <- trapz(s3_non$time, s3_non$Cc)
# Both solves are deterministic typical values, so comparing them is safe.
# Realised: no-anaesthesia -7.7%, ketamine -47.9% against the observed 4030.
stopifnot(abs(auc_non / 4030 - 1) < 0.25)
stopifnot(abs(auc_non / 4030 - 1) < abs(auc_ket / 4030 - 1))

The no-anaesthesia parameter set predicts the unanaesthetised Study 3 exposure to within 8%, while the ketamine set under-predicts it by roughly half. The packaged model is the ketamine set, and should be re-parameterised as above before being applied to unanaesthetised ferrets.

bind_rows(
  s3_ket |> mutate(Parameters = "Studies 1-2, ketamine"),
  s3_non |> mutate(Parameters = "Study 3, no anaesthesia")
) |>
  ggplot(aes(time, Cc, colour = Parameters)) +
  geom_line(linewidth = 1) +
  labs(x = "Time (h)", y = "OC concentration (ng/mL)",
       title = "Anaesthesia effect on typical-value OC PK") +
  theme(legend.position = "bottom")
Typical-value OC profiles after 3.8 mg/kg oseltamivir free base under the two Reddy 2015 Table 4 parameter sets.

Typical-value OC profiles after 3.8 mg/kg oseltamivir free base under the two Reddy 2015 Table 4 parameter sets.

Assumptions and deviations

  • Log-normal moment matching. Reddy 2015 states that parameters are log-normally distributed but reports means, variances and covariances on the natural scale (Table 6). The model stores omega2 = log(1 + var / mean^2) and cov_log(i,j) = log(1 + cov(i,j) / (mean(i) * mean(j))), so exp(l<param>) is the median of the log-normal and the natural-scale CV% and correlations round-trip exactly (verified above). Interpreting the printed means as medians rather than arithmetic means is the reading that reproduces the paper’s own Table 7 medians; treating them as arithmetic means would shift every exposure by a further 14%.
  • Per-kg parameters. The dose column of the modelling dataset is ug of oseltamivir free base per kg and there is no body-weight column, so the clearances and volumes printed as L/h and L are per-kg values. The model file and this vignette dose in ug of free base accordingly. Users with an absolute-amount dataset must divide the dose by body weight in kg.
  • fixed() on three variances. Methods, “Simulations” states that the CV% for Vc, CLd and Kt “were empirically reduced to 50%, and all related covariance terms were fixed to zero”. Those three variances are therefore wrapped in fixed() and sit outside the correlated block; the estimated values were 165%, 96.6% and 79.8% CV (Table 4). Covariances whose correlation was below 0.15 were also fixed to zero by the authors and are absent here.
  • No covariate effects. Anaesthesia and influenza inoculation strain were screened by Kruskal-Wallis tests on post-hoc parameter estimates rather than modelled as covariate coefficients, so neither appears in model(). Ketamine is documented under covariatesDataExcluded. The no-anaesthesia parameter set (Table 4, n = 8, no covariance matrix reported, never used for simulation) is reproduced above as a re-parameterisation rather than as a second model file.
  • Study 4 excluded upstream. Saffan-anaesthetised animals were left out of model development by the authors, so the packaged model makes no claim about them.
  • All parameter values come from the paper’s own tables and text. Nothing was digitised from a figure, supplied by correspondence, or carried from an upstream model.

Errata

  • Table 7 concentration units. The Cmax and Cmin columns of Table 7 are headed ug/L, but the values must be ug/mL. Internal evidence: the mean AUCss24h at 5.0 mg/kg is 5990 ug*h/L, so Cavg is 250 ug/L = 0.25 ug/mL, which lies between the printed Cmin 0.134 and Cmax 0.544 only if those two are read as ug/mL. Read as ug/L they would fall three orders of magnitude below Cavg. The AUCss24h column heading ug*h/L is correct. This vignette multiplies the Table 7 Cmax and Cmin values by 1000 to reach ng/mL.
  • The 5.08 mg/kg dose recommendation does not reproduce. Reddy 2015 concludes that 5.08 mg/kg oseltamivir phosphate every 12 h achieves the human steady-state AUCss12h of 3220 ngh/mL. The arithmetic the paper itself prescribes contradicts this: AUCss12h = dose / CLt gives 3870 ug / 1.52 L/h = 2546 ngh/mL, 21% below the target, and matching 3220 ng*h/mL would need 6.42 mg/kg oseltamivir phosphate. The simulated cohort above agrees with that closed form. Table 7’s own medians are consistent with the model to within 6%, so the discrepancy sits in the separate 100-ferret simulation that produced the 5.08 figure, not in the parameter set. No parameter was adjusted; the deviation is reported as measured.
  • Dose units in Methods, “Simulations”. The text gives the simulated doses as “1.0, 5.0, and 25.0 mg of oral OP”, while Table 7 heads the same column “OP dose (mg/kg)”. The dataset and the reproduced exposures both support mg/kg.
  • Equation (4) sign. pdftotext drops the operator glyphs in equations (1) to (5). The publisher’s equation images (pone.0138069.e001 to e005) were read directly to confirm the signs, in particular the + CLd * X5 / Vp return term in equation (4).
  • No erratum or corrigendum to this article was found on the PLOS ONE article page or in PubMed.