tests/testthat/test-event-sensitivities.R

In rxode2: Facilities for Simulating from ODE-Based Models

rxTest({
  # Event ("jump") sensitivities -- Phase A: build-time index map and the
  # symbolic total-derivatives of the dosing parameters (alag, F).
  # See ~/src/rxode2-event-sensitivities-plan.md.

  .mkMod <- function(modStr) {
    # depot+central first-order absorption with modeled alag + F, sensitivities
    # wrt two etas, full Jacobian (as eventSens="jump" will force).
    rxode2(modStr, calcSens = c("eta_ka", "eta_lag"), calcJac = TRUE)
  }

  .modConstF <- "
    ka <- exp(tka + eta_ka)
    cl <- exp(tcl)
    v  <- exp(tv)
    alag(depot) <- exp(tlag + eta_lag)
    f(depot)    <- expit(tf)
    d/dt(depot)   <- -ka * depot
    d/dt(central) <-  ka * depot - cl / v * central
  "

  .modStateF <- "
    ka <- exp(tka + eta_ka)
    cl <- exp(tcl)
    v  <- exp(tv)
    alag(depot) <- exp(tlag + eta_lag)
    f(depot)    <- 1 / (1 + exp(-(tf + 0.1 * central)))
    d/dt(depot)   <- -ka * depot
    d/dt(central) <-  ka * depot - cl / v * central
  "

  test_that(".rxEventSensMap relates sens compartments to states/params", {
    im <- .rxEventSensMap(.mkMod(.modConstF))
    expect_equal(im$states, c("depot", "central"))
    expect_equal(im$nState, 2L)
    expect_equal(unname(im$stateCmt[c("depot", "central")]), c(1L, 2L))
    expect_setequal(im$sensParams, c("eta_ka", "eta_lag"))
    # one map row per (state, param): 2 states x 2 params
    expect_equal(nrow(im$map), 4L)
    # rx__sens_depot_BY_eta_ka__ must be a real compartment mapped to depot
    .row <- im$map[im$map$state == "depot" & im$map$param == "eta_ka", ]
    expect_equal(.row$stateCmt, 1L)
    expect_true(.row$sensCmt > im$nState) # sens compartments come after states
    # depot carries modeled alag + F, not rate/dur
    expect_equal(im$lagCmt, 1L)
    expect_equal(im$fCmt, 1L)
    expect_length(im$rateCmt, 0L)
    expect_length(im$durCmt, 0L)
  })

  test_that(".rxEventSensMap returns NULL without sensitivities", {
    m <- rxode2("ka <- 1\nd/dt(depot) <- -ka * depot")
    expect_null(.rxEventSensMap(m))
  })

  test_that(".rxEventSensSplit handles underscored state/param names", {
    s <- c("rx__sens_depot_BY_eta_ka__", "rx__sens_central_BY_eta_lag__")
    sp <- .rxEventSensSplit(s, c("depot", "central"))
    expect_equal(sp$state, c("depot", "central"))
    expect_equal(sp$param, c("eta_ka", "eta_lag"))
    # second-order names are not first-order: state/param NA
    sp2 <- .rxEventSensSplit("rx__sens_depot_BY_eta_ka_BY_eta_lag__",
                             c("depot", "central"))
    expect_true(is.na(sp2$state))
  })

  test_that("alag total-derivative is the plain partial when state-independent", {
    d <- .rxEventSensDerivs(.mkMod(.modConstF))
    # alag depends only on eta_lag (no states) -> single non-zero row, no coupling
    expect_equal(nrow(d$lag), 1L)
    expect_equal(d$lag$param, "eta_lag")
    expect_false(any(grepl("rx__sens_", d$lag$expr)))
    # d(exp(eta_lag+tlag))/d(eta_lag) = exp(eta_lag+tlag)
    expect_equal(eval(parse(text = d$lag$expr),
                      list(eta_lag = 0.3, tlag = log(2))),
                 exp(0.3 + log(2)))
  })

  test_that("F total-derivative picks up the state-coupling term", {
    d <- .rxEventSensDerivs(.mkMod(.modStateF))
    # F depends on `central` -> d(F)/d(eta) is purely the coupling term
    # (no direct eta in F), nonzero for BOTH etas via S_central.
    expect_setequal(d$f$param, c("eta_ka", "eta_lag"))
    expect_true(all(grepl("rx__sens_central_BY_", d$f$expr)))
  })

  test_that(".rxEventSensMode resolves and validates the mode", {
    expect_equal(.rxEventSensMode("jump"), "jump")
    expect_equal(.rxEventSensMode("both"), "both")
    expect_equal(.rxEventSensMode("fd"), "fd")
    expect_equal(.rxEventSensMode("fdAll"), "fdAll")
    expect_error(.rxEventSensMode("bogus"), "eventSens")
    withr::with_options(list(rxode2.eventSens = "jump"), {
      expect_equal(.rxEventSensMode(NULL), "jump")
    })
    withr::with_options(list(rxode2.eventSens = NULL), {
      expect_equal(.rxEventSensMode(NULL), "jump")
    })
  })

  test_that("eventSens='jump' (default) is a no-op on rxode2()", {
    m <- rxode2("d/dt(depot) <- -ka*depot\nka <- 1",
                calcSens = "ka")
    expect_equal(m$eventSens, "jump")
    expect_null(m$eventSensInfo)
  })

  test_that("eventSens='jump' forces calcJac and attaches Phase-A info", {
    m <- rxode2(.modConstF, calcSens = c("eta_ka", "eta_lag"),
                eventSens = "jump")
    expect_equal(m$eventSens, "jump")
    expect_true(isTRUE(m$calcJac))
    info <- m$eventSensInfo
    expect_equal(info$mode, "jump")
    expect_equal(nrow(info$map$map), 4L)
    # constant F (expit(tf)) -> no eta-dependent F derivative; alag depends on
    # eta_lag -> exactly one lag derivative row.
    expect_equal(nrow(info$derivs$lag), 1L)
    expect_equal(nrow(info$derivs$f), 0L)
  })

  test_that("eventSens='jump' on a model without sensitivities is a no-op", {
    m <- rxode2("d/dt(depot) <- -depot", eventSens = "jump")
    expect_null(m$eventSensInfo)
  })

  test_that("eventSens on a linCmt()+mtime() model downgrades to fd", {
    m <- rxode2({
      C2 <- linCmt(CL, V)
      mtime(t1) <- mt1
      mtime(t2) <- mt2
    }, eventSens = "jump", linCmtSens = "linCmtA", linCmtSensType = "A")
    # linCmt models downgrade to finite differences for event sensitivities (#1119)
    expect_equal(m$eventSens, "fd")
    expect_null(m$eventSensInfo)
    e <- eventTable() |>
      add.dosing(dose = 3, nbr.doses = 2, dosing.interval = 8) |>
      add.sampling(0:16)
    s <- rxSolve(m, e, params = c(V = 20, CL = 25, mt1 = 0.5, mt2 = 1.75))
    expect_true(all(c(0.5, 1.75) %in% s$time))
  })

  test_that("mixed ODE+linCmt downgrades to fd for event sensitivities", {
    m <- rxode2({
      C2 <- linCmt(CL, V)
      alag(eff) <- exp(tlag + eta_lag)
      d/dt(eff) <- kin - kout * eff
    }, calcSens = "eta_lag", eventSens = "fd",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    # any linCmt model downgrades to FD (#1119): no jump metadata
    expect_equal(m$eventSens, "fd")
    expect_null(m$eventSensInfo)
  })

  test_that("eventSens='fdAll' keeps full finite-difference fallback", {
    m <- rxode2({
      C2 <- linCmt(CL, V)
      d/dt(eff) <- kin - kout * eff
    }, eventSens = "fdAll", linCmtSens = "linCmtA", linCmtSensType = "A")
    expect_equal(m$eventSens, "fd")
    expect_null(m$eventSensInfo)
  })

  test_that("eventSens='fdAll' works correctly in a population (multi-subject) solve", {
    # Regression for a previously-untested combination (plan Section 0.4 gap):
    # confirm fdAll's per-subject solving (each subject independently, no
    # jump machinery/eventSensInfo at all) isn't corrupted by the population
    # (multi-subject/parallel) solve path -- each subject's trajectory must
    # match an INDEPENDENT single-subject solve with that subject's own
    # modeled-lag parameter.
    mLin <- rxode2({
      C2 <- linCmt(CL, V)
      alag(central) <- exp(tlag + eta_lag)
    }, calcSens = "eta_lag", eventSens = "fdAll",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    expect_equal(mLin$eventSens, "fd")
    expect_null(mLin$eventSensInfo)

    th <- c(V = 20, CL = 25, tlag = 0)
    ev <- data.frame(eta_lag = c(0, 0.3, -0.2))
    e <- et(amt = 100, cmt = "central", id = 1:3) |> et(seq(0, 12, 0.5), id = 1:3)

    sPop <- rxSolve(mLin, e, th, iCov = ev)
    for (.id in 1:3) {
      e1 <- eventTable() |>
        add.dosing(dose = 100, cmt = "central") |>
        add.sampling(seq(0, 12, 0.5))
      s1 <- rxSolve(mLin, e1, c(th, eta_lag = ev$eta_lag[.id]))
      expect_equal(sPop$central[sPop$id == .id], s1$central,
        tolerance = 1e-8, info = paste0("id: ", .id)
      )
    }
  })

  test_that("pure linCmt with modeled alag downgrades to fd", {
    m <- rxode2({
      C2 <- linCmt(CL, V)
      alag(central) <- exp(tlag + eta_lag)
    }, calcSens = "eta_lag", eventSens = "fd",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    # pure linCmt downgrades to FD (#1119)
    expect_equal(m$eventSens, "fd")
    expect_null(m$eventSensInfo)
  })

  test_that("pure linCmt hybrid event path matches the ODE equivalent", {
    mLin <- rxode2({
      C2 <- linCmt(CL, V)
      alag(central) <- exp(tlag + eta_lag)
    }, calcSens = "eta_lag", eventSens = "fd",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    mOde <- rxode2({
      d/dt(central) <- -CL / V * central
      alag(central) <- exp(tlag + eta_lag)
    }, calcSens = "eta_lag", eventSens = "jump")

    e <- eventTable() |>
      add.dosing(dose = 100, cmt = "central") |>
      add.sampling(seq(0, 12, 0.5))
    p <- c(V = 20, CL = 25, tlag = 0, eta_lag = 0)

    sLin <- rxSolve(mLin, e, params = p)
    sOde <- rxSolve(mOde, e, params = p)

    expect_equal(sLin$central, sOde$central, tolerance = 1e-5)
  })

  test_that("mixed ODE+linCmt with explicit jump still downgrades to fd", {
    # Even eventSens="jump" downgrades to FD when linCmt() is present (#1119):
    # the linCmt-compartment moving-boundary jump is not implemented, so all
    # linCmt models use finite differences for event sensitivities.
    m <- rxode2({
      C2 <- linCmt(CL, V)
      alag(eff) <- exp(tlag + eta_lag)
      d/dt(eff) <- kin - kout * eff
    }, calcSens = "eta_lag", eventSens = "jump",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    expect_equal(m$eventSens, "fd")
    expect_null(m$eventSensInfo)
  })

  test_that("ODE d/dt() colliding with a linCmt reserved compartment name warns", {
    # Naming an ODE compartment `depot` alongside an oral/multi-cmt linCmt (which
    # reserves `depot`) conflates the two compartments: the ODE state loses its
    # sensitivity expansion, so its sensitivities are silently incorrect.  Warn.
    expect_warning(
      rxode2({
        ka <- exp(tka)
        alag(depot) <- 2 * exp(eta_lag)
        d/dt(depot)   <- -ka * depot
        d/dt(peri)    <-  ka * depot - 0.3 * peri
        C2 <- linCmt(CL, V)
      }, calcSens = "eta_lag", eventSens = "jump",
      linCmtSens = "linCmtA", linCmtSensType = "A"),
      "share a name with linCmt"
    )
    # a valid (non-colliding) linCmt model: no warning; downgrades to fd (#1119)
    mOk <- suppressWarnings(rxode2({
      ka <- exp(tka)
      alag(gut) <- 2 * exp(eta_lag)
      d/dt(gut) <- -ka * gut
      d/dt(eff) <-  ka * gut - 0.3 * eff
      C2 <- linCmt(CL, V)
    }, calcSens = "eta_lag", eventSens = "fd",
    linCmtSens = "linCmtA", linCmtSensType = "A"))
    expect_equal(mOk$eventSens, "fd")
    expect_null(mOk$eventSensInfo)
  })

  test_that("mixed ODE+linCmt hybrid event path matches the ODE equivalent", {
    mLin <- rxode2({
      C2 <- linCmt(CL, V)
      alag(central) <- exp(tlag + eta_lag)
      d/dt(eff) <- kin - kout * eff
    }, calcSens = "eta_lag", eventSens = "fd",
    linCmtSens = "linCmtA", linCmtSensType = "A")
    mOde <- rxode2({
      d/dt(central) <- -CL / V * central
      alag(central) <- exp(tlag + eta_lag)
      d/dt(eff) <- kin - kout * eff
    }, calcSens = "eta_lag", eventSens = "jump")

    e <- eventTable() |>
      add.dosing(dose = 100, cmt = "central") |>
      add.sampling(seq(0, 12, 0.5))
    p <- c(V = 20, CL = 25, kin = 1, kout = 0.2, tlag = 0, eta_lag = 0)

    sLin <- rxSolve(mLin, e, params = p)
    sOde <- rxSolve(mOde, e, params = p)

    expect_equal(sLin$central, sOde$central, tolerance = 1e-5)
    expect_equal(sLin$eff, sOde$eff, tolerance = 1e-5)
  })

  test_that("eventSens='jump' works with ODE+mtime() models", {
    .mod <- "
      d/dt(depot) <- -ka*depot
      d/dt(central) <- ka*depot-kel*central
      ka <- exp(tka)
      kel <- exp(tkel)
      mtime(t1) <- mt1
      mtime(t2) <- mt2
    "
    mj <- rxode2(.mod, calcSens = "tka", eventSens = "jump")
    mfd <- rxode2(.mod, calcSens = "tka", eventSens = "fd")
    e <- eventTable() |>
      add.dosing(dose = 100, cmt = "depot") |>
      add.sampling(0:12)
    .p <- c(tka = 0, tkel = log(0.2), mt1 = 0.5, mt2 = 1.75)
    sj <- rxSolve(mj, e, params = .p)
    sfd <- rxSolve(mfd, e, params = .p)
    expect_true(isTRUE(mj$calcJac))
    expect_equal(sj$central, sfd$central, tolerance = 1e-5)
  })

  test_that(".rxEventSensCLines emits correctly-indexed C assignment lines", {
    m <- rxode2(.modStateF, calcSens = c("eta_ka", "eta_lag"),
                eventSens = "jump")
    cl <- .rxEventSensCLines(m$eventSensInfo)
    expect_equal(cl$nSensParam, 2L)
    expect_equal(unname(cl$paramIdx[c("eta_ka", "eta_lag")]), c(0L, 1L))
    # alag(depot): cmt0=0, eta_lag idx 1 -> _dLagSave[0*2+1] = _dLagSave[1]
    expect_length(cl$lag, 1L)
    expect_match(cl$lag, "_dLagSave\\[1\\] = exp\\(eta_lag\\+tlag\\);", fixed = FALSE)
    # state-dependent F: both etas via S_central, indices 0 and 1
    expect_length(cl$f, 2L)
    expect_true(any(grepl("_dFSave[0]", cl$f, fixed = TRUE)))
    expect_true(any(grepl("_dFSave[1]", cl$f, fixed = TRUE)))
    # right-hand sides reference the exact local names the F function declares
    expect_true(all(grepl("rx__sens_central_BY_", cl$f)))
  })

  test_that(".rxEventSensCLines is NULL for fd / no-sens", {
    expect_null(.rxEventSensCLines(NULL))
  })

  test_that(".rxEventSensCExpr maps indexed vs plainly-declared params (1st + 2nd order)", {
    # nlmixr2's indexed THETA[n]/ETA[n] map to the codegen locals _THETA_n_/_ETA_n_
    # (populated from _PP[]); this is what the dF/d2F assignment lines reference.
    expect_equal(.rxEventSensCExpr("exp(THETA[2])"), "exp(_THETA_2_)")
    expect_equal(.rxEventSensCExpr("THETA[1]*ETA[3]"), "_THETA_1_*_ETA_3_")
    # But a model may declare the parameter under its own plain name THETA_n_/ETA_n_
    # (e.g. nlmixr2est's augmented outer-gradient model, whose direction params are
    # named THETA_n_/ETA_n_); the preamble then declares them verbatim, so the dose
    # (and 2nd-order dose) derivatives must reference the plain name -- otherwise the
    # emitted _THETA_n_ is undeclared and the model fails to compile.
    expect_equal(.rxEventSensCExpr("exp(THETA[2])", plainParams = "THETA_2_"), "exp(THETA_2_)")
    expect_equal(.rxEventSensCExpr("THETA[1]*ETA[3]", plainParams = c("THETA_1_", "ETA_3_")),
                 "THETA_1_*ETA_3_")
    # mixed: only the plainly-declared params drop the leading underscore
    expect_equal(.rxEventSensCExpr("THETA[1]*ETA[3]", plainParams = "THETA_1_"), "THETA_1_*_ETA_3_")
    # info$params carries the declared names so .rxEventSensCLines picks the right form
    m <- rxode2(.modStateF, calcSens = c("eta_ka", "eta_lag"), eventSens = "jump")
    expect_true("params" %in% names(m$eventSensInfo))
  })

  test_that("eventSens='jump' model with state-dependent F compiles and solves", {
    # Exercises the full lines channel: R-generated dLag/dF assignment lines are
    # pushed to codegen, spliced into the function bodies (with the state-coupling
    # term referencing populated sens locals), and the model compiles + solves.
    # (The exact emitted C is asserted at the R level by the .rxEventSensCLines
    # tests above; here we verify the end-to-end build/solve.)
    m <- rxode2(.modStateF, calcSens = c("eta_ka", "eta_lag"),
                eventSens = "jump")
    expect_s3_class(m, "rxode2")
    e <- et(amt = 100, cmt = "depot")
    e <- et(e, seq(0, 12, 4))
    ini <- c(tka = 0, tcl = 1, tv = 2, tlag = 0, tf = 0,
             eta_ka = 0, eta_lag = 0)
    s <- rxSolve(m, e, ini)
    expect_true(nrow(s) > 0L)
    # jump-mode solve matches fd-mode for the physical states: the jump injection
    # only touches the sensitivity compartments, never the physical state.  The
    # match is to solver tolerance (not bit-identical) because jump forces
    # calcJac=TRUE -> the analytic Jacobian gives LSODA a slightly different step
    # path than fd's numerical Jacobian (~1e-7 on the trajectory).
    mfd <- rxode2(.modStateF, calcSens = c("eta_ka", "eta_lag"),
                  eventSens = "fd")
    sfd <- rxSolve(mfd, e, ini)
    expect_equal(s$central, sfd$central, tolerance = 1e-5)
  })

  test_that("additive-bolus F jump (ddelta row) matches finite differences", {
    # Constant alag (no eta) + F = expit(tf + eta_f): the only event-sensitivity
    # contribution at the dose is the ddelta row, d(delta)/d(eta_f) = amt*dF.
    # The analytic jump-mode sensitivity must match a central difference of the
    # central state wrt eta_f; the fd-mode sens ODE alone (no jump) must NOT.
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- 1
      f(depot)    <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0, tcl = 1, tv = 2, tf = 0.3, eta_ka = 0, eta_f = 0)
    e <- et(amt = 100, cmt = "depot")
    e <- et(e, seq(0, 12, 2))
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = c("eta_ka", "eta_f"), eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")
    analytic <- sj[["rx__sens_central_BY_eta_f__"]]
    h <- 1e-4
    pp <- pars; pp["eta_f"] <- pars["eta_f"] + h
    pm <- pars; pm["eta_f"] <- pars["eta_f"] - h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    # analytic jump matches FD to FD precision
    expect_equal(analytic, fd, tolerance = 1e-4)
    # fd-mode sens ODE alone is missing the dF contribution -> clearly different
    sfd <- .central(pars, "fd")[["rx__sens_central_BY_eta_f__"]]
    expect_gt(max(abs(sfd - fd)), 1)
  })

  test_that("additive-bolus lag jump (dtau row) matches finite differences", {
    # alag = exp(tlag + eta_lag) depends on eta_lag, F constant: the jump is the
    # dtau row, -J[k][c]*delta*d(alag)/d(eta_lag), with J taken from a central
    # difference of dydt.  Observe AWAY from the (lagged) dose time, since at the
    # exact event time the sensitivity is discontinuous and a central FD straddles
    # the jump (the analytic value is the correct right-limit).
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- exp(tlag + eta_lag)
      f(depot)    <- expit(tf)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tf = 0.5, tlag = 0,
              eta_ka = 0, eta_lag = 0)
    e <- et(amt = 100, cmt = "depot")
    e <- et(e, seq(0.5, 12, 1)) # off the dose time (alag = 1)
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = c("eta_ka", "eta_lag"), eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")
    h <- 1e-4
    pp <- pars; pp["eta_lag"] <- h
    pm <- pars; pm["eta_lag"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    expect_equal(sj[["rx__sens_central_BY_eta_lag__"]], fd, tolerance = 1e-4)
    # eta_ka has no lag/F dependence -> no jump; must still match FD
    ppk <- pars; ppk["eta_ka"] <- h
    pmk <- pars; pmk["eta_ka"] <- -h
    fdK <- (.central(ppk, "fd")$central - .central(pmk, "fd")$central) / (2 * h)
    expect_equal(sj[["rx__sens_central_BY_eta_ka__"]], fdK, tolerance = 1e-4)
  })

  test_that("additive-bolus lag jump: coincident-time output uses one-sided limits", {
    # Regression for the doubled-sensitivity bug when an observation time exactly
    # equals a modeled-lag dose arrival.  alag(depot) = 2 * exp(eta_lag): with
    # eta_lag = 0 the dose (t = 0) arrives at EXACTLY t = 2, where an observation
    # is placed.  At that coincident time the sensitivity has a genuine jump
    # discontinuity; the reported value must match the STATE's reported side:
    #   - central (not dosed, continuous): reported state is the pre-arrival
    #     value, so its sensitivity must be the pre-jump (right/forward) limit.
    #   - depot (dosed): reported state is post-dose, so its sensitivity keeps the
    #     post-jump (left/backward) limit.
    # Before the fix the coincident central sensitivity was the post-jump
    # (backward) limit -- the wrong side, ~2x a symmetric FD -- which inflated the
    # FOCEi objective on models whose observations land on lagged arrivals.
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- 2 * exp(eta_lag)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = log(1.15), tcl = log(0.135), tv = log(8), eta_lag = 0)
    e <- et(amt = 100, cmt = "depot") |> et(c(1.9, 2.0, 2.1, 4))
    m <- rxode2(.mod, calcSens = "eta_lag", eventSens = "jump")
    sj <- rxSolve(m, e, pars)
    h <- 1e-5
    pp <- pars; pp["eta_lag"] <-  h
    pm <- pars; pm["eta_lag"] <- -h
    s0 <- rxSolve(m, e, pars)
    sP <- rxSolve(m, e, pp)
    sM <- rxSolve(m, e, pm)
    .at <- function(s, tm, col) s[[col]][which(abs(s$time - tm) < 1e-8)[1]]
    # one-sided finite differences of the STATE at the coincident time t = 2
    fwdCentral <- (.at(sP, 2, "central") - .at(s0, 2, "central")) / h   # right
    bwdCentral <- (.at(s0, 2, "central") - .at(sM, 2, "central")) / h   # left
    bwdDepot   <- (.at(s0, 2, "depot")   - .at(sM, 2, "depot"))   / h
    ajCentral <- .at(sj, 2, "rx__sens_central_BY_eta_lag__")
    ajDepot   <- .at(sj, 2, "rx__sens_depot_BY_eta_lag__")
    # central: analytic == pre-jump (forward/right) limit (~0 here), and clearly
    # NOT the old doubled (backward/left) value.
    expect_equal(ajCentral, fwdCentral, tolerance = 1e-3)
    expect_true(abs(ajCentral - bwdCentral) > 0.1 * abs(bwdCentral))
    # depot: dosed compartment keeps the post-jump (backward/left) limit.
    expect_equal(ajDepot, bwdDepot, tolerance = 1e-2)
    # off-dose observations still match a two-sided central difference.
    for (tm in c(1.9, 2.1, 4)) {
      fd <- (.at(sP, tm, "central") - .at(sM, tm, "central")) / (2 * h)
      expect_equal(.at(sj, tm, "rx__sens_central_BY_eta_lag__"), fd,
                   tolerance = 1e-3)
    }
  })

  test_that("additive-bolus jumps accumulate across a multi-dose regimen", {
    # alag(depot)=exp(tlag+eta_lag), F=expit(tf+eta_f): every dose contributes
    # both a dtau and a ddelta jump.  The rows use method "add", so the jumps
    # must accumulate across the regimen.  A multi-dose schedule must still match
    # a central difference of the solution wrt each eta.
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- exp(tlag + eta_lag)
      f(depot)    <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tf = 0.3, tlag = 0,
              eta_ka = 0, eta_lag = 0, eta_f = 0)
    e <- et(amt = 100, cmt = "depot", ii = 6, addl = 3) |>
      et(seq(0.5, 30, 1.5))
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = c("eta_lag", "eta_f"), eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")
    h <- 1e-4
    for (.eta in c("eta_lag", "eta_f")) {
      pp <- pars; pp[.eta] <- h
      pm <- pars; pm[.eta] <- -h
      fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
      expect_equal(sj[[paste0("rx__sens_central_BY_", .eta, "__")]], fd,
                   tolerance = 1e-3)
    }
  })

  test_that("replacement event jump (dp_j -> 0) matches finite differences", {
    # Replace central with a constant (50) at t=5 via an evid=5 event.  The
    # replaced state's value no longer depends on eta_ka, so d(central)/d(eta_ka)
    # must jump to 0 at the event and rebuild afterwards.  The analytic jump must
    # match a central difference of the solution; the fd-mode sensitivity ODE
    # (which never zeroes the sens compartment) must NOT.
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, eta_ka = 0)
    e <- et(amt = 100, cmt = "depot") |>
      et(time = 5, amt = 50, cmt = "central", evid = 5) |>
      et(seq(0.5, 12, 1))
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = "eta_ka", eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")
    h <- 1e-4
    pp <- pars; pp["eta_ka"] <- h
    pm <- pars; pm["eta_ka"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    # observations at/after the replace must match FD (analytic jump correct)
    .post <- sj$time >= 5
    expect_equal(sj[["rx__sens_central_BY_eta_ka__"]][.post], fd[.post],
                 tolerance = 1e-4)
    # fd-mode sens ODE alone never zeroes the replaced state's sens -> differs
    sfd <- .central(pars, "fd")[["rx__sens_central_BY_eta_ka__"]]
    expect_gt(max(abs(sfd[.post] - fd[.post])), 1)
  })

  test_that("multiplicative event jump (dp_j *= alpha) matches finite differences", {
    # Multiply central by alpha=0.5 at t=5 via an evid=6 event.  d(central)/deta
    # must be scaled by the same 0.5 at the event.  Analytic jump matches FD; the
    # fd-mode sens ODE (no scaling of the sens compartment) does not.
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, eta_ka = 0)
    e <- et(amt = 100, cmt = "depot") |>
      et(time = 5, amt = 0.5, cmt = "central", evid = 6) |>
      et(seq(0.5, 12, 1))
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = "eta_ka", eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")
    h <- 1e-4
    pp <- pars; pp["eta_ka"] <- h
    pm <- pars; pm["eta_ka"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    .post <- sj$time >= 5
    expect_equal(sj[["rx__sens_central_BY_eta_ka__"]][.post], fd[.post],
                 tolerance = 1e-4)
    sfd <- .central(pars, "fd")[["rx__sens_central_BY_eta_ka__"]]
    expect_gt(max(abs(sfd[.post] - fd[.post])), 0.1)
  })

  test_that("replace/multiply zero the 2nd-order (Hessian) compartment too", {
    # Extends Phase F's additive-bolus-only 2nd-order jump (the "dx1/dp_j"
    # row, NOT the dtau row) to replace/multiply: a constant replacement
    # value has an identically zero 2nd derivative wrt any parameter pair
    # (same reasoning as 1st order), and alpha (treated as parameter-fixed
    # for this row, same convention as 1st order) scales the 2nd-order
    # compartment the same way. Validated against a finite difference of the
    # analytic 1st-order sensitivity (the plan's standard 2nd-order
    # validation strategy -- a raw 2nd central difference of the solution is
    # too imprecise to catch a jump-sized correction).
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2)
    m1 <- rxode2(.mod, calcSens = "tka", eventSens = "jump")
    m2 <- rxode2(.mod, calcSens = "tka", calcSens2 = "tka", eventSens = "jump")
    h <- 1e-4
    pp <- pars; pp["tka"] <- pars["tka"] + h
    pm <- pars; pm["tka"] <- pars["tka"] - h

    e_replace <- et(amt = 100, cmt = "depot") |>
      et(time = 5, amt = 50, cmt = "central", evid = 5) |>
      et(seq(0.5, 12, 1))
    fd_replace <- (rxSolve(m1, e_replace, pp)$rx__sens_central_BY_tka__ -
      rxSolve(m1, e_replace, pm)$rx__sens_central_BY_tka__) / (2 * h)
    s_replace <- rxSolve(m2, e_replace, pars)
    expect_equal(s_replace$rx__sens_central_BY_tka_BY_tka__, fd_replace, tolerance = 1e-3)

    e_mult <- et(amt = 100, cmt = "depot") |>
      et(time = 5, amt = 0.5, cmt = "central", evid = 6) |>
      et(seq(0.5, 12, 1))
    fd_mult <- (rxSolve(m1, e_mult, pp)$rx__sens_central_BY_tka__ -
      rxSolve(m1, e_mult, pm)$rx__sens_central_BY_tka__) / (2 * h)
    s_mult <- rxSolve(m2, e_mult, pars)
    expect_equal(s_mult$rx__sens_central_BY_tka_BY_tka__, fd_mult, tolerance = 1e-3)
  })

  test_that("raw event-table replace with a modeled lag gets a dtau row (Phase B, B2)", {
    # Regression (Phase B "B2" gap, plan Section 4): a raw et(..., evid=5)
    # replace record on a compartment with a MODELED alag was found to
    # completely omit the dx1/dtau row (paper Table 1) -- the modeled lag
    # correctly shifts the event's OWN time (confirmed: alag applies to
    # evid=5/6 records exactly like normal doses), but the sensitivity was
    # silently exactly zero even though the true value is clearly nonzero
    # (confirmed via a hand-derived closed form and FD, ~20.2 here). The
    # in-model replace()/multiply() plugins were previously found NOT to need
    # this row (captured-dosing re-emits correctly) -- this gap is specific
    # to raw event-table evid=5/6 records.
    .mod <- "
      cl <- exp(tcl); v <- exp(tv)
      alag(central) <- exp(tlag + eta_lag)
      d/dt(central) <- -cl / v * central
    "
    pars <- c(tcl = 1, tv = 2, tlag = log(1.1), eta_lag = 0)
    e <- et(time = 5, amt = 50, cmt = "central", evid = 5) |> et(seq(0, 10, 0.05))
    mj <- rxode2(.mod, calcSens = "eta_lag", eventSens = "jump")
    mfd <- rxode2(.mod, calcSens = "eta_lag", eventSens = "fd")
    sj <- rxSolve(mj, e, pars)
    h <- 1e-4
    pp <- pars; pp["eta_lag"] <- pars["eta_lag"] + h
    pm <- pars; pm["eta_lag"] <- pars["eta_lag"] - h
    # tight atol/rtol on the FD reference solves: the default LSODA tolerance
    # is loose enough, relative to h, that its own step-to-step numerical
    # noise dominates a naive central difference here (~1% spurious offset
    # observed with default tolerances; confirmed to vanish under atol/rtol
    # 1e-12 -- same class of pitfall as too-small an h, not a real signal).
    fd <- (rxSolve(mfd, e, pp, atol = 1e-12, rtol = 1e-12)$central -
             rxSolve(mfd, e, pm, atol = 1e-12, rtol = 1e-12)$central) / (2 * h)
    # skip the sample landing exactly on the lagged event time (t=6.1): a
    # central FD straddling a true discontinuity is a known artifact, not a
    # correctness signal (documented in the plan's infusion-jump validation
    # note; the same subtlety applies here).
    .post <- sj$time >= 6.15 & sj$time < 10
    expect_equal(sj[["rx__sens_central_BY_eta_lag__"]][.post], fd[.post], tolerance = 1e-3)
    expect_gt(max(abs(sj[["rx__sens_central_BY_eta_lag__"]])), 15)
  })

  test_that("raw event-table multiply with a modeled lag gets a dtau row (Phase B, B2)", {
    # Same gap as above, for evid=6 (multiply): checked on a compartment OTHER
    # than the multiplied one (depot's multiply shifts central's future
    # absorption input), since a single-compartment homogeneous-linear decay
    # makes the paper's dx1/dtau row (J_cc*x1 - f_c) identically (and
    # correctly) zero -- not a useful regression target for THIS row.
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- exp(tlag + eta_lag)
      depot(0) <- 100
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.3, tcl = 1, tv = 2, tlag = log(1.1), eta_lag = 0)
    e <- et(time = 5, amt = 0.6, cmt = "depot", evid = 6) |> et(seq(0, 15, 0.1))
    mj <- rxode2(.mod, calcSens = "eta_lag", eventSens = "jump")
    mfd <- rxode2(.mod, calcSens = "eta_lag", eventSens = "fd")
    sj <- rxSolve(mj, e, pars)
    h <- 1e-3
    pp <- pars; pp["eta_lag"] <- pars["eta_lag"] + h
    pm <- pars; pm["eta_lag"] <- pars["eta_lag"] - h
    sp <- rxSolve(mfd, e, pp, atol = 1e-11, rtol = 1e-11)
    sm <- rxSolve(mfd, e, pm, atol = 1e-11, rtol = 1e-11)
    fd <- (sp$central - sm$central) / (2 * h)
    .post <- sj$time >= 6.15 & sj$time < 15
    expect_equal(sj[["rx__sens_central_BY_eta_lag__"]][.post], fd[.post], tolerance = 1e-3)
    expect_gt(max(abs(sj[["rx__sens_central_BY_eta_lag__"]])), 0.01)
  })

  test_that("constant-rate infusion needs no jump (sens ODE alone is correct)", {
    # An infusion rate that does not depend on any parameter (and a fixed start
    # time) enters dydt as a parameter-independent forcing, so the symbolic
    # sensitivity ODE already captures it -- no jump contribution is required.
    # Both modes must agree with each other and with a finite difference; this
    # also guards that future infusion jumps stay zero in the constant-rate case.
    .mod <- "
      ka <- exp(tka + eta_ka)
      cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, eta_ka = 0)
    e <- et(amt = 100, rate = 20, cmt = "depot") |> et(seq(0.5, 12, 1))
    .central <- function(p, mode) {
      m <- rxode2(.mod, calcSens = "eta_ka", eventSens = mode)
      rxSolve(m, e, p)
    }
    sj <- .central(pars, "jump")[["rx__sens_central_BY_eta_ka__"]]
    sf <- .central(pars, "fd")[["rx__sens_central_BY_eta_ka__"]]
    h <- 1e-4
    pp <- pars; pp["eta_ka"] <- h
    pm <- pars; pm["eta_ka"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    expect_equal(sj, sf)                    # jump adds nothing for constant rate
    expect_equal(sj, fd, tolerance = 1e-3)  # and both match the finite difference
  })

  test_that("modeled-rate infusion with NO eventSens/calcSens requested does not crash", {
    # Regression for a real C stack overflow: `_esInfusionDydtPre()`
    # (inst/include/rxode2parseHandleEvid.h) used to call `dydtEs()` (an
    # alias for the model's OWN compiled dydt()) UNCONDITIONALLY at every
    # MODEL_RATE_ON/MODEL_DUR_ON event, before any `_rxEsActive` check.
    # Since `handle_evid()` is itself invoked FROM WITHIN `dydt()` for
    # modeled-rate/duration events, this created unbounded re-entrant
    # recursion into `dydt()` on EVERY solve of a model with modeled
    # rate()/dur() -- entirely independent of whether eventSens/calcSens
    # was ever requested (found via `git bisect`, culprit commit
    # ce61e9a99, then reproduced standalone on a completely unrelated,
    # pre-existing test file with no calcSens= anywhere). This model
    # deliberately requests NEITHER calcSens NOR eventSens, so
    # `_rxEsActive` is 0 for the whole solve -- the crash fired
    # regardless.
    .mod <- "
      ka <- 0.5; cl <- 0.2; v <- 10
      rate(central) <- 20
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    m <- rxode2(.mod)
    e <- et(amt = 100, cmt = "central", rate = -1) |> et(seq(0.5, 12, 1))
    expect_no_error(rxSolve(m, e))
    expect_no_error(rxSolve(m, e, method = "lsode"))
  })

  test_that("in-model evid_() dosing plugins match finite differences (jump)", {
    # bolus(), replace(), multiply(), reset() and constant-rate infuse()/
    # infuseDur() pushed from inside the model must all yield jump-mode
    # sensitivities that match a central difference of the solution.  (Modeled
    # rate/dur -- where the infusion magnitude depends on a parameter -- is the
    # separate continuous-forcing case covered elsewhere.)
    .skel <- function(dose) paste0("
      ka <- exp(tka + eta_ka); cl <- exp(tcl); v <- exp(tv); rt <- 5
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
      if (t >= rt && t < rt + 0.4) { ", dose, " }")
    .doses <- c(
      bolus      = "bolus(50, central, 0, 0, 0)",
      infuse     = "infuse(100, 50, central, 0, 0, 0)",
      infuseDur  = "infuseDur(100, 2, central, 0, 0, 0)",
      replace    = "if (central > 0) replace(40, central)",
      multiply   = "if (central > 0) multiply(0.5, central)",
      reset      = "reset()"
    )
    pars <- c(tka = 0.2, tcl = 1, tv = 2, eta_ka = 0)
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.5, 14, 0.5))
    h <- 1e-4
    for (.nm in names(.doses)) {
      .mod <- .skel(.doses[[.nm]])
      .central <- function(p, mode) {
        rxSolve(rxode2(.mod, calcSens = "eta_ka", eventSens = mode), e, p)
      }
      sj <- .central(pars, "jump")
      pp <- pars; pp["eta_ka"] <- h
      pm <- pars; pm["eta_ka"] <- -h
      fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
      .post <- sj$time >= 5
      expect_equal(sj[["rx__sens_central_BY_eta_ka__"]][.post], fd[.post],
                   tolerance = 1e-3,
                   info = paste0("plugin: ", .nm))
    }
  })

  test_that("modeled-rate infusion continuous forcing matches finite differences", {
    # rate(central) = exp(tr + eta_r) with a rate=-1 dose: the infusion rate
    # depends on eta_r.  The rate is solver-applied forcing (not a symbolic term
    # in f), so the symbolic sens ODE misses d(rate)/dp -- the jump injects it as
    # a forcing on the sensitivity compartment over the infusion.  The analytic
    # jump must match a central difference; the fd-mode sens ODE must NOT.
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      rate(central) <- exp(tr + eta_r)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tr = 0.7, eta_r = 0)
    e <- et(amt = 100, cmt = "central", rate = -1) |> et(seq(0.25, 14, 0.25))
    .central <- function(p, mode) {
      rxSolve(rxode2(.mod, calcSens = "eta_r", eventSens = mode), e, p)
    }
    sj <- .central(pars, "jump")[["rx__sens_central_BY_eta_r__"]]
    sf <- .central(pars, "fd")[["rx__sens_central_BY_eta_r__"]]
    h <- 1e-4
    pp <- pars; pp["eta_r"] <- h
    pm <- pars; pm["eta_r"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    expect_equal(sj, fd, tolerance = 1e-3)
    expect_gt(max(abs(sf - fd)), 1)   # fd-mode misses the forcing sensitivity
  })

  test_that("modeled-duration infusion (forcing + moving boundary) matches FD", {
    # dur(central) = exp(tr + eta_r) with a rate=-2 dose: rate = amt/dur depends
    # on eta_r AND the infusion end tau2 = tau1 + dur(eta_r) moves with eta_r.
    # The jump needs BOTH the continuous d(rate)/dp forcing and the moving-
    # boundary state jump rate*d(dur)/dp at the infusion end.  Validated vs FD;
    # the fd-mode sens ODE (neither contribution) is badly wrong.
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      dur(central) <- exp(tr + eta_r)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tr = 0.7, eta_r = 0)
    e <- et(amt = 100, cmt = "central", rate = -2) |> et(seq(0.25, 14, 0.25))
    .central <- function(p, mode) {
      rxSolve(rxode2(.mod, calcSens = "eta_r", eventSens = mode), e, p)
    }
    sj <- .central(pars, "jump")[["rx__sens_central_BY_eta_r__"]]
    sf <- .central(pars, "fd")[["rx__sens_central_BY_eta_r__"]]
    h <- 1e-4
    pp <- pars; pp["eta_r"] <- h
    pm <- pars; pm["eta_r"] <- -h
    fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
    expect_equal(sj, fd, tolerance = 1e-2)
    expect_gt(max(abs(sf - fd)), 1)   # fd-mode misses forcing + boundary
  })

  test_that("modeled-dur infusion with estimated F and duration matches FD", {
    # dur(central)=exp(tr+eta_r) AND f(central)=expit(tf+eta_f): rate=F*amt/dur
    # depends on BOTH parameters.  The continuous forcing carries both
    # d(rate)/dp pieces ((amt/dur)*dF and -(rate/dur)*dDur) and the moving
    # boundary applies only to the duration (tau2 = tau1 + dur, independent of F).
    # Both eta sensitivities must match a central difference.
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      dur(central) <- exp(tr + eta_r)
      f(central)   <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tr = 0.7, tf = 0.3,
              eta_r = 0, eta_f = 0)
    e <- et(amt = 100, cmt = "central", rate = -2) |> et(seq(0.25, 14, 0.25))
    .central <- function(p, mode) {
      rxSolve(rxode2(.mod, calcSens = c("eta_r", "eta_f"), eventSens = mode), e, p)
    }
    sj <- .central(pars, "jump")
    h <- 1e-4
    for (.eta in c("eta_r", "eta_f")) {
      pp <- pars; pp[.eta] <- h
      pm <- pars; pm[.eta] <- -h
      fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
      expect_equal(sj[[paste0("rx__sens_central_BY_", .eta, "__")]], fd,
                   tolerance = 1e-2, info = paste0("param: ", .eta))
    }
  })

  test_that("infusion with a modeled lag matches FD (all infusion types)", {
    # alag(central) = exp(tl + eta_l): the whole infusion window [tau1, tau2]
    # shifts by d(alag)/dp (tau1 = t0 + alag, tau2 = tau1 + dur), giving moving
    # start/stop boundary jumps [S] = -/+ rate*d(alag)/dp.  Validated for fixed-
    # rate, modeled-rate and modeled-duration infusions.  An irregular grid keeps
    # observations off the (discontinuous) boundary times so the FD is well-posed.
    .obs <- seq(0.1, 14, 0.37)
    .lagLine <- "alag(central) <- exp(tl + eta_l)"
    .ode <- "d/dt(depot) <- -ka * depot\n      d/dt(central) <- ka * depot - cl / v * central"
    .cases <- list(
      fixed = list(extra = "", dose = et(amt = 100, rate = 30, cmt = "central"),
                   p = c(tka = 0.2, tcl = 1, tv = 2, tl = log(1.1), eta_l = 0)),
      mrate = list(extra = "rate(central) <- exp(tr)",
                   dose = et(amt = 100, rate = -1, cmt = "central"),
                   p = c(tka = 0.2, tcl = 1, tv = 2, tl = log(1.1), tr = log(30), eta_l = 0)),
      mdur  = list(extra = "dur(central) <- exp(td)",
                   dose = et(amt = 100, rate = -2, cmt = "central"),
                   p = c(tka = 0.2, tcl = 1, tv = 2, tl = log(1.1), td = log(3), eta_l = 0))
    )
    h <- 1e-4
    for (.nm in names(.cases)) {
      .cs <- .cases[[.nm]]
      .mod <- paste("ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)",
                    .lagLine, .cs$extra, .ode, sep = "\n      ")
      e <- .cs$dose |> et(.obs)
      .central <- function(p, mode) {
        rxSolve(rxode2(.mod, calcSens = "eta_l", eventSens = mode), e, p)
      }
      sj <- .central(.cs$p, "jump")[["rx__sens_central_BY_eta_l__"]]
      pp <- .cs$p; pp["eta_l"] <- .cs$p["eta_l"] + h
      pm <- .cs$p; pm["eta_l"] <- .cs$p["eta_l"] - h
      fd <- (.central(pp, "fd")$central - .central(pm, "fd")$central) / (2 * h)
      expect_equal(sj, fd, tolerance = 1e-2, info = paste0("infusion: ", .nm))
    }
  })

  test_that("jump sensitivities are correct per subject (population solve)", {
    # FOCEi solves many subjects at once (in parallel).  Each subject has its
    # own etas, so its own lagged dose time and bioavailability; the jump
    # injection must use per-subject state/dose.  Three subjects with distinct
    # etas must each match a per-subject central difference.
    .mod <- "
      ka <- exp(tka + eta_ka); cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- exp(tlag + eta_lag)
      f(depot)    <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    th <- c(tka = 0.2, tcl = 1, tv = 2, tf = 0.3, tlag = log(1.1))
    ev <- data.frame(eta_ka = c(0, 0.3, -0.2),
                     eta_lag = c(0, 0.1, -0.15),
                     eta_f = c(0, -0.2, 0.25))
    e <- et(amt = 100, cmt = "depot", id = 1:3) |> et(seq(0.13, 12, 0.41), id = 1:3)
    mj <- rxode2(.mod, calcSens = c("eta_lag", "eta_f"), eventSens = "jump")
    mfd <- rxode2(.mod, calcSens = c("eta_lag", "eta_f"), eventSens = "fd")
    sj <- rxSolve(mj, e, th, iCov = ev)
    h <- 1e-4
    for (.sp in c("eta_lag", "eta_f")) {
      evp <- ev; evp[[.sp]] <- evp[[.sp]] + h
      evm <- ev; evm[[.sp]] <- evm[[.sp]] - h
      fd <- (rxSolve(mfd, e, th, iCov = evp)$central -
             rxSolve(mfd, e, th, iCov = evm)$central) / (2 * h)
      expect_equal(sj[[paste0("rx__sens_central_BY_", .sp, "__")]], fd,
                   tolerance = 1e-2, info = paste0("param: ", .sp))
    }
  })

  test_that(".rxEventSensMap exposes a second-order (Hessian) compartment map", {
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      f(depot) <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    m <- rxode2(.mod, calcSens = "eta_f", calcSens2 = "eta_f")
    im <- .rxEventSensMap(m)
    expect_false(is.null(im$map2))
    # one (state, p, q) row per physical state for the single (eta_f, eta_f) pair
    expect_setequal(im$map2$state, c("depot", "central"))
    expect_true(all(im$map2$p == "eta_f" & im$map2$q == "eta_f"))
    # the depot 2nd-order compartment exists and comes after the 1st-order ones
    .row <- im$map2[im$map2$state == "depot", ]
    expect_true(.row$sensCmt > nrow(im$map))
  })

  test_that(".rxEventSensCLines indexes multi-parameter calcSens2 consistently with calcSens", {
    # Regression: .rxEventSensMap() used to sort `map2` alphabetically by
    # (p, q, stateCmt); .rxEventSensCLines()'s .qIdx (built from
    # unique(map2$q)) then inherited that alphabetical order instead of
    # calcSens2's *as-passed* order (matching the compiled compartment
    # layout rxExpandSens2_ actually used) -- for any calcSens2 with more
    # than one parameter where the as-passed order isn't alphabetical, the
    # 2nd-order derivative values were written into the WRONG compartment
    # (found while validating the 2nd-order infusion jump, 2026-06-30;
    # masked in every earlier test, which all used a single-parameter
    # calcSens2 where reordering is a no-op). trate/tlag deliberately spelled
    # so alphabetical order ("tlag","trate") differs from the as-passed order
    # ("trate","tlag").
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- exp(tlag)
      rate(depot) <- exp(trate)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    m <- rxode2(.mod, calcSens = c("trate", "tlag"), calcSens2 = c("trate", "tlag"),
                eventSens = "jump")
    info <- m$eventSensInfo
    expect_equal(names(info$map$sensParams), NULL) # sanity: sensParams unnamed char vec
    expect_equal(info$map$sensParams, c("trate", "tlag"))
    # .qIdx must be built in the SAME (as-passed) order as .pIdx/sensParams,
    # not alphabetically -- this is the actual bug: previously "tlag" (unique(map2$q)'s
    # alphabetically-first value) would come first here instead of "trate".
    expect_equal(unique(info$map$map2$q), c("trate", "tlag"))
    expect_equal(unique(info$map$map2$p), c("trate", "tlag"))
  })

  test_that(".rxEventSensD2Expr is the correct second total derivative", {
    # d2(F)/d(eta_f)^2 for F = expit(tf + eta_f) must equal expit''(tf+eta_f)
    # = p(1-p)(1-2p).  State-independent F -> only the direct second partial.
    m <- rxode2("ka<-exp(tka)\nf(depot)<-expit(tf+eta_f)\nd/dt(depot)<- -ka*depot\nd/dt(central)<-ka*depot-central",
                calcSens = "eta_f", calcSens2 = "eta_f")
    mdl <- .rxLoadPrune(m)
    e2 <- .rxEventSensD2Expr(get("rx_f_depot_", envir = mdl), "eta_f", "eta_f",
                             .rxEventSensMap(m)$states)
    expect_true(nzchar(e2) && e2 != "0")
    val <- eval(parse(text = e2),
                list(eta_f = 0, tf = 0.3, Rx_pow_di = function(a, b) a^b,
                     expit = function(x) 1 / (1 + exp(-x))))
    p <- 1 / (1 + exp(-0.3))
    expect_equal(val, p * (1 - p) * (1 - 2 * p), tolerance = 1e-8)
  })

  test_that("calcSens2 generates correct second-order sensitivity ODEs", {
    # rxode2(..., calcSens2=) emits the 2nd-order sensitivity compartments
    # rx__sens_<state>_BY_<p>_BY_<q>__ (the Hessian path).  For a parameter that
    # does not touch dosing there is no jump, so the 2nd-order sensitivity comes
    # purely from the continuous variational ODEs.  S^{pq} = d(S^p)/dq, so it must
    # match a (precise) first central difference of the first-order sensitivity.
    .mod <- "
      ka <- exp(tka + eta_ka); cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, eta_ka = 0)
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.5, 12, 0.5))
    m1 <- rxode2(.mod, calcSens = "eta_ka")
    m2 <- rxode2(.mod, calcSens = "eta_ka", calcSens2 = "eta_ka")
    expect_true("rx__sens_central_BY_eta_ka_BY_eta_ka__" %in% rxModelVars(m2)$state)
    s2 <- rxSolve(m2, e, pars)[["rx__sens_central_BY_eta_ka_BY_eta_ka__"]]
    h <- 1e-5
    .s1 <- function(p) rxSolve(m1, e, p)[["rx__sens_central_BY_eta_ka__"]]
    pp <- pars; pp["eta_ka"] <- h
    pm <- pars; pm["eta_ka"] <- -h
    fd <- (.s1(pp) - .s1(pm)) / (2 * h)
    expect_equal(s2, fd, tolerance = 1e-2)
  })

  test_that("calcSens3 generates correct third-order sensitivity ODEs (Phase H0)", {
    # rxode2(..., calcSens2=, calcSens3=) emits the 3rd-order sensitivity
    # compartments rx__sens_<state>_BY_<p>_BY_<q>_BY_<r>__ via
    # rxExpandSens3_() (PR #1092), exposed at the rxode2() build level for
    # the first time here (previously only reachable through the DDE-
    # specific delay() sensitivity path). No jump-sensitivity content is
    # involved -- this is pure continuous-ODE 3rd-order sensitivity, so it
    # is validated standalone against a nested finite difference of the
    # analytic 2nd-order sensitivity (the plan's H0 acceptance criterion),
    # with calcSens2 == calcSens3 == calcSens (the common "full Hessian/
    # third-order" case, matching how calcSens2 is used everywhere else).
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2)
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.5, 12, 1))
    m2 <- rxode2(.mod, calcSens = "tka", calcSens2 = "tka")
    m3 <- rxode2(.mod, calcSens = "tka", calcSens2 = "tka", calcSens3 = "tka")
    expect_true("rx__sens_central_BY_tka_BY_tka_BY_tka__" %in% rxModelVars(m3)$state)
    s3 <- rxSolve(m3, e, pars)[["rx__sens_central_BY_tka_BY_tka_BY_tka__"]]
    h <- 1e-4
    .s2 <- function(p) rxSolve(m2, e, p)[["rx__sens_central_BY_tka_BY_tka__"]]
    pp <- pars; pp["tka"] <- pars["tka"] + h
    pm <- pars; pm["tka"] <- pars["tka"] - h
    fd <- (.s2(pp) - .s2(pm)) / (2 * h)
    expect_equal(s3, fd, tolerance = 1e-3)

    # calcSens3 requires calcSens2
    expect_error(rxode2(.mod, calcSens = "tka", calcSens3 = "tka"))
  })

  test_that("second-order additive-bolus F jump matches finite differences", {
    # F = expit(tf + eta_f) feeds the dose only (ka/cl/v independent of eta_f),
    # so the second-order sensitivity S^{ff} of the dosed compartment is *entirely*
    # the explicit jump amt*d2F/d eta_f^2 -- without it S^{ff} is exactly 0.  The
    # jump-mode 2nd-order sensitivity must match d(S^f)/d eta_f (a precise first
    # central difference of the jump-correct first-order sensitivity).
    .mod <- "
      ka <- exp(tka); cl <- exp(tcl); v <- exp(tv)
      alag(depot) <- 1
      f(depot)    <- expit(tf + eta_f)
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(tka = 0.2, tcl = 1, tv = 2, tf = 0.3, eta_f = 0)
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.5, 12, 0.5))
    m1 <- rxode2(.mod, calcSens = "eta_f", eventSens = "jump")
    m2 <- rxode2(.mod, calcSens = "eta_f", calcSens2 = "eta_f", eventSens = "jump")
    s2 <- rxSolve(m2, e, pars)[["rx__sens_central_BY_eta_f_BY_eta_f__"]]
    h <- 1e-5
    .s1 <- function(p) rxSolve(m1, e, p)[["rx__sens_central_BY_eta_f__"]]
    pp <- pars; pp["eta_f"] <- h
    pm <- pars; pm["eta_f"] <- -h
    fd <- (.s1(pp) - .s1(pm)) / (2 * h)
    expect_equal(s2, fd, tolerance = 1e-2)
    # without the 2nd-order jump the contribution would be entirely missing
    expect_gt(max(abs(fd)), 1)
  })

  test_that("third-order additive-bolus F jump matches nested finite differences (Phase H1)", {
    # F = expit(tf) feeds the dose only, so the third-order sensitivity
    # S^{fff} of the dosed compartment is entirely the explicit jump
    # amt*d3F/dtf^3 -- validate against a central difference of the
    # jump-correct *second*-order sensitivity (mirrors the 2nd-order F-jump
    # test's own validation strategy one level deeper).
    ode_code_f <- "
      f(depot)    <- expit(tf)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    pars_f <- c(ka = 0.5, cl = 0.2, v = 10, tf = qlogis(0.7))
    e <- et(amt = 100, cmt = "depot") |> et(seq(0, 10, by = 0.5))
    m2 <- rxode2(ode_code_f, calcSens = c("tf", "ka"), calcSens2 = c("tf", "ka"),
                 eventSens = "jump")
    m3 <- rxode2(ode_code_f, calcSens = c("tf", "ka"), calcSens2 = c("tf", "ka"),
                 calcSens3 = "tf", eventSens = "jump")
    expect_true("rx__sens_central_BY_tf_BY_tf_BY_tf__" %in% rxModelVars(m3)$state)
    s3 <- rxSolve(m3, e, pars_f)[["rx__sens_central_BY_tf_BY_tf_BY_tf__"]]
    eps <- 1e-4
    p1 <- pars_f; p1["tf"] <- pars_f["tf"] + eps
    p2 <- pars_f; p2["tf"] <- pars_f["tf"] - eps
    .s2 <- function(p) rxSolve(m2, e, p)[["rx__sens_central_BY_tf_BY_tf__"]]
    fd3 <- (.s2(p1) - .s2(p2)) / (2 * eps)
    expect_equal(s3, fd3, tolerance = 1e-4)
    # without the 3rd-order jump the contribution would be entirely missing
    expect_gt(max(abs(fd3)), 1)
  })

  test_that("third-order additive-bolus F jump: replace/multiply scale/zero the 3rd-order compartment", {
    # Regression mirroring the 2nd-order replace/multiply fix: a raw
    # event-table replace() zeros every sensitivity order of the replaced
    # state (constant replacement value); a multiply() scales every order by
    # the same alpha.  Both validated against nested FD of the jump-correct
    # 2nd-order sensitivity.
    ode_code_f <- "
      f(depot)    <- expit(tf)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    pars_f <- c(ka = 0.5, cl = 0.2, v = 10, tf = qlogis(0.7))
    m2 <- rxode2(ode_code_f, calcSens = c("tf", "ka"), calcSens2 = c("tf", "ka"),
                 eventSens = "jump")
    m3 <- rxode2(ode_code_f, calcSens = c("tf", "ka"), calcSens2 = c("tf", "ka"),
                 calcSens3 = "tf", eventSens = "jump")
    for (.evid in c(5, 6)) {
      .amt <- if (.evid == 5) 50 else 0.5
      e <- et(amt = 100, cmt = "depot") |>
        et(time = 5, amt = .amt, cmt = "central", evid = .evid) |>
        et(seq(0.5, 12, 1))
      s3 <- rxSolve(m3, e, pars_f)[["rx__sens_central_BY_tf_BY_tf_BY_tf__"]]
      eps <- 1e-4
      p1 <- pars_f; p1["tf"] <- pars_f["tf"] + eps
      p2 <- pars_f; p2["tf"] <- pars_f["tf"] - eps
      .s2 <- function(p) rxSolve(m2, e, p)[["rx__sens_central_BY_tf_BY_tf__"]]
      fd3 <- (.s2(p1) - .s2(p2)) / (2 * eps)
      expect_equal(s3, fd3, tolerance = 1e-4, info = paste0("evid: ", .evid))
    }
  })

  test_that("second-order dtau/lag row: q unrelated to the modeled alag (Phase H1-dtau)", {
    # tlag drives alag; tf drives F only (independent of tlag).  Differentiating
    # the 1st-order dtau jump (-J[k][c]*delta*dLag_p[c]) by the product rule
    # wrt q=tf needs three terms: d(J)/dq (here 0, linear system), d(delta)/dq
    # = amt*dFQ[c][q] (nonzero, F depends on tf), and d2Lag[p][q] (here 0,
    # alag doesn't depend on tf).  Only the middle term survives, and it is
    # exactly what makes S^{tlag,tf} nonzero (without it, entirely missing).
    ode_code <- "
      alag(depot) <- exp(tlag)
      f(depot)    <- expit(tf)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl / v * central
    "
    pars <- c(ka = 0.5, cl = 0.2, v = 10, tlag = log(0.5), tf = qlogis(0.7))
    # avoid the observation landing exactly on the (parameter-dependent) lag
    # time (t=0.5): a known artifact source, see the plan's validation notes.
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.55, 10, by = 0.5))
    m1 <- rxode2(ode_code, calcSens = c("tlag", "tf"), eventSens = "jump")
    m2 <- rxode2(ode_code, calcSens = c("tlag", "tf"), calcSens2 = c("tlag", "tf"),
                 eventSens = "jump")
    s2 <- rxSolve(m2, e, pars, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag_BY_tf__"]]
    eps <- 1e-5
    p1 <- pars; p1["tf"] <- pars["tf"] + eps
    p2 <- pars; p2["tf"] <- pars["tf"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    expect_gt(max(abs(fd)), 1)
  })

  test_that("second-order dtau/lag row: d(J[k][c])/dq Jacobian-coupling term (nonlinear model)", {
    # Michaelis-Menten elimination makes the physical Jacobian state- and
    # parameter-dependent, exercising the d(J)/dq term (zero for the linear
    # models above) in isolation: J[central][central] = -Vm*Km/(Km+central)^2
    # depends explicitly on Vm, so d(J)/dtvm != 0.
    ode_code <- "
      alag(central) <- exp(tlag)
      vm <- exp(tvm)
      km <- exp(tkm)
      d/dt(central) = -vm*central/(km+central)
    "
    pars <- c(tvm = log(2), tkm = log(5), tlag = log(0.5))
    e <- et(amt = 10, cmt = "central") |> et(seq(0.55, 8, by = 0.5))
    m1 <- rxode2(ode_code, calcSens = c("tlag", "tvm"), eventSens = "jump")
    m2 <- rxode2(ode_code, calcSens = c("tlag", "tvm"), calcSens2 = c("tlag", "tvm"),
                 eventSens = "jump")
    info <- m2$eventSensInfo
    expect_true(nrow(info$derivs$lagJacQ) > 0L)
    s2 <- rxSolve(m2, e, pars, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag_BY_tvm__"]]
    eps <- 1e-5
    p1 <- pars; p1["tvm"] <- pars["tvm"] + eps
    p2 <- pars; p2["tvm"] <- pars["tvm"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    expect_gt(max(abs(fd)), 0.1)
  })

  test_that("second-order dtau/lag row: Leibniz/moving-boundary term when q ALSO drives the same alag (Phase H1-dtau)", {
    # Differentiating the 1st-order dtau row by the product rule alone is
    # PROVEN INCOMPLETE (by FD) whenever q ALSO shifts the same event's time
    # (e.g. the (tlag,tlag) diagonal, or two lag parameters combined
    # additively in one alag() expression): it misses a Leibniz/moving-
    # boundary term `-[dS^p_k/dt|post - dS^p_k/dt|pre] * dLag_q[c]`, the same
    # class of problem that originally blocked the infusion 2nd-order
    # attempt. `dS^p_k/dt` is read directly from the compiled dydt() output
    # (pre-jump state vs. a scratch post-jump copy), needing no separate
    # Jacobian-matrix construction. Validated against FD for both the
    # diagonal and a two-lag-parameter model.
    ode_code_diag <- "
      alag(depot) <- exp(tlag)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl / v * central
    "
    pars_diag <- c(ka = 0.5, cl = 0.2, v = 10, tlag = log(0.5))
    e_diag <- et(amt = 100, cmt = "depot") |> et(seq(0.55, 10, by = 0.5))
    m1d <- rxode2(ode_code_diag, calcSens = "tlag", eventSens = "jump")
    m2d <- rxode2(ode_code_diag, calcSens = "tlag", calcSens2 = "tlag", eventSens = "jump")
    s2d <- rxSolve(m2d, e_diag, pars_diag, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag_BY_tlag__"]]
    eps <- 1e-5
    p1 <- pars_diag; p1["tlag"] <- pars_diag["tlag"] + eps
    p2 <- pars_diag; p2["tlag"] <- pars_diag["tlag"] - eps
    .s1d <- function(p) rxSolve(m1d, e_diag, p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag__"]]
    fdd <- (.s1d(p1) - .s1d(p2)) / (2 * eps)
    expect_equal(s2d, fdd, tolerance = 1e-6)
    expect_gt(max(abs(fdd)), 1)

    ode_code_2p <- "
      alag(depot) <- exp(tlag1 + tlag2)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl / v * central
    "
    pars_2p <- c(ka = 0.5, cl = 0.2, v = 10, tlag1 = log(0.5) / 2, tlag2 = log(0.5) / 2)
    e_2p <- et(amt = 100, cmt = "depot") |> et(seq(0.55, 10, by = 0.25))
    m1p <- rxode2(ode_code_2p, calcSens = c("tlag1", "tlag2"), eventSens = "jump")
    m2p <- rxode2(ode_code_2p, calcSens = c("tlag1", "tlag2"), calcSens2 = c("tlag1", "tlag2"),
                  eventSens = "jump")
    s2p <- rxSolve(m2p, e_2p, pars_2p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag1_BY_tlag2__"]]
    p1p <- pars_2p; p1p["tlag2"] <- pars_2p["tlag2"] + eps
    p2p <- pars_2p; p2p["tlag2"] <- pars_2p["tlag2"] - eps
    .s1p <- function(p) rxSolve(m1p, e_2p, p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag1__"]]
    fdp <- (.s1p(p1p) - .s1p(p2p)) / (2 * eps)
    expect_equal(s2p, fdp, tolerance = 1e-6)
    expect_gt(max(abs(fdp)), 1)
  })

  test_that("second-order dtau/lag row: mirror compartment symmetry (Phase H1-dtau)", {
    # S^{p,q} and S^{q,p} must be identical (Schwarz).  Before the Leibniz
    # term was added, S^{tlag,tf} validated correctly (tlag drives alag, so
    # the outer p-loop ran) while S^{tf,tlag} stayed silently at 0 (tf does
    # not drive alag, so the outer p-loop skipped it entirely and never had
    # a chance to apply the Leibniz correction to S^tf's own post-ddelta-jump
    # value). The p-loop now also runs whenever ANY calcSens2 parameter
    # drives this compartment's alag, regardless of the current p's own
    # dLag_p, fixing this.
    ode_code <- "
      alag(depot) <- exp(tlag)
      f(depot)    <- expit(tf)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl / v * central
    "
    pars <- c(ka = 0.5, cl = 0.2, v = 10, tlag = log(0.5), tf = qlogis(0.7))
    e <- et(amt = 100, cmt = "depot") |> et(seq(0.55, 10, by = 0.5))
    m2 <- rxode2(ode_code, calcSens = c("tlag", "tf"), calcSens2 = c("tlag", "tf"),
                 eventSens = "jump")
    s2 <- rxSolve(m2, e, pars, atol = 1e-12, rtol = 1e-12)
    expect_equal(s2$rx__sens_central_BY_tlag_BY_tf__, s2$rx__sens_central_BY_tf_BY_tlag__,
                 tolerance = 1e-8)
    expect_gt(max(abs(s2$rx__sens_central_BY_tlag_BY_tf__)), 1)
  })

  test_that("second-order infusion-boundary jump (fixed-rate/duration infusion + modeled alag)", {
    # Re-derives the same "jump condition one level up" trick used for the
    # additive-bolus dtau row, applied to a forcing DISCONTINUITY (not a
    # state jump): a fixed-rate/duration infusion whose start/stop time is
    # shifted by a modeled alag jumps InfusionRate[cmt] by a fixed amount at
    # the lag-shifted boundary. The naive product-rule-only term
    # (tmp*d2Lag[p][q]) matched FD DURING the infusion window and diverged
    # right after it ended (documented in the plan as "found wrong, reverted"
    # in an earlier session) -- fixed here by adding the missing Leibniz
    # term, validated through BOTH the start and stop boundaries.
    ode_code <- "
      alag(central) <- exp(tlag)
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    pars <- c(ka = 0.5, cl = 0.2, v = 10, tlag = log(0.5))
    e <- et(amt = 100, cmt = "central", rate = 20) |> et(seq(0.55, 15, by = 0.5))
    m1 <- rxode2(ode_code, calcSens = "tlag", eventSens = "jump")
    m2 <- rxode2(ode_code, calcSens = "tlag", calcSens2 = "tlag", eventSens = "jump")
    s2 <- rxSolve(m2, e, pars, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag_BY_tlag__"]]
    eps <- 1e-5
    p1 <- pars; p1["tlag"] <- pars["tlag"] + eps
    p2 <- pars; p2["tlag"] <- pars["tlag"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-12, rtol = 1e-12)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    # crosses the stop boundary at tlag+amt/rate=5.5; matched only DURING the
    # window before the Leibniz term was added, so require post-window rows
    # to be present and nonzero (not just the pre-window ones).
    expect_gt(sum(e$time[e$evid == 0] > 6), 5)
    expect_gt(max(abs(fd)), 1)
  })

  test_that("second-order infusion-boundary jump propagates to a COUPLED (non-infused) compartment", {
    # The dosed/lag-shifted compartment's OWN sensitivity jump [S^p]_cmt is
    # the only *value* discontinuity at 1st order (other states only pick up
    # a kink -- their derivative, not their value, jumps, since it's a
    # forcing discontinuity, not a state jump like the additive-bolus case).
    # But at 2nd order, differentiating that kink's own (parameter-dependent)
    # strength wrt a SECOND parameter that ALSO shifts the boundary time
    # reintroduces a genuine Leibniz value-jump for the coupled compartment
    # too -- missed in an earlier version of this fix that only updated the
    # k=cmt (dosed) compartment, because every prior validation happened to
    # observe the SAME compartment being infused (an accidental blind spot).
    # This model infuses "depot" (coupled to "central" via ka) and checks
    # central's OWN 2nd-order sensitivity, exercising exactly that path.
    ode_code <- "
      alag(depot) <- tlag
      f(depot)    <- doseAmt
      dur(depot)  <- tinf
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    pars <- c(ka = 1, cl = 6, v = 60, tlag = 10, doseAmt = 200, tinf = 10)
    e <- et(amt = 1, cmt = "depot", rate = -2) |> et(seq(0.53, 60, by = 0.5))
    m1 <- rxode2(ode_code, calcSens = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    m2 <- rxode2(ode_code, calcSens = c("tlag", "doseAmt", "tinf"),
                 calcSens2 = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    s2 <- rxSolve(m2, e, pars, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag_BY_tlag__"]]
    eps <- 1e-4
    p1 <- pars; p1["tlag"] <- pars["tlag"] + eps
    p2 <- pars; p2["tlag"] <- pars["tlag"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    expect_gt(max(abs(fd)), 1)
  })

  test_that("additive-bolus F jump fires for indexed THETA[n]/ETA[n] params", {
    # The nlmixr2 FOCEi inner model writes sensitivities wrt the indexed
    # parameters ETA[n]/THETA[n] (compartments rx__sens_<state>_BY_ETA_n___).
    # rxFromSE renders these as `ETA[3]`, whose all.vars() collapses to "ETA" and
    # never matched the SE-mangled sensParam "ETA_3_" -- the dosing-parameter
    # derivatives were silently dropped and no jump was injected.  Guard both the
    # symbolic detection (free_symbols) and the codegen rewrite (ETA[n]->_ETA_n_).
    code <- paste(
      "param(THETA[1],THETA[2],THETA[3],THETA[4],ETA[1],ETA[3]);",
      "cmt(depot);",
      "cmt(central);",
      "f(depot)=1/(1+exp(-(ETA[3]+THETA[4])));",
      "d/dt(depot)=-exp(ETA[1]+THETA[1])*depot;",
      "d/dt(central)=exp(ETA[1]+THETA[1])*depot-exp(THETA[2]-THETA[3])*central;",
      "d/dt(rx__sens_depot_BY_ETA_3___)=-exp(ETA[1]+THETA[1])*rx__sens_depot_BY_ETA_3___;",
      "d/dt(rx__sens_central_BY_ETA_3___)=exp(ETA[1]+THETA[1])*rx__sens_depot_BY_ETA_3___-exp(THETA[2]-THETA[3])*rx__sens_central_BY_ETA_3___;",
      sep = "\n")
    mj <- rxode2(code, eventSens = "jump")
    # the F derivative wrt ETA[3] must be detected (non-empty derivs table)
    .df <- mj$eventSensInfo$derivs$f
    expect_true(nrow(.df) >= 1L)
    expect_true("ETA_3_" %in% .df$param)

    p <- c("THETA[1]" = 0.45, "THETA[2]" = 1, "THETA[3]" = 3.45,
           "THETA[4]" = 0.9, "ETA[1]" = 0, "ETA[3]" = 0.2)
    # sample immediately after the dose to capture the jump peak before decay
    e <- et(amt = 100, cmt = "depot") |> et(c(1e-4, seq(0.13, 24, length.out = 30)))
    s <- rxSolve(mj, p, e, atol = 1e-10, rtol = 1e-10)
    # the additive-bolus F jump makes the depot sensitivity wrt ETA[3] non-zero;
    # analytic value at t=0+ is amt*F*(1-F) with F = expit(0.2 + 0.9).
    .F <- 1 / (1 + exp(-(0.2 + 0.9)))
    expect_equal(max(s[["rx__sens_depot_BY_ETA_3___"]]), 100 * .F * (1 - .F),
                 tolerance = 1e-2)
  })

  test_that("eventSens mode is folded into the model cache key", {
    mfd <- rxode2(.modConstF, calcSens = c("eta_ka", "eta_lag"),
                  eventSens = "fd")
    mj <- rxode2(.modConstF, calcSens = c("eta_ka", "eta_lag"),
                 eventSens = "jump")
    # same model text, different mode -> distinct parsed md5 (distinct DLLs)
    expect_false(identical(
      unname(rxModelVars(mfd)$md5["parsed_md5"]),
      unname(rxModelVars(mj)$md5["parsed_md5"])
    ))
    # the cache key is reset after each build (no leak to later non-jump builds)
    expect_identical(.rxEventSensCacheKey, "")
  })

  test_that("reproduces the paper's PK/PD IMAX infusion example (dEffect/dtlag, ddose, dtinf)", {
    # Kaschek & Fidler, "Forward Sensitivity Equations in the Presence of
    # Events" (~/src/EventSensitivities/paper.tex, Sec. Example): a depot-
    # absorption model with a zero-order infusion into the depot compartment,
    # parameterized directly by lag time (tlag), infusion duration (tinf), and
    # dose (not etas) -- alag(depot)=tlag, dur(depot)=tinf, f(depot)=doseAmt
    # (bioavailability standing in for the dose, routed through the
    # duration-fixed F/amt convention: rate = F*amt/dur = doseAmt*1/tinf,
    # matching the paper's r1 = dose/tinf exactly). Effect is a downstream LHS
    # of central, so dEffect/dp = dEffect/dcentral * dcentral/dp (closed-form
    # chain rule below) using the jump-computed rx__sens_central_BY_p__.
    .mod <- "
      ka <- exp(lka)
      cl <- exp(lcl)
      v  <- exp(lv)
      imax <- exp(limax)
      ic50 <- exp(lic50)
      e0 <- exp(le0)
      alag(depot) <- tlag
      dur(depot)  <- tinf
      f(depot)    <- doseAmt
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
      effect <- e0 * (1 - (central / v) * imax / (ic50 + (central / v)))
    "
    pars <- c(
      lka = log(1), lcl = log(6), lv = log(60), limax = log(1), lic50 = log(1),
      le0 = log(15), tlag = 10, doseAmt = 200, tinf = 10
    )
    mj <- rxode2(.mod, calcSens = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    e <- et(amt = 1, cmt = "depot", rate = -2) |> et(seq(0, 60, by = 0.5))
    sj <- rxSolve(mj, e, pars)

    dEffect_dCentral <- with(sj, -e0 * imax * ic50 / (v * (ic50 + central / v)^2))
    d_tlag <- dEffect_dCentral * sj$rx__sens_central_BY_tlag__
    d_tinf <- dEffect_dCentral * sj$rx__sens_central_BY_tinf__
    d_dose <- dEffect_dCentral * sj$rx__sens_central_BY_doseAmt__

    # (1) matches a true finite difference of the solved "effect" output
    # (the plan's primary correctness gate, Section 5) for all three params.
    .h <- 1e-4
    .fd <- function(pname) {
      .pp <- pars; .pp[pname] <- pars[pname] + .h
      .pm <- pars; .pm[pname] <- pars[pname] - .h
      (rxSolve(mj, e, .pp)$effect - rxSolve(mj, e, .pm)$effect) / (2 * .h)
    }
    expect_equal(d_tlag, .fd("tlag"), tolerance = 1e-3)
    expect_equal(d_tinf, .fd("tinf"), tolerance = 1e-3)
    expect_equal(d_dose, .fd("doseAmt"), tolerance = 1e-3)

    # (2) qualitative claims from the paper (Sec. Example):
    # "higher doses lead to stronger inhibition" -- Effect decreases (relative
    # to E0) as dose increases, i.e. dEffect/ddose <= 0 everywhere.
    expect_true(all(d_dose <= 1e-8))
    # "a longer infusion time will decrease the inhibition during the
    # infusion and increase it afterwards. The same holds for [tlag]" --
    # Effect increases (dEffect/dp > 0) during [tlag, tlag+tinf] and
    # decreases (dEffect/dp < 0) after tlag+tinf, for both tlag and tinf.
    .during <- sj$time > pars["tlag"] + 0.5 & sj$time < pars["tlag"] + pars["tinf"] - 0.5
    .after <- sj$time > pars["tlag"] + pars["tinf"] + 0.5 & sj$time < 40
    expect_true(all(d_tinf[.during] > 0))
    expect_true(all(d_tinf[.after] < 0))
    expect_true(all(d_tlag[.during] > 0))
    expect_true(all(d_tlag[.after] < 0))
    # "the impact of the lag time is five times higher than the impact of the
    # infusion time" -- loose bound around the paper's approximate reading.
    .ratio <- max(abs(d_tlag)) / max(abs(d_tinf))
    expect_true(.ratio > 3 && .ratio < 10)
  })

  test_that("reproduces the paper's PK/PD IMAX infusion example at 2nd order (S^{tlag,tlag})", {
    # Same model as the 1st-order reproduction above, extended to calcSens2:
    # exercises the infusion-boundary Leibniz term together with the
    # additive-bolus dtau row's own Leibniz term in the paper's actual
    # headline scenario (depot infused, central observed, alag+F+dur all
    # modeled and all parameter-dependent). Offset time grid avoids landing
    # an observation exactly on the tlag/tlag+tinf boundaries (a known FD
    # artifact source, not a model error).
    .mod <- "
      ka <- exp(lka)
      cl <- exp(lcl)
      v  <- exp(lv)
      alag(depot) <- tlag
      dur(depot)  <- tinf
      f(depot)    <- doseAmt
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(lka = log(1), lcl = log(6), lv = log(60), tlag = 10, doseAmt = 200, tinf = 10)
    m1 <- rxode2(.mod, calcSens = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    m2 <- rxode2(.mod, calcSens = c("tlag", "doseAmt", "tinf"),
                 calcSens2 = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    e <- et(amt = 1, cmt = "depot", rate = -2) |> et(seq(0.53, 60, by = 0.5))
    s2 <- rxSolve(m2, e, pars, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag_BY_tlag__"]]
    eps <- 1e-4
    p1 <- pars; p1["tlag"] <- pars["tlag"] + eps
    p2 <- pars; p2["tlag"] <- pars["tlag"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    expect_gt(max(abs(fd)), 1)
  })

  test_that("second-order modeled-DUR boundary: combined alag+dur cross term S^{tlag,tinf}", {
    # alag(depot)+dur(depot), NO F: isolates the MODEL_DUR_ON/OFF boundary's
    # combined derivation (tau2 = tau1(alag) + dur(tinf) shifts with BOTH
    # mechanisms at once) from the additive-bolus dtau row's own Leibniz
    # term. An earlier implementation handled alag and dur as two SEPARATE
    # per-mechanism terms and got this cross term wrong (not just missing) --
    # confirmed by FD, ~8% relative error -- because dtau2/dp and dtau2/dq
    # must be SUMMED (dLag_p+dDur_p, dLagQ_q+dDurQ_q) before use, not treated
    # independently. Checks BOTH the infused (depot) and downstream
    # (central) compartments through and after the boundary.
    .mod <- "
      alag(depot) <- tlag
      dur(depot)  <- tinf
      d/dt(depot)   = -ka * depot
      d/dt(central) =  ka * depot - cl/v * central
    "
    pars <- c(ka = 1, cl = 6, v = 60, tlag = 10, tinf = 10)
    m1 <- rxode2(.mod, calcSens = c("tlag", "tinf"), eventSens = "jump")
    m2 <- rxode2(.mod, calcSens = c("tlag", "tinf"), calcSens2 = c("tlag", "tinf"),
                 eventSens = "jump")
    e <- et(amt = 1, cmt = "depot", rate = -2) |> et(seq(0.53, 60, by = 0.5))
    r2 <- rxSolve(m2, e, pars, atol = 1e-11, rtol = 1e-11)
    eps <- 1e-4
    p1 <- pars; p1["tinf"] <- pars["tinf"] + eps
    p2 <- pars; p2["tinf"] <- pars["tinf"] - eps
    r1a <- rxSolve(m1, e, p1, atol = 1e-11, rtol = 1e-11)
    r1b <- rxSolve(m1, e, p2, atol = 1e-11, rtol = 1e-11)
    fd_central <- (r1a$rx__sens_central_BY_tlag__ - r1b$rx__sens_central_BY_tlag__) / (2 * eps)
    fd_depot <- (r1a$rx__sens_depot_BY_tlag__ - r1b$rx__sens_depot_BY_tlag__) / (2 * eps)
    expect_equal(r2$rx__sens_central_BY_tlag_BY_tinf__, fd_central, tolerance = 1e-6)
    expect_equal(r2$rx__sens_depot_BY_tlag_BY_tinf__, fd_depot, tolerance = 1e-6)
    expect_gt(max(abs(fd_central)), 0.05)
  })

  test_that("second-order modeled-DUR boundary: alag x F cross term S^{tlag,doseAmt}", {
    # Full alag+dur+F IMAX model: exercises the SAME combined boundary
    # derivation as the alag+dur-only test above, but with q (doseAmt)
    # driving F instead of dur. A first fix (validated against the
    # alag+dur-only isolate) still had this pair wrong by a full sign flip:
    # the dRateQq helper (d(rate)/dq used in both the MODEL_DUR_ON tau1-
    # boundary term and the MODEL_DUR_OFF tau2-boundary term) had its dur
    # component's sign backwards relative to the already-validated 1st-order
    # d(rate)/dp = (amt*dF + tmp*dDur)/dur formula -- it happened to cancel
    # out for a pure-dur q (like tinf above) but flipped the answer for a
    # pure-F q (doseAmt), so both cross terms are needed to catch this class
    # of bug.
    .mod <- "
      ka <- exp(lka)
      cl <- exp(lcl)
      v  <- exp(lv)
      alag(depot) <- tlag
      dur(depot)  <- tinf
      f(depot)    <- doseAmt
      d/dt(depot)   <- -ka * depot
      d/dt(central) <-  ka * depot - cl / v * central
    "
    pars <- c(lka = log(1), lcl = log(6), lv = log(60), tlag = 10, doseAmt = 200, tinf = 10)
    m1 <- rxode2(.mod, calcSens = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    m2 <- rxode2(.mod, calcSens = c("tlag", "doseAmt", "tinf"),
                 calcSens2 = c("tlag", "doseAmt", "tinf"), eventSens = "jump")
    e <- et(amt = 1, cmt = "depot", rate = -2) |> et(seq(0.53, 60, by = 0.5))
    s2 <- rxSolve(m2, e, pars, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag_BY_doseAmt__"]]
    eps <- 1e-4
    p1 <- pars; p1["doseAmt"] <- pars["doseAmt"] + eps
    p2 <- pars; p2["doseAmt"] <- pars["doseAmt"] - eps
    .s1 <- function(p) rxSolve(m1, e, p, atol = 1e-11, rtol = 1e-11)[["rx__sens_central_BY_tlag__"]]
    fd <- (.s1(p1) - .s1(p2)) / (2 * eps)
    expect_equal(s2, fd, tolerance = 1e-6)
    expect_gt(max(abs(fd)), 0.05)
  })

  test_that("modeled dur() STEADY STATE forward jump-sensitivities match FD (solveSSinf moving boundary)", {
    # solveSSinf re-expresses a modeled dur() ss infusion as a fixed-rate window;
    # its classic OFF used to skip the MODEL_DUR_OFF sens forcing-removal + moving-
    # boundary jump, so post-boundary sensitivities were ~30% wrong.  handleSS now
    # arms a marker so handle_evid runs that logic at the re-expressed OFF.
    mt <- "d/dt(depot)=-ka*depot\nd/dt(central)=ka*depot-(cl/v)*central\ndur(central)=9*cl/3.5\ncp=central/(v/1000)"
    cs <- c("ka", "cl"); p <- c(ka = 1.2, cl = 3.5, v = 25)
    mb <- rxode2(mt); mj <- rxode2(mt, calcSens = cs, eventSens = "jump")
    ev <- et(amt = 100, rate = -2, cmt = "central", ss = 1, ii = 24) |> et(c(4, 8, 9.5, 12, 16, 24, 30))
    f <- as.data.frame(suppressWarnings(rxSolve(mj, ev, params = p, atol = 1e-11, rtol = 1e-11)))
    mx <- 0
    for (pn in cs) {
      hh <- abs(p[[pn]]) * 1e-6; pp <- p; pm <- p; pp[pn] <- pp[pn] + hh; pm[pn] <- pm[pn] - hh
      sp <- as.data.frame(suppressWarnings(rxSolve(mb, ev, params = pp, atol = 1e-11, rtol = 1e-11)))
      sm <- as.data.frame(suppressWarnings(rxSolve(mb, ev, params = pm, atol = 1e-11, rtol = 1e-11)))
      for (st in c("depot", "central")) {
        fd <- (sp[[st]] - sm[[st]]) / (2 * hh)
        mx <- max(mx, max(abs(f[[sprintf("rx__sens_%s_BY_%s__", st, pn)]] - fd), na.rm = TRUE))
      }
    }
    # 1e-4 keeps power against the ~30% skip bug; 1e-5 flakes on macOS ARM64 (3.5e-5)
    expect_lt(mx, 1e-4)
  })

  test_that("modeled rate() moving-boundary forward jump-sensitivities match FD (non-ss AND ss)", {
    # The modeled rate() OFF boundary tau2 = tau1 + F*amt/rate(p) moves with the
    # parameters; that transversality jump [S] = amt*dF - (F*amt/rate)*d(rate)/dp
    # was previously DEFERRED (MODEL_RATE_OFF only removed the continuous forcing),
    # so both regular AND steady-state modeled-rate sensitivities were ~30% wrong.
    cs <- c("ka", "cl"); p <- c(ka = 1.2, cl = 3.5, v = 25)
    chk <- function(mt, ev, drop = NA) {
      st0 <- rxode2::.rxAdjointExpand(mt, cs)$st
      mb <- rxode2(mt); mj <- rxode2(mt, calcSens = cs, eventSens = "jump")
      a0 <- list(mj, ev, params = p, atol = 1e-11, rtol = 1e-11); if (!is.na(drop)) a0$addlDropSs <- drop
      f <- as.data.frame(suppressWarnings(do.call(rxSolve, a0)))
      mx <- 0
      for (pn in cs) {
        hh <- abs(p[[pn]]) * 1e-6; pp <- p; pm <- p; pp[pn] <- pp[pn] + hh; pm[pn] <- pm[pn] - hh
        a1 <- list(mb, ev, params = pp, atol = 1e-11, rtol = 1e-11); a2 <- list(mb, ev, params = pm, atol = 1e-11, rtol = 1e-11)
        if (!is.na(drop)) { a1$addlDropSs <- drop; a2$addlDropSs <- drop }
        sp <- as.data.frame(suppressWarnings(do.call(rxSolve, a1)))
        sm <- as.data.frame(suppressWarnings(do.call(rxSolve, a2)))
        for (st in st0) { fd <- (sp[[st]] - sm[[st]]) / (2 * hh)
          mx <- max(mx, max(abs(f[[sprintf("rx__sens_%s_BY_%s__", st, pn)]] - fd), na.rm = TRUE)) }
      }
      expect_lt(mx, 5e-5)   # FD-truncation-safe; the deferred-boundary bug was ~1e1
    }
    mt  <- "d/dt(depot)=-ka*depot\nd/dt(central)=ka*depot-(cl/v)*central\nrate(central)=11*cl/3.5\ncp=central/(v/1000)"
    mtF <- "d/dt(depot)=-ka*depot\nd/dt(central)=ka*depot-(cl/v)*central\nf(central)=0.61\nrate(central)=11*cl/3.5\ncp=central/(v/1000)"
    chk(mt,  et(amt = 100, rate = -1, cmt = "central") |> et(c(1, 4, 8, 12, 16, 20)))                              # non-ss
    chk(mtF, et(amt = 100, rate = -1, cmt = "central") |> et(c(1, 4, 8, 12, 16, 20)))                              # non-ss, F != 1
    chk(mt,  et(amt = 100, rate = -1, cmt = "central", ss = 1, ii = 24) |> et(c(4, 8, 12, 16, 24, 30)), TRUE)      # ss
  })
})