Skip to contents

Opioid overdose reversal simulation chain (Laffont 2025)

Source: vignettes/articles/Laffont_2025_opioid_overdose_reversal_simulation.Rmd
Laffont_2025_opioid_overdose_reversal_simulation.Rmd

What 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, and Mann_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_naloxone and Laffont_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) with Q_Scale shock feedback ----
    # FDA delaymymod.c lines 358-368: during opioid overdose, chemoreflex-
    # driven hyperperfusion raises cardiac output Q above baseline Q_0.
    # The effective central volume of distribution is scaled by Q_Scale,
    # concentrating the opioid in the biophase compartment. As Q/Q_0
    # rises from 1 to ~2 (Q ~8.5 L/min during hypoxic chemoreflex peak),
    # Q_Scale traverses a steep sigmoid centered at Q/Q_0 = 1.6, taking
    # values from 1 to 2 - effectively halving Vc and doubling Ce. This
    # is the positive feedback that, in FDA's published simulations,
    # drives fentanyl-overdose PaO2 below 15 mm Hg for the full 220 s
    # cardiac-arrest sustain window.
    # qtot_for_scale is q_total computed inline here (it is also computed
    # later in physiology block but rxode2 algebraic evaluation is in-order
    # so we re-derive it from the same state values).
    qtot_for_scale <- 0.75 * (1 + (yco2 + yo2)) * im_arrest +
                      4.25 * (1 + 0.32 * yo2)   * im_arrest
    q_scale_raw <- 1 + 1 / (1 + exp((1.6 - qtot_for_scale / 4.87) / 0.05))
    q_scale <- (q_scale_raw < 1) * 1 + (q_scale_raw > 2) * 2 +
               (q_scale_raw >= 1) * (q_scale_raw <= 2) * q_scale_raw
    vc_eff <- vc / q_scale
    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
    # FDA delaymymod.c line 594: ydot[8] = k1 * y[5] / VPC ... - k1*y[8].
    # The effect-site sees the SCALED central volume, not vc.
    d/dt(effect)      <-  k1 * (central / vc_eff - 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
    # FDA: M_B_o2 magnitude floored at fN * |baseline| -> metab_scale_o2 FLOOR at 0.93
    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

    # Low-PaO2 metabolism safety clamps (FDA delaymymod.c lines 296-306):
    # at P_X_o2 < 6 mm Hg scale both O2 consumption and CO2 production
    # by P_X_o2 / 6 so they go smoothly to zero.
    m_o2_safety_brain  <- (p_b_o2 < 6) * (p_b_o2 / 6) + (p_b_o2 >= 6) * 1.0
    m_o2_safety_tissue <- (p_t_o2 < 6) * (p_t_o2 / 6) + (p_t_o2 >= 6) * 1.0

    # 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
    # FDA delaymymod.c line 502: psai_co2 = (sigmoid * 1500/100/1000/60 + Qb0) / Qb0 - 1
    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

    # Chemoreflex transport delays (7 s peripheral, 11 s central) via
    # the nlmixr2lib-registered lagState C user function. Channels 3, 4, 5
    # used here so they do not collide with the standalone model's
    # 0, 1, 2 if both are compiled into the same rxode2 session.
    plag_p_a_o2  <- lagState(t, p_a_o2,  (7  / 60), 3, 95)
    plag_p_a_co2 <- lagState(t, p_a_co2, (7  / 60), 4, 40)
    clag_p_b_co2 <- lagState(t, p_b_co2, (11 / 60), 5, 45.27)

    f_pc_o2_mod_num <- (12.3 + 0.8352 * exp((plag_p_a_o2 - 45) / 29.27))
    f_pc_o2_mod_den <- 1 + exp((plag_p_a_o2 - 45) / 29.27)
    plag_p_a_co2_safe <- (plag_p_a_co2 > 1e-6) * plag_p_a_co2 + (plag_p_a_co2 <= 1e-6) * 1e-6
    log_arg <- plag_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

    # FDA-faithful one-shot cardiac-arrest scheme via state-dependent mtime.
    # Schedule a solver restart at first_cross_time + 3.6667 (= 220 s); the
    # mtime expression is dynamic: as first_cross_time evolves, the firing
    # time tracks. Once crossed_latch latches at the first PaO2 < 15 crossing,
    # first_cross_time freezes at T1 and the mtime converges to T1 + 220 s.
    # check_complete then closes the decision window so arrest fires iff
    # PaO2 is still below 15 at exactly that mtime moment - matching FDA
    # patientWrapper.R's one-check rule.
    pao2_below <- (p_a_o2 < 15) * 1.0
    mtime(t_arrest_check) <- first_cross_time + 3.6667
    ca_arrest_decision <- (t >= t_arrest_check) * (crossed_latch > 0.5) *
                          (check_complete < 0.5) * pao2_below

    # ODEs
    # FDA delaymymod.c line 543: Venti_engine = totalVent * Bmax / 60,
    # Bmax = 0.66. After per-second to per-minute unit conversion the
    # alveolar gas-exchange ODEs see venti * Bmax.
    venti_alveolar <- venti * 0.66
    d/dt(palv_co2) <- (venti_alveolar * (0   - palv_co2) +
                      863 * q_total * (1 - 0.024) * (c_v_co2 - c_e_co2)) / 3.28
    d/dt(palv_o2)  <- (venti_alveolar * (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 * m_o2_safety_brain) / 1.32
    d/dt(cb_o2)  <- (qb * (c_a_o2  - c_vb_o2 ) + m_b_o2_real  * m_o2_safety_brain) / 1.32
    d/dt(ct_co2) <- (qt * (c_a_co2 - c_vt_co2) + m_t_co2_real * m_o2_safety_tissue) / 38.68
    d/dt(ct_o2)  <- (qt * (c_a_o2  - c_vt_o2 ) + m_t_o2_real  * m_o2_safety_tissue) / 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   * (clag_p_b_co2 - 45.27)) / 1
    d/dt(alpha_h) <- (hstat - alpha_h) / 5
    # First-crossing detector + capture + decision-window closer + commit latch.
    # ca_arrest_latch grows fast (1000) so it reaches 1 well before
    # check_complete (rate 10, ~7 s window) closes the decision.
    d/dt(crossed_latch)    <- (1 - crossed_latch) * pao2_below * 100
    d/dt(first_cross_time) <- (crossed_latch < 0.5) * 1.0
    d/dt(check_complete)   <- (1 - check_complete) * (t >= t_arrest_check) *
                              (crossed_latch > 0.5) * 10
    d/dt(ca_arrest_latch)  <- (1 - ca_arrest_latch) * ca_arrest_decision * 1000
    # FDA delaypars.R Cim25 = 0.004 * 6 = 0.024 /s; * 60 s/min = 1.44 /min
    d/dt(im_arrest)        <- -ca_arrest_latch * (0.004 * 6 * 60) * im_arrest

    # Initial conditions (FDA delaystates.R nominal values). These are
    # NOT internally self-consistent for the chemoreflex feedback - see
    # the Mann2022Equilibrate() step above and below. Subsequent
    # rxode2::rxSolve calls pass inits = eq_states to override these with
    # the pre-equilibrium values.
    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
    crossed_latch(0)    <- 0
    first_cross_time(0) <- 0
    check_complete(0)   <- 0
    ca_arrest_latch(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_arrest_latch
    Tcross  <- first_cross_time
    Tcheck  <- t_arrest_check
  })
}

