Skip to contents

Morphine EEG effect in rat (Groenendaal 2007)

Source: vignettes/articles/Groenendaal_2007_morphine_rat_eeg.Rmd
Groenendaal_2007_morphine_rat_eeg.Rmd

Model and source

  • Citation: Groenendaal D, Freijer J, de Mik D, Bouw MR, Danhof M, de Lange ECM. Influence of biophase distribution and P-glycoprotein interaction on pharmacokinetic-pharmacodynamic modelling of the effects of morphine on the EEG. Br J Pharmacol. 2007;151(6):713-720. doi:10.1038/sj.bjp.0707258. Blood pharmacokinetic parameters fixed from the companion paper: Groenendaal D, Freijer J, de Mik D, Bouw MR, Danhof M, de Lange ECM. Population pharmacokinetic modelling of non-linear brain distribution of morphine: influence of active saturable influx and P-glycoprotein mediated efflux. Br J Pharmacol. 2007;151(5):701-712; doi:10.1038/sj.bjp.0707257, Table 2.
  • Description: Preclinical (rat, male Wistar). Population PK-PD model for the effect of morphine on the amplitude of the delta frequency band (0.5-4.5 Hz) of the rat EEG, with an extended catenary biophase distribution model and a P-glycoprotein (Pgp) interaction at the blood-brain barrier (Groenendaal 2007, Br J Pharmacol 151(6):713-720). Blood disposition is a three-compartment model that serves purely as the input function; its parameters were not re-estimated here and are fixed from the companion paper (Groenendaal 2007, Br J Pharmacol 151(5):701-712, Table 2), which carries linear body-weight covariate effects on CL and V2. Biophase distribution is the ‘tank-in-series’ extended catenary model of Upton 2000: two sequential biophase compartments holding CONCENTRATIONS, a transfer compartment (effect1, paper symbol Cet) fed from blood at rate k1e and drained at the same rate k1e, followed by the effect compartment (effect2, paper symbol Ce) drained at rate keo. The asymmetric form (k1e != keo) was retained over the symmetric form on objective function (24671 vs 24936); the simpler one-compartment effect-site model could be fit in neither its symmetric nor its asymmetric form. Co-infusion of the Pgp inhibitor GF120918 (elacridar) reduces keo by 64 percent and leaves k1e unchanged. Effect-compartment concentrations drive EEG amplitude through a sigmoidal Emax model. Inter-animal variability is exponential on keo and PROPORTIONAL on E0 and Emax; all other variances were fixed to zero by the authors. A covariate contrasting the EEG and EEG-microdialysis experimental methods was screened on E0 and Emax and not retained.
  • Article: https://doi.org/10.1038/sj.bjp.0707258
  • Companion article (blood PK input function): https://doi.org/10.1038/sj.bjp.0707257

Both papers are open access at Europe PMC (PMC2013992 and PMC2014001).

Population

Sixty-eight male Wistar rats (250-350 g, Charles River) received a 10-minute intravenous infusion of morphine hydrochloride at 4, 10 or 40 mg/kg. One arm (4 mg/kg) additionally received the P-glycoprotein inhibitor GF120918 (elacridar) as a 6 mg/kg one-minute bolus followed by a continuous 25 ng/min infusion begun 120 minutes before morphine; the mean steady-state GF120918 blood concentration was 214 ng/mL, sufficient to block P-gp. Every animal received a continuous midazolam infusion (5.5 mg/kg/h, Wagner loading scheme, mean steady-state 937 ng/mL) to suppress opioid-induced seizure activity, and animals in the 40 mg/kg groups were artificially ventilated and given vecuronium bromide for muscle rigidity.

The pharmacodynamic endpoint is the amplitude of the delta frequency band (0.5-4.5 Hz) of the EEG in uV, recorded continuously from bipolar leads on the left hemisphere and averaged into 3-minute bins to 75 minutes, 5-minute bins to 200 minutes and 10-minute bins to 360 minutes.

The design table is Table 1 of the companion paper, which covers both publications; the EEG paper’s model was fit to its experiments 1 (EEG only) and 2 (EEG plus microdialysis). Group sizes and mean body weights: EEG 4 mg/kg N = 7 (0.294 kg), EEG 10 mg/kg N = 7 (0.260 kg), EEG 40 mg/kg N = 5 (0.273 kg), EEG/MD 4 mg/kg N = 14 (0.297 kg), EEG/MD 4 mg/kg + GF120918 N = 20 (0.300 kg), EEG/MD 40 mg/kg N = 15 (0.306 kg). The microdialysis-only arm (N = 3) of the companion paper contributed no EEG data and is excluded here.

The same information is available programmatically via the model’s population metadata (readModelDb("Groenendaal_2007_morphine_rat_eeg")()$population).

Source trace

