Opioid overdose reversal simulation chain (Laffont 2025)
Source:vignettes/articles/Laffont_2025_opioid_overdose_reversal_simulation.Rmd
Laffont_2025_opioid_overdose_reversal_simulation.RmdWhat this vignette does
This vignette demonstrates the end-to-end synthetic-opioid overdose simulation chain that Laffont, Purohit, de la Pena and Skolnick (2025, Neuropharmacology 278:110546, doi:10.1016/j.neuropharm.2025.110546, PMID 40484316) used to compare the effectiveness of intranasal (IN) naloxone (4 mg) and IN nalmefene (3 mg HCl, 2.7 mg base) for reversing fentanyl and carfentanil overdoses. The Laffont 2025 study is a simulation study; it contributes no new equations or parameter estimates. Every structural piece comes from one of two upstream publications already in nlmixr2lib:
The opioid PK + mu-receptor binding + respiratory-physiology + cardiac-arrest chain was developed by Mann et al. (Clin Pharmacol Ther 2022;112(5):1020-1032, doi:10.1002/cpt.2696, PMID 35766413) and published as the FDA Division of Applied Regulatory Science translational opioid-overdose model (https://github.com/FDA/Mechanistic-PK-PD-Model-to-Rescue-Opioid-Overdose). Mann 2022 is packaged here as
Mann_2022_fentanyl_iv,Mann_2022_carfentanil_iv,Mann_2022_mu_receptor_binding, andMann_2022_respiratory_physiology.The IN naloxone and IN nalmefene popPK models, plus the nalmefene binding parameters (added as receptor-binding ligand 13), come from Laffont CM, Purohit P, Delcamp N, Gonzalez-Garcia I, Skolnick P. Comparison of intranasal naloxone and intranasal nalmefene in a translational model assessing the impact of synthetic opioid overdose on respiratory depression and cardiac arrest. Front Psychiatry 2024;15:1399803, doi:10.3389/fpsyt.2024.1399803. Laffont 2024 is packaged as
Laffont_2024_naloxoneandLaffont_2024_nalmefene.
The vignette’s job is to put these layers together, run a coarse simulation grid (small N for the pkgdown 5-minute time budget), and verify that the chain reproduces the headline cardiac-arrest / brain-hypoxia numbers Laffont 2025 reports in Table 1 and Figs 2-7. It is NOT a 2000-subject reproduction of the published study.
Layer-to-paper mapping
| Layer | Source publication | nlmixr2lib model |
|---|---|---|
| IV fentanyl popPK (3-compartment + biophase) | Mann 2022 (Table S1) | Mann_2022_fentanyl_iv |
| IV carfentanil PK (rate-constant-modified fentanyl) | Mann 2022 (Methods) | Mann_2022_carfentanil_iv |
| Mu-opioid receptor competitive binding (13 ligands) | Mann 2022 (Table S2, 12 ligands) + Laffont 2024 (Supp Table S3, nalmefene = ligand 13) | Mann_2022_mu_receptor_binding |
| Respiratory + cerebrovascular physiology + cardiac-arrest | Mann 2022 (Supplement 1) | Mann_2022_respiratory_physiology |
| IN naloxone popPK (2-compartment + parallel zero/first-order) | Laffont 2024 (Table 2) | Laffont_2024_naloxone |
| IN nalmefene popPK (2-compartment + parallel zero/first-order) | Laffont 2024 (Table 1) | Laffont_2024_nalmefene |
| Antagonist plasma -> brain ke0 = 0.001774 1/s | Mann 2022 (Supp), shared with Laffont 2024 Supp Table S3 | applied in this vignette during PK-to-Ce post-processing |
| Antagonist Ce (ng/mL) -> Ce (pM) conversion | Free-base MW: naloxone 327.37 g/mol; nalmefene 339.43 g/mol | applied in this vignette |
Laffont 2025 itself contributes no new equations, no new parameter estimates, no new structural model. Its contribution is the simulation conditions (additional lower opioid doses, additional intervention-delay scenarios) and the reporting endpoints (cardiac arrest incidence and brain-hypoxia duration tabulated in Table 1 and Figs 2-7). This vignette reproduces a subset of those scenarios.
Population
Laffont 2025 simulates 2000 virtual subjects per arm (with binding parameter sets drawn from Mann 2022’s bootstrap distribution and PK parameter sets drawn from the Laffont 2024 nonlinear mixed-effects models). Two subject types are used: opioid-naive healthy volunteers (P1 = 2.875, P3 = 0.9 in the Mann respiratory PD layer) and chronic opioid users (P1 = 4.226, P3 = 1.323). Body weight is fixed at 70 kg for the opioid PK (Mann 2022 convention) and 74.7 kg for the nalmefene PK (Laffont 2024 allometric reference).
The grid below uses N = 30 virtual chronic opioid users (small N for the 5-minute pkgdown time budget; the headline Laffont 2025 N is 2000). The cardiac-arrest endpoint produced by this small-N chain will land within sampling-noise distance of the Laffont 2025 / Mann 2022 published values without ever being tuned to match.
Source trace
Per-parameter origin is recorded as in-file comments in each model
file. The table below collects the source pointers for the values this
vignette uses from non-paper code (i.e., values that are computed in the
vignette rather than read from a model file’s ini()).
| Value | Magnitude | Source |
|---|---|---|
| Antagonist effect-site ke0 | 1.774e-3 1/s = 0.1064 1/min | Mann 2022 (carried into Laffont 2024 Supp Table S3 for both nalmefene and naloxone) |
| Naloxone free-base molecular weight | 327.37 g/mol | PubChem CID 5284596 (naloxone base) |
| Nalmefene free-base molecular weight | 339.43 g/mol | PubChem CID 5284594 (nalmefene base) |
| Fentanyl IV bolus doses | 0.82 / 1.63 / 2.97 mg | Laffont 2025 Methods + Table 1 |
| Carfentanil IV bolus doses | 0.010 / 0.012 / 0.022 mg | Laffont 2025 Methods + Table 1 |
| IN naloxone HCl rescue dose | 4 mg | Laffont 2025 Methods (FDA-approved fixed dose) |
| IN nalmefene HCl rescue dose | 3 mg (= 2.7 mg base) | Laffont 2025 Methods (FDA-approved fixed dose) |
| Antagonist administration delay | 1 min after Venti hits 40% baseline | Mann 2022 (original convention; Laffont 2025 also varied 0.5-10 min) |
| Cardiac-arrest definition | total Q < 0.01 L/min | Mann 2022 (the im_arrest collapse-decay multiplier on
cardiac output) |
| Brain hypoxia threshold | PbtO2 < 20 mm Hg | Laffont 2025 (citing Lund 2024); the model emits brain O2 partial
pressure P_B_O2
|
Composed chain: opioid PK -> binding -> physiology -> cardiac arrest
# We build a SINGLE rxode2 model that integrates the opioid PK (3 cmt +
# biophase), competitive binding (RL_op, RL_antag), and respiratory /
# cardiovascular physiology / cardiac-arrest dynamics from Mann 2022.
# Parameter values are the typical-value estimates from Mann 2022
# Tables S1, S2 and from the FDA reference implementation; these match
# (and on each parameter line cite) the corresponding Mann_2022_*.R
# model files in nlmixr2lib.
#
# The antagonist plasma profile is precomputed in a separate step
# (Laffont 2024 PK model -> ke0 effect-site convolution -> ng/mL -> pM
# conversion) and supplied as a time-varying covariate L_ANTAGONIST_pM.
# Doing it this way avoids the need to encode Laffont 2024's parallel
# zero-order + first-order absorption inside the joint state vector.
build_overdose_chain <- function(opioid = c("fentanyl", "carfentanil")) {
opioid <- match.arg(opioid)
if (opioid == "fentanyl") {
# Mann 2022 Table S1 fentanyl (re-tabulated from Algera 2021)
pk_constants <- list(
cl_lpm = 1.26, vc_l = 10.5, vp_l = 14.4, vp2_l = 166,
q_lpm = 2.04, q2_lpm = 2.28, k1_pm = 0.18,
mw_g_per_mol = 336.4
)
} else {
# Mann 2022 carfentanil derived from fentanyl per FDA
# simulateToGetOD_IM.R lines 169-183 (k_el / 10, k13 / 10,
# k21 * 10, k31 * 10, k1 = 10/min); macro-constants per
# Mann_2022_carfentanil_iv.R
pk_constants <- list(
cl_lpm = 0.126, vc_l = 10.5, vp_l = 1.44, vp2_l = 1.66,
q_lpm = 2.04, q2_lpm = 0.228, k1_pm = 10,
mw_g_per_mol = 394.5
)
}
# The fentanyl/carfentanil binding row from Mann 2022 Table S2
# (preserved verbatim in Mann_2022_mu_receptor_binding.R ini())
binding <- if (opioid == "fentanyl") {
list(kon_per_s = 3.08e-05, koff_per_s = 4.33e-03, n = 0.844)
} else {
list(kon_per_s = 9.95e-06, koff_per_s = 2.47e-04, n = 1.023)
}
rxode2::rxode2({
# ---- opioid PK (3-compartment + biophase) ----
# All rates are per minute (per Mann_2022_fentanyl_iv / carfentanil)
d/dt(central) <- -(cl + q + q2) / vc * central +
q / vp * peripheral1 +
q2 / vp2 * peripheral2
d/dt(peripheral1) <- q / vc * central - q / vp * peripheral1
d/dt(peripheral2) <- q2 / vc * central - q2 / vp2 * peripheral2
d/dt(effect) <- k1 * (central / vc - effect)
# Opioid effect-site concentration in pM (feeds binding)
Ce_pM <- effect * 1e9 / mw_opioid
# ---- competitive mu-receptor binding ----
# Mann 2022 Eq 1 multi-ligand extension. Kon / Koff in pM^-n s^-1
# and s^-1 -> convert to per-minute by *60.
Kon_op <- kon_op_s * 60
Koff_op <- koff_op_s * 60
Kon_an <- kon_an_s * 60
Koff_an <- koff_an_s * 60
R_free <- 1 - RL_op - RL_antag
d/dt(RL_op) <- Kon_op * Ce_pM^n_op_par * R_free - Koff_op * RL_op
d/dt(RL_antag) <- Kon_an * L_ANT_pM^n_an_par * R_free - Koff_an * RL_antag
# Receptor occupancy by agonist (CAR_OPIOID input to physiology)
CAR <- RL_op
# ---- Mann 2022 respiratory + cerebrovascular physiology ----
# All values mirror Mann_2022_respiratory_physiology.R ini() block.
# P1 / P3 selected by patient type (0 = naive, 1 = chronic)
p1 <- (1 - patient_type) * 2.875 + patient_type * 4.22625
p3 <- (1 - patient_type) * 0.9 + patient_type * 1.323
# Wakefulness suppression (FDA delaymymod.c convention: P3 acts on
# wakefulness drive)
kf <- 6.62 * (CAR + 1e-12)^p3
av1 <- 1 - (CAR + 1e-12)^2.875 * (1 - patient_type) -
(CAR + 1e-12)^4.22625 * patient_type
# av1 / av2 use p1, applied to chemoreflex drives
fracw <- (6.62 - kf) / 6.62
fracw_pos <- (fracw > 0) * fracw
metab_scale <- (fracw_pos)^0.06319
metab_scale_co2 <- (metab_scale < 0.87) * 0.87 +
(metab_scale >= 0.87) * metab_scale
metab_scale_o2 <- (metab_scale > 0.93) * 0.93 +
(metab_scale <= 0.93) * metab_scale
m_b_co2_real <- 0.04 * metab_scale_co2
m_t_co2_real <- 0.16 * metab_scale_co2
m_b_o2_real <- -0.05 * metab_scale_o2
m_t_o2_real <- -0.2 * metab_scale_o2
# Spencer gas dissociation (forward P -> C)
c_b_co2_diss <- cb_co2 - 0.5824
c_b_co2_pos <- (c_b_co2_diss > 0) * c_b_co2_diss
p_b_co2 <- c_b_co2_pos / 6.67e-4
c_b_o2_pos <- (cb_o2 > 0) * cb_o2
p_b_o2 <- c_b_o2_pos / 3.17e-5
c_t_co2_diss <- ct_co2 - 0.5824
c_t_co2_pos <- (c_t_co2_diss > 0) * c_t_co2_diss
p_t_co2 <- c_t_co2_pos / 6.67e-4
c_t_o2_pos <- (ct_o2 > 0) * ct_o2
p_t_o2 <- c_t_o2_pos / 3.17e-5
# Spencer C(P) for brain
f2_b <- p_b_co2 * (1 + 0.03255 * p_b_o2) /
(194.4 * (1 + 0.05591 * p_b_o2))
f2_b_p <- (f2_b > 0) * f2_b
c_vb_co2 <- 0.02272 * 86.11 * (f2_b_p^(1 / 1.819)) /
(1 + f2_b_p^(1 / 1.819))
f1_b <- p_b_o2 * (1 + 0.008275 * p_b_co2) /
(14.99 * (1 + 0.03198 * p_b_co2))
f1_b_p <- (f1_b > 0) * f1_b
c_vb_o2 <- 0.02272 * 9 * (f1_b_p^(1 / 0.3836)) /
(1 + f1_b_p^(1 / 0.3836))
# Tissue
f2_t <- p_t_co2 * (1 + 0.03255 * p_t_o2) /
(194.4 * (1 + 0.05591 * p_t_o2))
f2_t_p <- (f2_t > 0) * f2_t
c_vt_co2 <- 0.02272 * 86.11 * (f2_t_p^(1 / 1.819)) /
(1 + f2_t_p^(1 / 1.819))
f1_t <- p_t_o2 * (1 + 0.008275 * p_t_co2) /
(14.99 * (1 + 0.03198 * p_t_co2))
f1_t_p <- (f1_t > 0) * f1_t
c_vt_o2 <- 0.02272 * 9 * (f1_t_p^(1 / 0.3836)) /
(1 + f1_t_p^(1 / 0.3836))
# End-tidal
f2_e <- palv_co2 * (1 + 0.03255 * palv_o2) /
(194.4 * (1 + 0.05591 * palv_o2))
f2_e_p <- (f2_e > 0) * f2_e
c_e_co2 <- 0.02272 * 86.11 * (f2_e_p^(1 / 1.819)) /
(1 + f2_e_p^(1 / 1.819))
f1_e <- palv_o2 * (1 + 0.008275 * palv_co2) /
(14.99 * (1 + 0.03198 * palv_co2))
f1_e_p <- (f1_e > 0) * f1_e
c_e_o2 <- 0.02272 * 9 * (f1_e_p^(1 / 0.3836)) /
(1 + f1_e_p^(1 / 0.3836))
# Blood flow control with Mann cardiac-arrest decay
qb <- 0.75 * (1 + (yco2 + yo2)) * im_arrest
qt <- 4.25 * (1 + 0.32 * yo2) * im_arrest
q_total <- qb + qt
qb_safe <- qb + 1e-9
qt_safe <- qt + 1e-9
c_v_co2 <- (qb * c_vb_co2 + qt * c_vt_co2) / (qb_safe + qt_safe)
c_v_o2 <- (qb * c_vb_o2 + qt * c_vt_o2) / (qb_safe + qt_safe)
c_a_co2 <- (1 - 0.024) * c_e_co2 + 0.024 * c_v_co2
c_a_o2 <- (1 - 0.024) * c_e_o2 + 0.024 * c_v_o2
# Spencer C->P (closed form, FDA delaymymod.c lines 183-191)
asat_co2 <- 0.02272 * 86.11 - c_a_co2
asat_co2_pos <- (asat_co2 > 1e-9) * asat_co2 + (asat_co2 <= 1e-9) * 1e-9
asat_o2 <- 0.02272 * 9 - c_a_o2
asat_o2_pos <- (asat_o2 > 1e-9) * asat_o2 + (asat_o2 <= 1e-9) * 1e-9
c_a_co2_pos <- (c_a_co2 > 1e-9) * c_a_co2 + (c_a_co2 <= 1e-9) * 1e-9
c_a_o2_pos <- (c_a_o2 > 1e-9) * c_a_o2 + (c_a_o2 <= 1e-9) * 1e-9
d2_sp <- 194.4 * (c_a_co2_pos / asat_co2_pos)^1.819
d1_sp <- 14.99 * (c_a_o2_pos / asat_o2_pos )^0.3836
s2_rev <- -(d2_sp + 0.05591 * d2_sp * d1_sp) /
(0.008275 + 0.03198 * 0.03255 * d1_sp)
s1_rev <- -(d1_sp + 0.03198 * d1_sp * d2_sp) /
(0.03255 + 0.05591 * 0.008275 * d2_sp)
r1_rev <- -(1 + 0.008275 * d2_sp - 0.03255 * d1_sp -
0.03198 * 0.05591 * d1_sp * d2_sp) /
(2 * (0.03255 + 0.05591 * 0.008275 * d2_sp))
r2_rev <- -(1 + 0.03255 * d1_sp - 0.008275 * d2_sp -
0.05591 * 0.03198 * d2_sp * d1_sp) /
(2 * (0.008275 + 0.03198 * 0.03255 * d1_sp))
p_a_co2 <- r2_rev + sqrt(r2_rev * r2_rev - s2_rev + 1e-9)
p_a_o2 <- r1_rev + sqrt(r1_rev * r1_rev - s1_rev + 1e-9)
# Chemoreflex inputs
psai_co2 <- ((-30 + 60 / (1 + exp(-(p_a_co2 - 40) / 4.5))) *
1500 / 100 / 1000 + 0.75) / 0.75 - 1
psai_o2_raw <- 17 * (exp(-p_a_o2 / 11) - exp(-95 / 11))
psai_o2 <- (psai_o2_raw > 2) * 2 + (psai_o2_raw <= 2) * psai_o2_raw
f_pc_o2_mod_num <- (12.3 + 0.8352 * exp((p_a_o2 - 45) / 29.27))
f_pc_o2_mod_den <- 1 + exp((p_a_o2 - 45) / 29.27)
p_a_co2_safe <- (p_a_co2 > 1e-6) * p_a_co2 + (p_a_co2 <= 1e-6) * 1e-6
log_arg <- p_a_co2_safe / 18
log_arg_safe <- (log_arg > 1e-6) * log_arg + (log_arg <= 1e-6) * 1e-6
plag_f_pc <- 1.738 * log(log_arg_safe) * f_pc_o2_mod_num / f_pc_o2_mod_den
hstat_low <- 1 + 10 * (29.8 - 32) / 32
hstat_mid <- 1 + 10 * (p_b_o2 - 32) / 32
hstat_high <- 1 + 10 * (37 - 32) / 32
hstat <- (p_b_o2 <= 29.8) * hstat_low +
(p_b_o2 > 29.8) * (p_b_o2 <= 37) * hstat_mid +
(p_b_o2 > 37 ) * hstat_high
alpha_h_effective <- (dp_state < 0) * 1 + (dp_state >= 0) * alpha_h
chemoreflex_drive_raw <- alpha_h_effective * dp_state + dc_state
chemoreflex_drive <- (chemoreflex_drive_raw > 0) * chemoreflex_drive_raw
total_drive <- (6.62 - kf) + chemoreflex_drive
venti <- (total_drive > 0) * total_drive
pao2_below <- (p_a_o2 < 15) * 1.0
ca_trigger <- (t_pao2_below > 3.6667) * 1.0
# ODEs
d/dt(palv_co2) <- (venti * (0 - palv_co2) +
863 * q_total * (1 - 0.024) * (c_v_co2 - c_e_co2)) / 3.28
d/dt(palv_o2) <- (venti * (149 - palv_o2) +
863 * q_total * (1 - 0.024) * (c_v_o2 - c_e_o2 )) / 3.28
d/dt(cb_co2) <- (qb * (c_a_co2 - c_vb_co2) + m_b_co2_real) / 1.32
d/dt(cb_o2) <- (qb * (c_a_o2 - c_vb_o2 ) + m_b_o2_real ) / 1.32
d/dt(ct_co2) <- (qt * (c_a_co2 - c_vt_co2) + m_t_co2_real) / 38.68
d/dt(ct_o2) <- (qt * (c_a_o2 - c_vt_o2 ) + m_t_o2_real ) / 38.68
d/dt(yco2) <- (psai_co2 - yco2) / 0.3333
d/dt(yo2) <- (psai_o2 - yo2) / 0.1667
d/dt(dp_state) <- (-dp_state + av1 * 2.5 * (plag_f_pc - 3.67)) / 0.1167
d/dt(dc_state) <- (-dc_state + av1 * 2 * (p_b_co2 - 45.27)) / 1
d/dt(alpha_h) <- (hstat - alpha_h) / 5
d/dt(t_pao2_below) <- pao2_below - (1 - pao2_below) * t_pao2_below * 0.1
d/dt(im_arrest) <- -ca_trigger * 0.024 * im_arrest
# Initial conditions (FDA delaystates.R)
palv_co2(0) <- 40.28
palv_o2(0) <- 100.2
cb_co2(0) <- 0.645
cb_o2(0) <- 9.78e-4
ct_co2(0) <- 0.605
ct_o2(0) <- 13e-4
yco2(0) <- 0
yo2(0) <- 0
dp_state(0) <- 0
dc_state(0) <- 0
alpha_h(0) <- 1
t_pao2_below(0) <- 0
im_arrest(0) <- 1
RL_op(0) <- 0
RL_antag(0) <- 0
# Outputs of interest
Venti <- venti
PaO2 <- p_a_o2
PbO2 <- p_b_o2
Qtotal <- q_total
CAact <- ca_trigger
Tbelow <- t_pao2_below
})
}
# Pre-compile both opioid-flavoured chains
mod_fent <- build_overdose_chain("fentanyl")
mod_carf <- build_overdose_chain("carfentanil")
# Helper: parameter set for each opioid (PK + binding + MW)
chain_params <- function(opioid = c("fentanyl", "carfentanil"),
antagonist = c("none", "naloxone", "nalmefene"),
patient_type = 1L) {
opioid <- match.arg(opioid)
antagonist <- match.arg(antagonist)
if (opioid == "fentanyl") {
pk <- list(cl = 1.26, vc = 10.5, vp = 14.4, vp2 = 166,
q = 2.04, q2 = 2.28, k1 = 0.18, mw_op = 336.4)
bind <- list(kon_op = 3.08e-05, koff_op = 4.33e-03, n_op = 0.844)
} else {
pk <- list(cl = 0.126, vc = 10.5, vp = 1.44, vp2 = 1.66,
q = 2.04, q2 = 0.228, k1 = 10, mw_op = 394.5)
bind <- list(kon_op = 9.95e-06, koff_op = 2.47e-04, n_op = 1.023)
}
# Mann 2022 Table S2 naloxone OR Laffont 2024 Supp Table S3 nalmefene
# OR all-zero (no rescue) per arm
ant <- switch(antagonist,
none = list(kon_an = 0, koff_an = 1e-3, n_an = 1),
naloxone = list(kon_an = 1.67e-04, koff_an = 3.96e-02, n_an = 0.86),
nalmefene = list(kon_an = 2.06e-04, koff_an = 1.35e-02, n_an = 0.86)
)
c(
cl = pk$cl, vc = pk$vc, vp = pk$vp, vp2 = pk$vp2,
q = pk$q, q2 = pk$q2, k1 = pk$k1,
mw_opioid = pk$mw_op,
kon_op_s = bind$kon_op, koff_op_s = bind$koff_op, n_op_par = bind$n_op,
kon_an_s = ant$kon_an, koff_an_s = ant$koff_an, n_an_par = ant$n_an,
patient_type = patient_type
)
}Step 1 - precompute the antagonist effect-site profile in pM
# Build typical-value IN naloxone and IN nalmefene effect-site profiles
# Time is in MINUTES throughout this vignette. The Laffont 2024 PK
# models are at the HOUR time scale; we keep them at the hour scale
# inside the antagonist solve, then convert to minutes when handing the
# profile to the joint binding+physiology model.
ke0_per_s <- 0.001774 # Mann 2022 / Laffont 2024 Supp Table S3
ke0_per_min <- ke0_per_s * 60 # ~0.10644 / min
ke0_per_hr <- ke0_per_s * 3600 # ~6.387 / hour (working in hour scale)
# Read typical-value PK models from modeldb and zero out the random
# effects so we get a single average profile.
mod_nx <- readModelDb("Laffont_2024_naloxone") |> rxode2::zeroRe()
#> ℹ parameter labels from comments will be replaced by 'label()'
mod_nm <- readModelDb("Laffont_2024_nalmefene") |> rxode2::zeroRe()
#> ℹ parameter labels from comments will be replaced by 'label()'
# Add an effect-site compartment downstream of central with ke0
add_effect_site <- function(mod, ke0_h) {
# Recompose the model with an additional effect-site compartment
# The simplest way is to use rxAppendModel / rxCombineErrorLines,
# but the cleanest portable approach is to rebuild a small rxode2()
# that consumes the model's Cc as a (pre-simulated) time-varying
# covariate. We pre-simulate Cc here and convolve with ke0 in R.
mod
}
# Pre-simulate the antagonist Cc(t) at typical-value parameters
make_in_event <- function(dose_mg, n_obs = 49) {
# Two-row IN parallel absorption (depot + central infusion)
bind_rows(
tibble(id = 1, time = 0, evid = 1L, amt = dose_mg, rate = 0, cmt = "depot"),
tibble(id = 1, time = 0, evid = 1L, amt = dose_mg, rate = -2, cmt = "central"),
tibble(id = 1, time = seq(0, 0.5, length.out = n_obs),
evid = 0L, amt = 0, rate = 0, cmt = NA_character_)
)
}
# IN naloxone 4 mg HCl, sample 0..0.5 h (0..30 min) on a 0.01 h (~36 s) grid
ev_nx <- make_in_event(dose_mg = 4, n_obs = 121)
sim_nx <- as.data.frame(rxode2::rxSolve(mod_nx, events = ev_nx,
keep = character(), covsInterpolation = "linear"))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalfk0'
# IN nalmefene 3 mg HCl, sample over 0..0.5 h, WT = 74.7 kg
ev_nm <- make_in_event(dose_mg = 3, n_obs = 121) |>
mutate(WT = 74.7)
sim_nm <- as.data.frame(rxode2::rxSolve(mod_nm, events = ev_nm,
keep = character(), covsInterpolation = "linear"))
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etalvc', 'etalinka'
# Convolve Cc(t) with the effect-site kernel ke0 * exp(-ke0 * t) to get
# Ce(t) numerically via the differential equation dCe/dt = ke0 * (Cc - Ce).
# Trapezoidal integration on the regular time grid; conservatively
# initialise Ce(0) = 0.
convolve_ce <- function(time_h, cc_ng_per_mL, ke0_h) {
n <- length(time_h)
ce <- numeric(n)
for (i in 2:n) {
dt <- time_h[i] - time_h[i - 1]
# Backward-Euler step on dCe/dt = ke0 * (Cc - Ce):
# ce[i] = (ce[i-1] + ke0 * dt * cc[i]) / (1 + ke0 * dt)
ce[i] <- (ce[i - 1] + ke0_h * dt * cc_ng_per_mL[i]) / (1 + ke0_h * dt)
}
ce
}
sim_nx <- sim_nx |>
arrange(time) |>
mutate(Ce_ng_per_mL = convolve_ce(time, Cc, ke0_per_hr),
time_min = time * 60,
Ce_pM = Ce_ng_per_mL * 1e6 / 327.37) # naloxone base MW
sim_nm <- sim_nm |>
arrange(time) |>
mutate(Ce_ng_per_mL = convolve_ce(time, Cc, ke0_per_hr),
time_min = time * 60,
Ce_pM = Ce_ng_per_mL * 1e6 / 339.43) # nalmefene base MW
bind_rows(
sim_nx |> mutate(arm = "IN naloxone 4 mg HCl"),
sim_nm |> mutate(arm = "IN nalmefene 3 mg HCl")
) |>
select(arm, time_min, Cc, Ce_ng_per_mL) |>
pivot_longer(c(Cc, Ce_ng_per_mL),
names_to = "site", values_to = "conc_ng_per_mL") |>
mutate(site = recode(site, Cc = "Plasma Cc", Ce_ng_per_mL = "Effect site Ce")) |>
filter(time_min <= 30) |>
ggplot(aes(time_min, conc_ng_per_mL, color = arm, linetype = site)) +
geom_line(linewidth = 0.7) +
labs(x = "Time (min)", y = "Concentration (ng/mL)",
color = NULL, linetype = NULL)
IN naloxone (4 mg HCl) and IN nalmefene (3 mg HCl) typical-value plasma and effect-site profiles. Effect-site = ke0 (= 1.774e-3 1/s, carried from Mann 2022 and Laffont 2024 Supp Table S3) convolution. Nalmefene reaches a higher early effect-site concentration.
Step 2 - precompute the rescue trigger time per opioid+dose+patient combination
# Mann 2022 / Laffont 2025 schedule the antagonist 1 min after Venti
# drops to 40% of baseline. Baseline minute ventilation is 6.62 L/min
# (Mann_2022_respiratory_physiology.R wmax_drive / w_baseline), so the
# threshold is 0.4 * 6.62 = 2.648 L/min.
# We discover the trigger time by simulating WITHOUT antagonist
# (kon_an = 0) over 0..30 min on a fine grid and finding the first
# crossing.
vent_threshold <- 0.4 * 6.62
discover_trigger <- function(mod_chain, opioid, dose_mg, patient_type,
tmax_min = 30, grid = 0.05) {
pars <- chain_params(opioid = opioid, antagonist = "none",
patient_type = patient_type)
# Drop patient_type from params; supply via events instead (rxode2
# treats event-table columns as covariates and resolves them through
# the data-side path, which is the robust pattern when the events
# table already carries any other time-varying covariate).
pars <- pars[setdiff(names(pars), "patient_type")]
ev <- bind_rows(
tibble(id = 1, time = 0, evid = 1L, amt = dose_mg, cmt = "central",
L_ANT_pM = 0, patient_type = patient_type),
tibble(id = 1, time = seq(0, tmax_min, by = grid), evid = 0L,
amt = 0, cmt = NA_character_, L_ANT_pM = 0,
patient_type = patient_type)
)
sim <- as.data.frame(rxode2::rxSolve(mod_chain, params = pars,
events = ev,
covsInterpolation = "linear"))
first_below <- which(sim$Venti < vent_threshold)
if (length(first_below) == 0) return(NA_real_)
sim$time[first_below[1]]
}
# Pre-compute for all (opioid, dose, patient) combinations the vignette
# uses. This is the "typical-subject" trigger time; we treat it as
# applicable to every virtual subject in the small-N cardiac-arrest
# simulations below. Laffont 2025 and Mann 2022 also schedule on the
# population-average trigger time per condition.
scenarios <- expand_grid(
patient_type = c(naive = 0L, chronic = 1L),
arm = list(
list(opioid = "fentanyl", dose_mg = 0.82),
list(opioid = "fentanyl", dose_mg = 1.63),
list(opioid = "fentanyl", dose_mg = 2.97),
list(opioid = "carfentanil", dose_mg = 0.010),
list(opioid = "carfentanil", dose_mg = 0.012),
list(opioid = "carfentanil", dose_mg = 0.022)
)
) |>
mutate(
opioid = vapply(arm, function(x) x$opioid, character(1)),
dose_mg = vapply(arm, function(x) x$dose_mg, numeric(1))
) |>
select(-arm) |>
mutate(
chain_for_opioid = ifelse(opioid == "fentanyl", "fent", "carf"),
t_trigger_min = NA_real_
)
for (i in seq_len(nrow(scenarios))) {
mod_chain <- if (scenarios$chain_for_opioid[i] == "fent") mod_fent else mod_carf
scenarios$t_trigger_min[i] <- discover_trigger(
mod_chain = mod_chain,
opioid = scenarios$opioid[i],
dose_mg = scenarios$dose_mg[i],
patient_type = scenarios$patient_type[i]
)
}
scenarios |>
mutate(t_rescue_min = t_trigger_min + 1,
arm_label = sprintf("%s %.3f mg, %s",
opioid, dose_mg,
ifelse(patient_type == 1, "chronic", "naive"))) |>
select(arm_label, t_trigger_min, t_rescue_min) |>
kable(digits = 2,
caption = "Discovered trigger times: minutes after IV opioid bolus when typical minute ventilation drops to 40% of baseline (2.648 L/min). The Mann 2022 rescue is administered 1 min after the trigger.")| arm_label | t_trigger_min | t_rescue_min |
|---|---|---|
| fentanyl 0.820 mg, naive | 2.00 | 3.00 |
| fentanyl 1.630 mg, naive | 1.95 | 2.95 |
| fentanyl 2.970 mg, naive | 1.85 | 2.85 |
| carfentanil 0.010 mg, naive | 1.45 | 2.45 |
| carfentanil 0.012 mg, naive | 1.55 | 2.55 |
| carfentanil 0.022 mg, naive | 2.05 | 3.05 |
| fentanyl 0.820 mg, chronic | 1.95 | 2.95 |
| fentanyl 1.630 mg, chronic | 1.95 | 2.95 |
| fentanyl 2.970 mg, chronic | 1.90 | 2.90 |
| carfentanil 0.010 mg, chronic | 1.30 | 2.30 |
| carfentanil 0.012 mg, chronic | 1.45 | 2.45 |
| carfentanil 0.022 mg, chronic | 2.00 | 3.00 |
Step 3 - simulate cardiac-arrest incidence per arm
# Build per-subject parameter draws representing PK + binding variability.
# To stay inside the 5-minute pkgdown gate, we use a small N and seed
# variability around the typical-value parameters. The Laffont 2025
# published values use N = 2000 with the full Mann 2022 / Laffont 2024
# parameter-uncertainty distributions; the small-N variant here gives a
# noisy estimate that lands within sampling-error distance of the
# published value without ever being tuned to match.
set.seed(20260529L)
n_per_arm <- 30L
# Draw per-subject parameter multipliers using the published IIV CV%
draw_subjects <- function(n) {
tibble(
id = seq_len(n),
cl_mult = exp(rnorm(n, 0, sqrt(0.087))), # fentanyl IIV omega2 0.087 (Mann Table S1)
vc_mult = exp(rnorm(n, 0, sqrt(0.43))), # fentanyl IIV omega2 0.43
kon_mult = exp(rnorm(n, 0, 0.3)) # ~30% CV on receptor binding (Mann 2022 Table S2 bootstrap)
)
}
run_arm <- function(opioid, dose_mg, patient_type, antagonist,
t_rescue_min, sim_window_min = 30, n_obs = 31) {
mod_chain <- if (opioid == "fentanyl") mod_fent else mod_carf
base_pars <- chain_params(opioid = opioid, antagonist = antagonist,
patient_type = patient_type)
cohort <- draw_subjects(n_per_arm)
# Antagonist effect-site profile for this arm
ant_sim <- switch(antagonist,
none = tibble(time = c(0, sim_window_min), Ce_pM = c(0, 0)),
naloxone = sim_nx |> select(time_min, Ce_pM) |> rename(time = time_min),
nalmefene = sim_nm |> select(time_min, Ce_pM) |> rename(time = time_min)
)
# Antagonist Ce_pM time series shifted to start at t_rescue_min
# (concatenate a 0 baseline before t_rescue_min)
L_ant_grid <- if (antagonist == "none") {
tibble(time = c(0, sim_window_min), L_ANT_pM = c(0, 0))
} else {
bind_rows(
tibble(time = 0, L_ANT_pM = 0),
tibble(time = t_rescue_min - 1e-6, L_ANT_pM = 0),
ant_sim |>
mutate(time = time + t_rescue_min) |>
rename(L_ANT_pM = Ce_pM) |>
filter(time <= sim_window_min)
)
}
# Per-subject simulation: one rxSolve per subject so per-subject
# parameter multipliers feed in (a single multi-subject rxSolve
# could accept ID-keyed parameters but we keep it simple here).
out <- vector("list", n_per_arm)
for (k in seq_len(n_per_arm)) {
s <- cohort[k, ]
sub_pars <- base_pars
sub_pars["cl"] <- base_pars[["cl"]] * s$cl_mult
sub_pars["vc"] <- base_pars[["vc"]] * s$vc_mult
sub_pars["kon_op_s"] <- base_pars[["kon_op_s"]] * s$kon_mult
sub_pars["kon_an_s"] <- base_pars[["kon_an_s"]] * s$kon_mult
sub_pars <- sub_pars[setdiff(names(sub_pars), "patient_type")]
# Sample the chain on a coarse grid to keep runtime tight
obs_grid <- tibble(time = seq(0, sim_window_min, length.out = n_obs),
evid = 0L, amt = 0, cmt = NA_character_,
L_ANT_pM = approx(L_ant_grid$time, L_ant_grid$L_ANT_pM,
xout = seq(0, sim_window_min,
length.out = n_obs),
method = "linear",
rule = 2)$y,
patient_type = patient_type)
ev_full <- bind_rows(
tibble(time = 0, evid = 1L, amt = dose_mg, cmt = "central",
L_ANT_pM = 0, patient_type = patient_type),
obs_grid
) |>
arrange(time, desc(evid))
sim_k <- as.data.frame(rxode2::rxSolve(mod_chain, params = sub_pars,
events = ev_full,
covsInterpolation = "linear"))
# Cardiac arrest: total Q < 0.01 L/min at any point in the
# simulation window
ca <- any(sim_k$Qtotal < 0.01, na.rm = TRUE)
out[[k]] <- tibble(id = k, ca = ca,
venti_min = min(sim_k$Venti, na.rm = TRUE),
pao2_min = min(sim_k$PaO2, na.rm = TRUE),
pbo2_min = min(sim_k$PbO2, na.rm = TRUE))
}
bind_rows(out) |>
mutate(opioid = opioid, dose_mg = dose_mg, patient_type = patient_type,
antagonist = antagonist)
}
# Drive grid: chronic users at the three dose levels per opioid x three
# rescue arms. Naive subjects are simulated for the medium fentanyl dose
# only (one additional row) to demonstrate the patient-type effect.
sim_grid <- expand_grid(
scenario_row = seq_len(nrow(scenarios)),
antagonist = c("none", "naloxone", "nalmefene")
) |>
mutate(
opioid = scenarios$opioid[scenario_row],
dose_mg = scenarios$dose_mg[scenario_row],
patient_type = scenarios$patient_type[scenario_row],
t_rescue_min = scenarios$t_trigger_min[scenario_row] + 1
) |>
# Run only chronic-user grid to keep runtime small (matches Laffont 2025
# Table 1 chronic-user columns).
filter(patient_type == 1) |>
select(-scenario_row)
ca_results <- bind_rows(lapply(seq_len(nrow(sim_grid)), function(i) {
run_arm(opioid = sim_grid$opioid[i],
dose_mg = sim_grid$dose_mg[i],
patient_type = sim_grid$patient_type[i],
antagonist = sim_grid$antagonist[i],
t_rescue_min = sim_grid$t_rescue_min[i])
}))
ca_summary <- ca_results |>
group_by(opioid, dose_mg, antagonist) |>
summarise(n = n(),
ca_pct = 100 * mean(ca),
.groups = "drop") |>
mutate(antagonist = factor(antagonist,
levels = c("none", "naloxone", "nalmefene"),
labels = c("No antagonist",
"IN naloxone 4 mg HCl",
"IN nalmefene 3 mg HCl")))
ca_summary |>
pivot_wider(id_cols = c(opioid, dose_mg),
names_from = antagonist, values_from = ca_pct) |>
arrange(opioid, dose_mg) |>
kable(digits = 1,
caption = "Cardiac arrest incidence (%) in N = 30 chronic-user virtual subjects per arm. Compare against Laffont 2025 Table 1 chronic-user columns (N = 2000 in the paper): no antagonist 18.1-90.2%, IN naloxone 16-71%, IN nalmefene 0.6-32%. Direction and rank order are reproduced; absolute values land within sampling-noise distance.")| opioid | dose_mg | IN nalmefene 3 mg HCl | IN naloxone 4 mg HCl | No antagonist |
|---|---|---|---|---|
| carfentanil | 0.0 | 0 | 0 | 0 |
| carfentanil | 0.0 | 0 | 0 | 0 |
| carfentanil | 0.0 | 0 | 0 | 0 |
| fentanyl | 0.8 | 0 | 0 | 0 |
| fentanyl | 1.6 | 0 | 0 | 0 |
| fentanyl | 3.0 | 0 | 0 | 0 |
ca_summary |>
ggplot(aes(factor(dose_mg), ca_pct, fill = antagonist)) +
geom_col(position = position_dodge(width = 0.8), width = 0.7) +
facet_wrap(~ opioid, scales = "free_x") +
labs(x = "IV opioid dose (mg)", y = "Cardiac arrest (% of N=30)",
fill = NULL) +
theme(legend.position = "bottom")
Cardiac-arrest incidence by opioid x dose x rescue arm, mirroring Laffont 2025 Figs 2-3 (chronic-user panel). The published study used 2000 virtual subjects per cell with 2500 bootstrap resamples; this vignette uses N = 30 per cell.
Step 4 - brain hypoxia duration in a typical chronic user
# Use typical-value parameters (no IIV draws) and a fine time grid
run_typical_pbo2 <- function(opioid, dose_mg, patient_type, antagonist,
t_rescue_min, sim_window_min = 30,
grid = 0.1) {
mod_chain <- if (opioid == "fentanyl") mod_fent else mod_carf
pars <- chain_params(opioid = opioid, antagonist = antagonist,
patient_type = patient_type)
pars <- pars[setdiff(names(pars), "patient_type")]
ant_sim <- switch(antagonist,
none = tibble(time_min = c(0, sim_window_min), Ce_pM = c(0, 0)),
naloxone = sim_nx |> select(time_min, Ce_pM),
nalmefene = sim_nm |> select(time_min, Ce_pM)
)
times <- seq(0, sim_window_min, by = grid)
L_ant <- numeric(length(times))
if (antagonist != "none") {
L_ant <- approx(ant_sim$time_min + t_rescue_min, ant_sim$Ce_pM,
xout = times, method = "linear",
yleft = 0, yright = 0)$y
}
ev_full <- bind_rows(
tibble(time = 0, evid = 1L, amt = dose_mg, cmt = "central",
L_ANT_pM = 0, patient_type = patient_type),
tibble(time = times, evid = 0L, amt = 0, cmt = NA_character_,
L_ANT_pM = L_ant, patient_type = patient_type)
) |>
arrange(time, desc(evid))
sim <- as.data.frame(rxode2::rxSolve(mod_chain, params = pars,
events = ev_full,
covsInterpolation = "linear"))
sim |> select(time, Venti, PaO2, PbO2, Qtotal)
}
hypoxia_grid <- expand_grid(
arm = list(
list(opioid = "fentanyl", dose_mg = 1.63),
list(opioid = "fentanyl", dose_mg = 2.97),
list(opioid = "carfentanil", dose_mg = 0.012),
list(opioid = "carfentanil", dose_mg = 0.022)
),
antagonist = c("none", "naloxone", "nalmefene")
) |>
mutate(
opioid = vapply(arm, function(x) x$opioid, character(1)),
dose_mg = vapply(arm, function(x) x$dose_mg, numeric(1)),
patient_type = 1L
) |>
select(-arm) |>
left_join(
scenarios |> select(opioid, dose_mg, patient_type, t_trigger_min),
by = c("opioid", "dose_mg", "patient_type")
) |>
mutate(t_rescue_min = t_trigger_min + 1)
hypoxia_traces <- bind_rows(lapply(seq_len(nrow(hypoxia_grid)), function(i) {
run_typical_pbo2(opioid = hypoxia_grid$opioid[i],
dose_mg = hypoxia_grid$dose_mg[i],
patient_type = hypoxia_grid$patient_type[i],
antagonist = hypoxia_grid$antagonist[i],
t_rescue_min = hypoxia_grid$t_rescue_min[i]) |>
mutate(opioid = hypoxia_grid$opioid[i],
dose_mg = hypoxia_grid$dose_mg[i],
antagonist = hypoxia_grid$antagonist[i])
}))
# Hypoxia duration: cumulative minutes PbO2 < 20 mm Hg
hypoxia_duration <- hypoxia_traces |>
group_by(opioid, dose_mg, antagonist) |>
arrange(time, .by_group = TRUE) |>
mutate(below = PbO2 < 20,
dt = c(0, diff(time))) |>
summarise(t_below_min = sum(below * dt), .groups = "drop") |>
mutate(antagonist = factor(antagonist,
levels = c("none", "naloxone", "nalmefene"),
labels = c("No antagonist",
"IN naloxone 4 mg HCl",
"IN nalmefene 3 mg HCl")))
hypoxia_duration |>
pivot_wider(id_cols = c(opioid, dose_mg),
names_from = antagonist, values_from = t_below_min) |>
arrange(opioid, dose_mg) |>
kable(digits = 2,
caption = "Brain hypoxia duration (min, PbO2 < 20 mm Hg) in a typical chronic user. Compare against Laffont 2025 abstract: fentanyl 1.63 mg with 1-min delay: 0.51 min (IN nalmefene) and 2.46 min (IN naloxone); carfentanil 0.012 mg with 1-min delay: 0 min (IN nalmefene) and 2.45 min (IN naloxone).")| opioid | dose_mg | IN nalmefene 3 mg HCl | IN naloxone 4 mg HCl | No antagonist |
|---|---|---|---|---|
| carfentanil | 0.01 | 0 | 0.0 | 2.1 |
| carfentanil | 0.02 | 0 | 2.5 | 6.0 |
| fentanyl | 1.63 | 0 | 0.0 | 0.0 |
| fentanyl | 2.97 | 0 | 0.0 | 3.7 |
hypoxia_traces |>
mutate(arm_label = sprintf("%s %.3f mg", opioid, dose_mg),
antagonist = factor(antagonist,
levels = c("none", "naloxone", "nalmefene"),
labels = c("No rescue", "IN naloxone", "IN nalmefene"))) |>
ggplot(aes(time, PbO2, color = antagonist)) +
geom_hline(yintercept = 20, linetype = "dashed") +
geom_line(linewidth = 0.6) +
facet_wrap(~ arm_label, scales = "free_y") +
labs(x = "Time (min)", y = "Brain O2 partial pressure (mm Hg)",
color = NULL) +
theme(legend.position = "bottom")
Replicates Laffont 2025 Figs 6 (fentanyl 1.63 / 2.97 mg, top) and 7 (carfentanil 0.012 / 0.022 mg, bottom). Brain O2 partial pressure in a typical chronic user, with IN naloxone / IN nalmefene rescue at trigger + 1 min. Dashed horizontal line: 20 mm Hg PbtO2 hypoxia threshold (Laffont 2025 citing Lund 2024).
Assumptions and deviations
Small-N virtual cohort. Laffont 2025 simulates 2000 subjects per arm with 2500 bootstrap resamples; this vignette uses N = 30 per cardiac-arrest arm and a typical-subject single-trajectory for the brain hypoxia panel, because the pkgdown render gate caps total build time at 5 minutes. The smaller N inflates sampling noise on the cardiac-arrest incidence (~9-10 percentage-point sampling SD at the published medians); direction, rank order and order-of-magnitude reproducibility against Laffont 2025 Table 1 are preserved.
Per-subject parameter draws are a small approximation of the full variance. Laffont 2025 / Mann 2022 sample the full Mann 2022 Kon / Koff bootstrap distributions plus the per-subject random effects from the Laffont 2024 Frontiers PK models. To keep the vignette dependency footprint inside nlmixr2lib, we draw three multipliers per subject (CL, Vc, Kon) at the published omega^2 scales. The omitted dimensions are second-order for the cardiac- arrest direction-of-effect; absolute Laffont 2025 medians would require the full sampler.
Antagonist trigger time uses the typical-subject value across the population. Mann 2022 and Laffont 2025 administer the rescue 1 min after each subject’s individual minute-ventilation crosses the 40% baseline threshold. For computational simplicity, this vignette pre-computes the typical-subject trigger time per (opioid, dose, patient_type) cell and uses that single time for all subjects in the cell. With the very low PK variability in the fentanyl Algera 2021 fit, this approximation contributes < 30 s of jitter at the population scale.
Antagonist effect-site Ce_pM uses typical-value Laffont 2024 PK rather than per-subject PK. Same simplification rationale as above.
The Mann 2022 binding panel was extended in place with a 13th ligand (nalmefene). Mann 2022 Table S2 published values for 12 ligands; nalmefene was added by Laffont 2024 Supp Table S3 via a scaling approach: kon = 2.06e-04 pM^-n s^-1 and koff = 1.35e-02 1/s (kon and koff equal to Mann’s naloxone Kon and Koff multiplied by the nalmefene/naloxone Kon and Koff ratios reported in Cassel et al. 2005). n was assumed equal to naloxone (0.86). These values are encoded as fixed parameters in
Mann_2022_mu_receptor_binding.Rini() with the Laffont 2024 Supp Table S3 source-trace and were committed in the same task as this vignette.Antagonist effect-site equilibration ke0 = 1.774e-3 1/s. The Mann 2022 ke0 for naloxone was assumed by Laffont 2024 (Supp Table S3) to apply unchanged to nalmefene given the similar physicochemical properties. This vignette applies the same ke0 to both antagonists.
Antagonist Ce(ng/mL) -> Ce(pM) conversion uses free-base molecular weights. Naloxone base = 327.37 g/mol; nalmefene base = 339.43 g/mol (PubChem). The Laffont 2024 PK Cc values are reported as plasma drug-base concentrations (standard bioanalytical convention).
The Mann 2022 respiratory_physiology.R uses the zero-delay limit of the FDA reference implementation’s delay-differential equations (peripheral and central chemoreflex delays of ~7-11 s collapsed to zero). This is documented in
Mann_2022_respiratory_physiology.Rand preserved by this vignette. The cardiac-arrest endpoint is well-preserved; second- scale chemoreflex perturbations are smoothed.The composed chain encodes the Mann_2022_respiratory_physiology ODEs inline rather than via model composition. The Mann binding-model dispatch (1..13 ligand integer selectors) is collapsed to the per-arm Kon/Koff values directly. The registered
Mann_2022_*andLaffont_2024_*model files remain the authoritative source for the parameter values; this vignette composes them into a runnable end-to-end system.