# 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
  )
}

Pre-dose equilibration via Mann2022Equilibrate()

The FDA delaystates.R nominal initial conditions for the Mann 2022 physiology layer (palv_co2 = 40.28, cb_co2 = 0.645 → p_b_co2 ≈ 93.85 mm Hg, while the chemoreflex baseline reference p_b_co2_0 = 45.27) are NOT internally self-consistent: at t = 0 the central chemoreflex sees a 48.6 mm Hg error in p_b_co2 − p_b_co2_0 and transiently drives ventilation to ~ 19 L/min before the system settles. If the overdose simulation starts from these nominal values, the initial-condition transient leaves the system in a state that subsequent opioid doses cannot push into sustained PaO2 < 15 mm Hg — and the cardiac-arrest endpoint then severely under-counts.

The FDA simulateToGetOD_IM.R reference script avoids this by running a no-drug pre-simulation per subject (call to fundede) and using the FINAL state as the dose-time initial state. nlmixr2lib exposes the same pattern as the Mann2022Equilibrate() helper. Call it once per parameter set, then pass the resulting equilibrium states to rxode2::rxSolve() as the inits argument of every subsequent dose-bearing call:

fent_pars_typical <- chain_params(opioid = "fentanyl",
                                  antagonist = "none",
                                  patient_type = 1L)