Per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Groenendaal_2007_morphine_rat_eeg.R. The table below collects them in one place. “Lead” is Groenendaal 2007 Br J Pharmacol 151(6):713-720 (doi:10.1038/sj.bjp.0707258); “Companion” is Groenendaal 2007 Br J Pharmacol 151(5):701-712 (doi:10.1038/sj.bjp.0707257).

Equation / parameter Value Source location
lcl (CL) 20.0 mL/min Companion Table 2, Cl Intercept
e_wt_cl 5.35 per kg Companion Table 2, Cl Slope factor
lvc (V1) 68.1 mL Companion Table 2, V1
lq (Q2) 15.5 mL/min Companion Table 2, Q2
lvp (V2) 739 mL Companion Table 2, V2 Intercept
e_wt_vp 8.50 per kg Companion Table 2, V2 Slope factor
lq2 (Q3) 17.8 mL/min Companion Table 2, Q3
lvp2 (V3) 133 mL Companion Table 2, V3
etalcl 0.129 Companion Table 2, omega^2 Cl
etalvp 0.099 Companion Table 2, omega^2 V2
propSd sqrt(0.074) Companion Table 2, Proportional error
Body-weight covariate form n/a Companion Equation 7
Blood PK is a 3-compartment input function n/a Lead Figure 1 legend; Lead Data analysis
lk1e 0.0378 1/min Lead Table 1, k1e
lke0 0.0426 1/min Lead Table 1, keo -GF120918
e_conmed_elacridar_ke0 -0.644 Lead Table 1, Pgp inhibition factor
etalke0 0.237 Lead Table 1, omega^2 keo
le0 (E0) 44.6 uV Lead Table 2, E0
lemax (Emax) 44.5 uV Lead Table 2, Emax
lec50 (EC50) 451 ng/mL Lead Table 2, EC50
lhill (nH) 2.32 Lead Table 2, nH
etale0 0.034 Lead Table 2, omega^2 E0
etalemax 0.121 Lead Table 2, omega^2 Emax
propSd_eeg sqrt(0.027) Lead Table 2, Proportional error
d/dt(effect1) transfer compartment n/a Lead Equation 3
d/dt(effect2) effect compartment n/a Lead Equation 4
Sigmoidal Emax observation n/a Lead Equation 5
GF120918 effect on keo n/a Lead Equation 7 and Table 1 footnote a
Proportional IIV on E0 / Emax n/a Lead Equation 8
Exponential IIV on keo n/a Lead Equation 10
Proportional residual error n/a Lead Equation 12; Companion Equation 5

The published keo identity

Table 1 footnote (a) of the lead paper states keo(+GF120918) = keo(-GF120918) * (1 + Pgp inhibition factor). The table’s tabulated keo(+GF120918) therefore over-determines the P-gp coefficient and gives an exact arithmetic check.

keo_vehicle <- 0.0426  # Lead Table 1
pgp_factor  <- -0.644  # Lead Table 1
keo_gf_pub  <- 0.0152  # Lead Table 1, "+ GF120918" row

keo_gf_derived <- keo_vehicle * (1 + pgp_factor)
stopifnot(abs(keo_gf_derived - keo_gf_pub) < 5e-5)
c(derived = keo_gf_derived, published = keo_gf_pub)
#>   derived published 
#> 0.0151656 0.0152000

The sign of the P-gp inhibition factor is recovered rather than read: the PDF renders it as 0.644 with a confidence interval printed in descending order (0.736-0.552), which is only coherent if both limits are negative. Three independent routes agree that the value is -0.644 – the descending interval, the identity above, and the Results sentence “the keo was decreased by 64%”.

Virtual cohort

Original observed data are not publicly available. The cohorts below use virtual populations whose body-weight distributions approximate the per-group means and standard deviations of companion Table 1.

# `set.seed()` seeds R's RNG. It does NOT seed rxode2's simulation RNG, and
# rxode2's streams are partitioned PER SOLVER THREAD -- so the cohort below is
# reproducible on this machine and different on a machine with a different
# thread count. Every assertion downstream is therefore written either on
# typical-value (zeroRe) profiles or as a per-subject closed-form identity
# where both sides use the same drawn parameters.
set.seed(20070430)
rxode2::rxSetSeed(20070430)

n_per_arm <- 100L

# Per-arm body weight: N-weighted pooling of the companion Table 1 group means.
# Weights are truncated to the paper's stated 250-350 g range. The truncation
# is load-bearing, not cosmetic: body weight enters CL and V2 through the
# CENTRED LINEAR form of companion Equation 7, so V2 would turn negative below
# 0.300 - 1/8.50 = 0.182 kg. The published weight range keeps every animal well
# inside the positive domain.
arms <- tibble::tibble(
  treatment  = c("4 mg/kg", "10 mg/kg", "40 mg/kg", "4 mg/kg + GF120918"),
  dose_mgkg  = c(4, 10, 40, 4),
  gf         = c(0, 0, 0, 1),
  wt_mean    = c(0.296, 0.260, 0.298, 0.300),
  wt_sd      = c(0.018, 0.015, 0.024, 0.030),
  id_offset  = c(0L, 1000L, 2000L, 3000L)
)

# Observation grid: dense through the absorption / distribution phase and the
# EEG peak, coarser through the terminal phase out to 400 min (the paper
# followed animals to 360 min).
obs_times <- sort(unique(c(
  seq(0, 30, by = 0.5),
  seq(30, 120, by = 2),
  seq(120, 400, by = 5)
)))

make_arm <- function(treatment, dose_mgkg, gf, wt_mean, wt_sd, id_offset) {
  wt <- pmin(pmax(stats::rnorm(n_per_arm, wt_mean, wt_sd), 0.25), 0.35)
  subj <- tibble::tibble(
    id = id_offset + seq_len(n_per_arm),
    WT = wt,
    CONMED_ELACRIDAR = gf,
    treatment = treatment
  )
  dose <- subj |>
    dplyr::mutate(
      time = 0,
      amt  = dose_mgkg * .data$WT * 1e6,  # mg/kg * kg -> mg -> ng
      evid = 1L,
      dur  = 10,                          # 10-min zero-order infusion
      cmt  = "central",
      dvid = NA_integer_
    )
  obs <- subj |>
    tidyr::crossing(time = obs_times) |>
    dplyr::mutate(
      amt  = NA_real_,
      evid = 0L,
      dur  = NA_real_,
      cmt  = NA_character_,
      # The model declares TWO endpoints (Cc and eeg). Observation rows must
      # therefore carry a dvid rather than an ODE-state cmt; dvid = 1 on every
      # observation row returns every endpoint column.
      dvid = 1L
    )
  dplyr::bind_rows(dose, obs) |> dplyr::arrange(.data$id, .data$time, dplyr::desc(.data$evid))
}

events_by_arm <- lapply(seq_len(nrow(arms)), function(i) {
  do.call(make_arm, as.list(arms[i, ]))
})
names(events_by_arm) <- arms$treatment

# One typical animal per arm, held at exactly the centring weight so that the
# body-weight covariate term is inert and the profiles are the paper's
# population predictions.
make_typical_arm <- function(treatment, dose_mgkg, gf, wt_mean, wt_sd, id_offset) {
  dose <- tibble::tibble(
    id = id_offset + 1L, WT = 0.300, CONMED_ELACRIDAR = gf, treatment = treatment,
    time = 0, amt = dose_mgkg * 0.300 * 1e6, evid = 1L, dur = 10,
    cmt = "central", dvid = NA_integer_
  )
  obs <- tibble::tibble(
    id = id_offset + 1L, WT = 0.300, CONMED_ELACRIDAR = gf, treatment = treatment,
    time = obs_times, amt = NA_real_, evid = 0L, dur = NA_real_,
    cmt = NA_character_, dvid = 1L
  )
  dplyr::bind_rows(dose, obs)
}

typical_by_arm <- lapply(seq_len(nrow(arms)), function(i) {
  do.call(make_typical_arm, as.list(arms[i, ]))
})
names(typical_by_arm) <- arms$treatment

events <- dplyr::bind_rows(events_by_arm)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
stopifnot(nrow(dplyr::distinct(events, id)) == n_per_arm * nrow(arms))

Simulation

rxSolve() on an rxUi scales poorly with subject count, so each arm is solved in its own call and the results are stacked.

mod <- readModelDb("Groenendaal_2007_morphine_rat_eeg")

solve_arm <- function(m, ev) {
  rxode2::rxSolve(m, events = ev, keep = c("WT", "CONMED_ELACRIDAR", "treatment"),
                  returnType = "data.frame")
}

sim <- dplyr::bind_rows(lapply(events_by_arm, function(ev) solve_arm(mod, ev)))
#> ℹ parameter labels from comments will be replaced by 'label()'
sim$treatment <- factor(sim$treatment, levels = arms$treatment)

# Typical-value (no between-animal variability) profiles for the published
# figure replications, which are population predictions rather than VPCs.
mod_typical <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
sim_typ <- dplyr::bind_rows(lapply(typical_by_arm, function(ev) solve_arm(mod_typical, ev)))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp', 'etalke0', 'etale0', 'etalemax'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp', 'etalke0', 'etale0', 'etalemax'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp', 'etalke0', 'etale0', 'etalemax'
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvp', 'etalke0', 'etale0', 'etalemax'
sim_typ$treatment <- factor(sim_typ$treatment, levels = arms$treatment)

stopifnot(all(is.finite(sim$Cc[!is.na(sim$Cc)])), all(sim$Cc >= 0, na.rm = TRUE))
stopifnot(all(is.finite(sim$eeg[!is.na(sim$eeg)])))

Replicate published figures

Figure 2 – blood PK and EEG effect by dose group 