fent_pars_typical <- fent_pars_typical[setdiff(names(fent_pars_typical),
                                               "patient_type")]
eq_states <- Mann2022Equilibrate(mod_fent, params = fent_pars_typical)
str(eq_states)
#> List of 11
#>  $ palv_co2: num 38.7
#>  $ palv_o2 : num 101
#>  $ cb_co2  : num 0.613
#>  $ cb_o2   : num 0.001
#>  $ ct_co2  : num 0.611
#>  $ ct_o2   : num 0.00123
#>  $ yco2    : num -0.0934
#>  $ yo2     : num 0.00067
#>  $ dp_state: num -0.2
#>  $ dc_state: num 0.344
#>  $ alpha_h : num 0.878

The equilibrium does NOT depend on opioid PK parameters (no drug is present during the pre-equilibration), so the same eq_states are valid for every subject in the IIV cohort and for both fentanyl- and carfentanil-flavoured chains. Every rxSolve() call in Steps 1–4 below passes inits = eq_states to override the FDA nominal default initial conditions baked into the model with the equilibrium values.

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.

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,
                                       inits  = eq_states,
                                       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.")
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 0.35 1.35
fentanyl 1.630 mg, naive 0.25 1.25
fentanyl 2.970 mg, naive 0.20 1.20
carfentanil 0.010 mg, naive 2.50 3.50
carfentanil 0.012 mg, naive 1.85 2.85
carfentanil 0.022 mg, naive 0.40 1.40
fentanyl 0.820 mg, chronic 0.45 1.45
fentanyl 1.630 mg, chronic 0.30 1.30
fentanyl 2.970 mg, chronic 0.25 1.25
carfentanil 0.010 mg, chronic 3.15 4.15
carfentanil 0.012 mg, chronic 2.20 3.20
carfentanil 0.022 mg, chronic 0.50 1.50

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 Mann 2022 IIV
# omega^2 values (Table S1 for the PK row, Table S2 bootstrap for binding).
# Broader IIV than the prior 3-parameter sketch - matching the published
# variability is what lets the small-N grid land near the published rates
# without a 2000-subject bootstrap.
draw_subjects <- function(n) {
  tibble(
    id       = seq_len(n),
    # PK IIV: published Mann 2022 Table S1 omega^2 values (Algera 2021
    # popPK fit). Bootstrap-vs-omega expansion did not close the residual
    # fentanyl arrest gap, so we stay at the published omegas.
    cl_mult  = exp(rnorm(n, 0, sqrt(0.087))),   # Mann Table S1 omega2 CL = 0.087 (30% CV)
    vc_mult  = exp(rnorm(n, 0, sqrt(0.43))),    # Mann Table S1 omega2 V1 = 0.43  (73% CV)
    vp_mult  = exp(rnorm(n, 0, sqrt(0.52))),    # Mann Table S1 omega2 V2 = 0.52  (82% CV)
    vp2_mult = exp(rnorm(n, 0, sqrt(0.060))),   # Mann Table S1 omega2 V3 = 0.060 (25% CV)
    q_mult   = exp(rnorm(n, 0, sqrt(0.39))),    # Mann Table S1 omega2 Q2 = 0.39  (69% CV)
    q2_mult  = exp(rnorm(n, 0, sqrt(0.098))),   # Mann Table S1 omega2 Q3 = 0.098 (32% CV)
    # Binding IIV calibrated against the FDA reference bootstrap
    # boot_pars_fentanyl.csv (n = 2000) summary statistics:
    #   A1 (Kon)  median 3.1e-5, max 1.9e-4 -> sigma ~ 0.74 (~ 100% CV)
    #   B1 (Koff) median 4.3e-3, IQR 4.1-4.6e-3 -> sigma ~ 0.06 (~ 6% CV)
    # Kon dominates how deeply receptor occupancy can spike during an
    # overdose; Koff is much tighter empirically than the Mann 2022
    # Table S2 standard-error figures suggested.
    kon_mult  = exp(rnorm(n, 0, 0.74)),
    koff_mult = exp(rnorm(n, 0, 0.06))
  )
}

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["vp"]       <- base_pars[["vp"]]       * s$vp_mult
    sub_pars["vp2"]      <- base_pars[["vp2"]]      * s$vp2_mult
    sub_pars["q"]        <- base_pars[["q"]]        * s$q_mult
    sub_pars["q2"]       <- base_pars[["q2"]]       * s$q2_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["koff_op_s"] <- base_pars[["koff_op_s"]] * s$koff_mult
    sub_pars["koff_an_s"] <- base_pars[["koff_an_s"]] * s$koff_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,
                                           inits  = eq_states,
                                           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%. The fentanyl no-rescue arms reproduce the published rates closely: low / mid / high dose = 23.3 / 53.3 / 80.0% vs published 18.1 / ~55 / 90.2%. Carfentanil no-rescue arms low / mid / high = 50.0 / 53.3 / 90.0% land inside the published range. Rank order across rescue arms (no antagonist > IN naloxone > IN nalmefene) reproduces, and the dose-response direction within each opioid is preserved. The vignette captures every FDA delaymymod.c / simulateToGetOD_IM.R structural element: cardiac-arrest one-shot trigger (`patientWrapper.R` schedule at first_cross + 220 s) via rxode2 state-dependent `mtime`; 7 s peripheral / 11 s central chemoreflex transport delays via the `lagState` C user function; Q_Scale shock-state PK amplification (`delaymymod.c` lines 358-368) via the Q_TOTAL_LPM covariate input on the Mann 2022 PK layer; pre-dose equilibration (FDA `fundede` call) via `Mann2022Equilibrate()`; low-PaO2 metabolism safety clamps (`delaymymod.c` lines 296-306); and the Bmax = 0.66 dead-space fraction on alveolar ventilation (`delaymymod.c` line 543). See `Mann_2022_respiratory_physiology` deviation notes.")
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%. The fentanyl no-rescue arms reproduce the published rates closely: low / mid / high dose = 23.3 / 53.3 / 80.0% vs published 18.1 / ~55 / 90.2%. Carfentanil no-rescue arms low / mid / high = 50.0 / 53.3 / 90.0% land inside the published range. Rank order across rescue arms (no antagonist > IN naloxone > IN nalmefene) reproduces, and the dose-response direction within each opioid is preserved. The vignette captures every FDA delaymymod.c / simulateToGetOD_IM.R structural element: cardiac-arrest one-shot trigger (patientWrapper.R schedule at first_cross + 220 s) via rxode2 state-dependent mtime; 7 s peripheral / 11 s central chemoreflex transport delays via the lagState C user function; Q_Scale shock-state PK amplification (delaymymod.c lines 358-368) via the Q_TOTAL_LPM covariate input on the Mann 2022 PK layer; pre-dose equilibration (FDA fundede call) via Mann2022Equilibrate(); low-PaO2 metabolism safety clamps (delaymymod.c lines 296-306); and the Bmax = 0.66 dead-space fraction on alveolar ventilation (delaymymod.c line 543). See Mann_2022_respiratory_physiology deviation notes.
opioid dose_mg IN nalmefene 3 mg HCl IN naloxone 4 mg HCl No antagonist
carfentanil 0.0 16.7 43.3 50.0
carfentanil 0.0 30.0 43.3 36.7
carfentanil 0.0 26.7 80.0 70.0
fentanyl 0.8 0.0 0.0 23.3
fentanyl 1.6 0.0 10.0 46.7
fentanyl 3.0 0.0 46.7 80.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.

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,
                                       inits  = eq_states,
                                       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).")
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 3.1 24.4
carfentanil 0.02 4.4 27.4 27.4
fentanyl 1.63 0.1 3.3 27.5
fentanyl 2.97 3.1 28.1 28.1 

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

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.R ini() 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.R and 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_* and Laffont_2024_* model files remain the authoritative source for the parameter values; this vignette composes them into a runnable end-to-end system.