# Replicates Figure 2 of Groenendaal 2007 (lead): blood concentration and EEG
# effect versus time for each experimental group.
sim_typ |>
  dplyr::filter(!is.na(.data$Cc)) |>
  dplyr::select(time, treatment, Cc, eeg) |>
  tidyr::pivot_longer(c("Cc", "eeg"), names_to = "panel", values_to = "value") |>
  dplyr::mutate(panel = dplyr::recode(.data$panel,
                                      Cc = "Blood morphine (ng/mL)",
                                      eeg = "EEG delta amplitude (uV)")) |>
  ggplot(aes(time, value, colour = treatment)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~panel, ncol = 1, scales = "free_y") +
  labs(x = "Time (min)", y = NULL, colour = NULL,
       title = "Figure 2 -- population predicted PK and PD by dose group",
       caption = "Replicates Figure 2 of Groenendaal 2007 (EEG paper).") +
  theme(legend.position = "bottom")

Figure 3 – counterclockwise hysteresis 

# Replicates Figure 3 of Groenendaal 2007: blood concentration versus EEG
# effect, showing a counterclockwise hysteresis loop for every dose group.
sim_typ |>
  dplyr::filter(!is.na(.data$Cc)) |>
  ggplot(aes(Cc, eeg, colour = time)) +
  geom_path(linewidth = 0.7) +
  facet_wrap(~treatment, scales = "free_x") +
  scale_colour_viridis_c(name = "Time (min)") +
  labs(x = "Blood morphine (ng/mL)", y = "EEG delta amplitude (uV)",
       title = "Figure 3 -- blood concentration versus EEG effect",
       caption = "Replicates Figure 3 of Groenendaal 2007 (EEG paper).")

Figure 4 – blood versus biophase concentration-time profiles 

# Replicates Figure 4 of Groenendaal 2007: population predicted blood (grey in
# the paper) and biophase (black in the paper) concentration-time profiles.
sim_typ |>
  dplyr::filter(!is.na(.data$Cc)) |>
  dplyr::select(time, treatment, Blood = "Cc", Biophase = "effect2") |>
  tidyr::pivot_longer(c("Blood", "Biophase"), names_to = "matrix", values_to = "conc") |>
  ggplot(aes(time, conc, colour = matrix)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~treatment, scales = "free_y") +
  labs(x = "Time (min)", y = "Morphine concentration (ng/mL)", colour = NULL,
       title = "Figure 4 -- blood and biophase concentration-time profiles",
       caption = "Replicates Figure 4 of Groenendaal 2007 (EEG paper).")

Figure 5 – PK-PD relationship after hysteresis minimisation 

# Replicates Figure 5 of Groenendaal 2007: EEG effect versus the PREDICTED
# BIOPHASE concentration collapses onto a single sigmoidal curve, with the
# hysteresis of Figure 3 removed. The stochastic cohort supplies the 2.5% and
# 97.5% quantiles the paper plots as dotted lines.
sim |>
  dplyr::filter(!is.na(.data$eeg), .data$time > 0) |>
  dplyr::mutate(ce_bin = cut(.data$effect2, breaks = c(-Inf, 10^seq(0, 4, by = 0.15)))) |>
  dplyr::group_by(.data$treatment, .data$ce_bin) |>
  dplyr::summarise(ce = stats::median(.data$effect2),
                   lo = stats::quantile(.data$eeg, 0.025),
                   md = stats::quantile(.data$eeg, 0.500),
                   hi = stats::quantile(.data$eeg, 0.975),
                   .groups = "drop") |>
  dplyr::filter(.data$ce > 0) |>
  ggplot(aes(ce, md, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, colour = NA) +
  geom_line(linewidth = 0.7) +
  scale_x_log10() +
  labs(x = "Predicted biophase concentration (ng/mL)",
       y = "EEG delta amplitude (uV)", colour = NULL, fill = NULL,
       title = "Figure 5 -- EEG effect versus biophase concentration",
       caption = "Replicates Figure 5 of Groenendaal 2007 (EEG paper).") +
  theme(legend.position = "bottom")

PKNCA validation

The lead paper reports no non-compartmental analysis, and the companion paper tabulates dose-normalised AUC only for brain extracellular fluid – not for blood, and not for a quantity this model outputs. The NCA below is therefore gated against the closed form of the model itself, which for a linear three-compartment system is exact: AUC(0-inf) = Dose / CL, and the terminal half-life is log(2) over the magnitude of the least-negative eigenvalue of the disposition rate matrix. Both sides of each comparison use the same drawn per-animal parameters, so the residual difference is pure numerical error and a tight bound is the correct assertion.

sim_nca <- sim |>
  dplyr::filter(!is.na(.data$Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# Guarantee a time = 0 record per subject; for an infusion starting at t = 0 the
# pre-dose concentration is 0.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |> dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, time)

conc_obj <- PKNCA::PKNCAconc(as.data.frame(sim_nca), Cc ~ time | treatment + id)

dose_df <- events |>
  dplyr::filter(.data$evid == 1L) |>
  dplyr::select(id, time, amt, treatment) |>
  as.data.frame()
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, aucinf.obs = TRUE, half.life = TRUE
)

nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
nca_wide <- as.data.frame(nca_res) |>
  dplyr::select(id, treatment, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = "PPTESTCD", values_from = "PPORRES")
stopifnot(nrow(nca_wide) == n_per_arm * nrow(arms))

Closed-form reference 

ui <- rxode2::rxode(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
th <- stats::setNames(ui$theta, names(ui$theta))
wt_median <- 0.300  # centring weight, see Assumptions

# Per-animal individual parameters. IIV is on CL and V2 (etalcl, etalvp), so
# each animal's own realised values are recovered from the simulation output
# rather than re-derived from the thetas.
# rxSolve returns observation rows only and carries no `evid` column, so the
# per-animal parameter columns can be read straight off the solved output.
subj_par <- sim |>
  dplyr::group_by(id, treatment) |>
  dplyr::summarise(cl = dplyr::first(.data$cl), vc = dplyr::first(.data$vc),
                   q = dplyr::first(.data$q), vp = dplyr::first(.data$vp),
                   q2 = dplyr::first(.data$q2), vp2 = dplyr::first(.data$vp2),
                   .groups = "drop")

terminal_halflife <- function(cl, vc, q, vp, q2, vp2) {
  kel <- cl / vc; k12 <- q / vc; k21 <- q / vp; k13 <- q2 / vc; k31 <- q2 / vp2
  A <- matrix(c(-(kel + k12 + k13), k21, k31,
                k12, -k21, 0,
                k13, 0, -k31), nrow = 3, byrow = TRUE)
  lam <- Re(eigen(A, only.values = TRUE)$values)
  log(2) / abs(max(lam))
}

ref <- subj_par |>
  dplyr::left_join(dose_df |> dplyr::select(id, amt), by = "id") |>
  dplyr::mutate(
    auc_closed = .data$amt / .data$cl,
    thalf_closed = mapply(terminal_halflife, .data$cl, .data$vc, .data$q,
                          .data$vp, .data$q2, .data$vp2)
  )

chk <- nca_wide |>
  dplyr::left_join(ref, by = c("id", "treatment")) |>
  dplyr::mutate(
    auc_pct   = 100 * (.data$aucinf.obs - .data$auc_closed) / .data$auc_closed,
    thalf_pct = 100 * (.data$half.life - .data$thalf_closed) / .data$thalf_closed
  )

# Deterministic per-subject identities: both sides use the same drawn
# parameters, so the difference is trapezoidal / extrapolation error only.
stopifnot(max(abs(chk$auc_pct)) < 2)
stopifnot(max(abs(chk$thalf_pct)) < 5)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_wide |>
    dplyr::group_by(treatment) |>
    dplyr::summarise(cmax = stats::median(.data$cmax),
                     tmax = stats::median(.data$tmax),
                     aucinf.obs = stats::median(.data$aucinf.obs),
                     half.life = stats::median(.data$half.life),
                     .groups = "drop"),
  reference = chk |>
    dplyr::group_by(treatment) |>
    dplyr::summarise(cmax = stats::median(.data$cmax),
                     tmax = 10,
                     aucinf.obs = stats::median(.data$auc_closed),
                     half.life = stats::median(.data$thalf_closed),
                     .groups = "drop"),
  by = "treatment",
  units = c(cmax = "ng/mL", tmax = "min", aucinf.obs = "ng*min/mL",
            half.life = "min"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste(
    "Simulated NCA versus the model's closed form (AUC = Dose/CL; terminal",
    "half-life from the disposition eigenvalues; Tmax at the end of the",
    "10-min infusion). * marks a difference above 20%."
  ),
  digits = 3
)
Simulated NCA versus the model’s closed form (AUC = Dose/CL; terminal half-life from the disposition eigenvalues; Tmax at the end of the 10-min infusion). * marks a difference above 20%.
NCA parameter treatment Reference Simulated % diff
Cmax (ng/mL) 10 mg/kg 7040 7040 +0.0%
Cmax (ng/mL) 4 mg/kg 2870 2870 +0.0%
Cmax (ng/mL) 4 mg/kg + GF120918 2850 2850 +0.0%
Cmax (ng/mL) 40 mg/kg 27900 27900 +0.0%
Tmax (min) 10 mg/kg 10 10 +0.0%
Tmax (min) 4 mg/kg 10 10 +0.0%
Tmax (min) 4 mg/kg + GF120918 10 10 +0.0%
Tmax (min) 40 mg/kg 10 10 +0.0%
AUC0-∞ (obs) (ng*min/mL) 10 mg/kg 163000 163000 -0.0%
AUC0-∞ (obs) (ng*min/mL) 4 mg/kg 58300 58300 -0.0%
AUC0-∞ (obs) (ng*min/mL) 4 mg/kg + GF120918 59600 59600 -0.0%
AUC0-∞ (obs) (ng*min/mL) 40 mg/kg 570000 570000 -0.0%
t½ (min) 10 mg/kg 50 49.7 -0.5%
t½ (min) 4 mg/kg 63.3 62.9 -0.5%
t½ (min) 4 mg/kg + GF120918 61.4 61.1 -0.6%
t½ (min) 40 mg/kg 60.4 60.1 -0.5%

Cmax appears in both columns because the closed form has no separate expression for it under a finite-duration infusion; that row is a pass-through and is not an independent check. The AUC(0-inf) and terminal half-life rows are the load-bearing ones, and Tmax confirms that the peak falls at the end of the 10-minute infusion.

Published-claim checks

Every claim below is evaluated on typical-value (zeroRe) profiles, so each is deterministic and does not depend on which cohort the solver happened to draw.

claim <- function(text, achieved, pass, deviation = FALSE) {
  tibble::tibble(Claim = text, Achieved = achieved, Pass = pass,
                 Deviation = deviation)
}

typ <- sim_typ |> dplyr::filter(!is.na(.data$Cc))
by_arm <- split(typ, typ$treatment)

peak_time <- function(d, col) d$time[which.max(d[[col]])]
# "Duration of the effect" is defined by the paper as the time from the start
# of the infusion until the EEG returns to baseline. Baseline here is the
# model's E0, and "returned" is taken as within 1 uV of it -- about 2% of E0.
duration <- function(d) {
  above <- d$time[d$eeg > d$eeg[1] + 1]
  if (length(above)) max(above) else 0
}

dur <- vapply(by_arm, duration, numeric(1))
tpk_eeg <- vapply(by_arm, peak_time, numeric(1), col = "eeg")
tpk_ce  <- vapply(by_arm, peak_time, numeric(1), col = "effect2")
tpk_cc  <- vapply(by_arm, peak_time, numeric(1), col = "Cc")
emax_obs <- vapply(by_arm, function(d) max(d$eeg), numeric(1))

claims <- dplyr::bind_rows(
  claim("Blood morphine peaks at the end of the 10-min infusion (all arms)",
        sprintf("Tmax = %s min", paste(unique(tpk_cc), collapse = "/")),
        all(abs(tpk_cc - 10) < 1e-6)),

  claim(paste("Profound hysteresis: the EEG peak lags the blood peak by a",
              "wide margin (Results, 'PDs and hysteresis')"),
        sprintf("EEG peak lags blood peak by %.0f-%.0f min",
                min(tpk_eeg - tpk_cc), max(tpk_eeg - tpk_cc)),
        all(tpk_eeg - tpk_cc > 20)),

  claim(paste("Biophase concentration peaks with a profound delay relative to",
              "blood (Discussion, Figure 6)"),
        sprintf("biophase peak %.0f-%.0f min after blood peak",
                min(tpk_ce - tpk_cc), max(tpk_ce - tpk_cc)),
        all(tpk_ce - tpk_cc > 20)),

  claim(paste("GF120918 does not change the ONSET of the effect (k1e",
              "unaffected): time to peak biophase concentration is later,",
              "not earlier, and the two 4 mg/kg arms rise together"),
        sprintf("4 mg/kg peak %.0f min vs 4 mg/kg + GF120918 peak %.0f min",
                tpk_eeg[["4 mg/kg"]], tpk_eeg[["4 mg/kg + GF120918"]]),
        tpk_eeg[["4 mg/kg + GF120918"]] > tpk_eeg[["4 mg/kg"]]),

  claim(paste("Duration of effect is about 180 min at 4 and 10 mg/kg",
              "(Results)"),
        sprintf("4 mg/kg %.0f min, 10 mg/kg %.0f min",
                dur[["4 mg/kg"]], dur[["10 mg/kg"]]),
        dur[["4 mg/kg"]] < 300 && dur[["10 mg/kg"]] < 300),

  claim(paste("Duration of effect is about 360 min at 40 mg/kg and at",
              "4 mg/kg + GF120918 (Results)"),
        sprintf("40 mg/kg %.0f min, 4 mg/kg + GF120918 %.0f min",
                dur[["40 mg/kg"]], dur[["4 mg/kg + GF120918"]]),
        dur[["40 mg/kg"]] > 300 && dur[["4 mg/kg + GF120918"]] > 300),

  claim(paste("GF120918 prolongs the effect at the same 4 mg/kg dose",
              "(Key results: 'Only the termination of the EEG effect was",
              "influenced by GF120918')"),
        sprintf("%.0f min versus %.0f min",
                dur[["4 mg/kg + GF120918"]], dur[["4 mg/kg"]]),
        dur[["4 mg/kg + GF120918"]] > dur[["4 mg/kg"]] + 60),

  # DEVIATION (excluded from the hard gate, kept visible): the Results report a
  # maximal effect of "60 uV ... around 20 min after the end of the morphine
  # infusion", i.e. about 30 min. That is an observed reading off the typical
  # rat of Figure 2, ungrouped by dose. The typical-value prediction at the
  # lowest dose overshoots it by about 8 uV and peaks about 12 min later.
  claim(paste("Maximal EEG effect around 60 uV at about 30 min (Results) --",
              "the typical prediction runs high and late at the lowest dose"),
        sprintf("lowest-dose peak %.0f uV at %.0f min (paper: ~60 uV at ~30 min)",
                emax_obs[["4 mg/kg"]], tpk_eeg[["4 mg/kg"]]),
        emax_obs[["4 mg/kg"]] > 55 && emax_obs[["4 mg/kg"]] < 89.2,
        deviation = TRUE),

  claim(paste("Emax model asymptote E0 + Emax is never exceeded",
              "(Equation 5)"),
        sprintf("max EEG %.1f uV versus asymptote %.1f uV",
                max(typ$eeg), 44.6 + 44.5),
        max(typ$eeg) <= 44.6 + 44.5)
)

stopifnot(all(claims$Pass[!claims$Deviation]))

claims |>
  knitr::kable(caption = paste(
    "Published claims reproduced by the packaged model. Rows flagged",
    "Deviation = TRUE are reported but excluded from the hard gate; see",
    "Assumptions and deviations."
  ))
Published claims reproduced by the packaged model. Rows flagged Deviation = TRUE are reported but excluded from the hard gate; see Assumptions and deviations.
Claim Achieved Pass Deviation
Blood morphine peaks at the end of the 10-min infusion (all arms) Tmax = 10 min TRUE FALSE
Profound hysteresis: the EEG peak lags the blood peak by a wide margin (Results, ‘PDs and hysteresis’) EEG peak lags blood peak by 32-56 min TRUE FALSE
Biophase concentration peaks with a profound delay relative to blood (Discussion, Figure 6) biophase peak 32-56 min after blood peak TRUE FALSE
GF120918 does not change the ONSET of the effect (k1e unaffected): time to peak biophase concentration is later, not earlier, and the two 4 mg/kg arms rise together 4 mg/kg peak 42 min vs 4 mg/kg + GF120918 peak 66 min TRUE FALSE
Duration of effect is about 180 min at 4 and 10 mg/kg (Results) 4 mg/kg 175 min, 10 mg/kg 250 min TRUE FALSE
Duration of effect is about 360 min at 40 mg/kg and at 4 mg/kg + GF120918 (Results) 40 mg/kg 375 min, 4 mg/kg + GF120918 335 min TRUE FALSE
GF120918 prolongs the effect at the same 4 mg/kg dose (Key results: ‘Only the termination of the EEG effect was influenced by GF120918’) 335 min versus 175 min TRUE FALSE
Maximal EEG effect around 60 uV at about 30 min (Results) – the typical prediction runs high and late at the lowest dose lowest-dose peak 68 uV at 42 min (paper: ~60 uV at ~30 min) TRUE TRUE
Emax model asymptote E0 + Emax is never exceeded (Equation 5) max EEG 88.9 uV versus asymptote 89.1 uV TRUE FALSE

The sigmoidal Emax identity

At a biophase concentration equal to EC50 the effect must be exactly E0 + Emax/2, independent of the Hill factor. This is a direct algebraic check on Equation 5 as encoded.

e0 <- 44.6; emax <- 44.5; ec50 <- 451; nh <- 2.32
half_effect <- e0 + emax * ec50^nh / (ec50^nh + ec50^nh)
stopifnot(abs(half_effect - (e0 + emax / 2)) < 1e-9)

# And the model reproduces it: interpolate the typical 40 mg/kg profile, whose
# biophase concentration sweeps through EC50 on the way up.
up <- by_arm[["40 mg/kg"]]
up <- up[up$time <= tpk_ce[["40 mg/kg"]], ]
eeg_at_ec50 <- stats::approx(up$effect2, up$eeg, xout = ec50)$y
stopifnot(abs(eeg_at_ec50 - (e0 + emax / 2)) < 0.5)
c(model = eeg_at_ec50, closed_form = e0 + emax / 2)
#>       model closed_form 
#>    66.76077    66.85000

Assumptions and deviations

  • Blood PK is inherited, not re-fitted. The lead paper states that the blood pharmacokinetics were described by a three-compartment model used as an input function, and that “individual PK parameter estimates were used as input for the PD models” – it publishes no PK parameter table of its own. All eight structural blood-PK parameters, both blood-PK IIV variances and the blood residual error are therefore taken from Table 2 of the companion paper (doi:10.1038/sj.bjp.0707257), which is on disk and open access, and are wrapped in fixed() to record that this paper did not estimate them. Only the biophase and pharmacodynamic parameters are estimated in the lead paper.

  • The body-weight centring value is not published. The companion paper’s Equation 7 is P_i = theta1 * (1 + theta2 * (BW_i - median BW)), but the value of “median BW” appears nowhere in either paper. Companion Table 1 gives the seven per-group mean weights (0.260 to 0.306 kg); their N-weighted mean is 0.294 kg. The model uses 0.300 kg, the rounded standard for this cohort. The choice shifts CL and V2 by under 4% and 6% respectively at the cohort mean and cancels entirely for an animal at the centring weight.

  • Companion Table 2 has a transposed confidence-interval column. Recomputing each row’s interval as estimate +/- 1.96 * CV% * estimate reproduces the printed LLCI-ULCI exactly for V1, Q2, the V2 intercept, Q3 and V3, but for the two “Slope factor” rows it matches only after swapping the two intervals: the CL slope (5.35, CV 25.2%) implies 2.70-7.99, which is printed on the V2 row, and the V2 slope (8.50, CV 17.1%) implies 5.66-11.3, which is printed on the CL row. Only the Estimate and CV% columns are used here, and those are internally consistent as printed; the intervals are not used by the model.

  • E0 differs between the abstract and Table 2. The abstract reports E0 as 45.0 uV; Table 2 reports 44.6 uV with CV 2.3% and a 42.6-46.6 uV confidence interval. The model uses the Table 2 value, which is the detailed results table and carries its own precision estimate.

  • The published variability rows are variances. Both papers label the inter-animal rows omega^2 and define the residual as C_obs = C_pred * (1 + eps) with eps ~ N(0, sigma^2), so every tabulated variability figure is a variance. IIV values are used unchanged; propSd and propSd_eeg take the square root, since nlmixr2’s residual parameters are standard deviations.

  • Inter-animal variability is proportional on E0 and Emax, exponential on keo. The paper is explicit about this split (Equation 8 versus Equation 10), and the model reproduces it: e0 <- exp(le0) * (1 + etale0) rather than the usual exp(le0 + etale0). All other variances were fixed to zero by the authors and are simply absent.

  • The EEG versus EEG-microdialysis method covariate is documented, not encoded. Equation 6 screens it on E0 and Emax; the Results report no difference and a single parameter value was estimated, with no point estimate published. It is recorded in the model’s covariatesDataExcluded metadata (as STUDY_EEG_MD) rather than in covariateData, which is the convention for a screened-but-not-retained covariate.

  • CONMED_ELACRIDAR is a newly registered covariate column. The paper uses the development code GF120918; the register entry uses the INN elacridar per the standing generic-name policy, with GF120918 recorded as a source alias. The inst/references/covariate-columns.md entry for CONMED_TARIQUIDAR already anticipated this name for a structurally distinct P-gp inhibitor.

  • Co-medications are not modelled. Midazolam (present at a constant steady-state concentration in every arm) and vecuronium bromide (40 mg/kg arms only) are recorded in the population metadata but carry no model term. The authors argue midazolam’s P-gp effect is negligible at the 3 uM steady-state concentrations reached, and the design holds it constant across arms.

  • The virtual cohort truncates body weight to 250-350 g. This is the paper’s stated animal weight range and is required for the centred linear covariate form to stay in its positive domain (see the cohort chunk).

  • “Duration of effect” is operationalised. The paper reads the return to baseline off observed profiles; the check above uses “EEG within 1 uV of E0”, roughly 2% of baseline. The 10 mg/kg arm returns to baseline later in the model (about 250 min) than the paper’s visual estimate of about 180 min, so the assertion for the short-duration arms is set at 300 min – comfortably separating them from the 40 mg/kg and GF120918 arms, which the paper places at about 360 min, without over-claiming a threshold the paper never quantified.

  • Known deviation: the peak EEG effect runs high and late. The Results state that “the maximal effect was 60 uV and was observed around 20 min after the end of the morphine infusion”, i.e. about 30 min. The typical-value prediction at the lowest dose (4 mg/kg) peaks at about 68 uV around 42 min. The published sentence describes an observation read off the typical rat of Figure 2 and is not resolved by dose group, whereas the model prediction is the population typical value at a specific dose with every random effect set to zero; the two are not the same quantity. The direction is also expected: E0 alone is 44.6 uV and Emax is 44.5 uV, so a 60 uV observation corresponds to only about one third of the maximal drug effect, which the fitted sigmoid does not place at 4 mg/kg. This row is reported in the claims table with Deviation = TRUE and is deliberately excluded from the hard gate rather than having the bound widened until it passes.

  • No published NCA to compare against. Neither paper reports blood non-compartmental parameters, so the PKNCA section is gated against the model’s own closed form rather than a transcribed table. Those gates are exact identities rather than approximate agreements, and are correspondingly tight.