tests/testthat/test-state-dep-rate.R

In rxode2: Facilities for Simulating from ODE-Based Models

rxTest({
  # -- Analytical helper ------------------------------------------------------
  # depot(t) for a single constant-rate infusion into a one-compartment depot
  # that drains with first-order rate constant ka.
  #
  #   R  = infusion rate (amount / time)
  #   D  = duration      = amt / R   (infusion starts at t = 0)
  #   ka = depot elimination rate constant
  #
  # During infusion (0 <= t <= D):
  #   depot(t) = (R/ka) * (1 - exp(-ka*t))
  #
  # After infusion (t > D):
  #   depot(t) = depot(D) * exp(-ka*(t-D))
  #            = (R/ka) * (1-exp(-ka*D)) * exp(-ka*(t-D))
  .depotAnalytic <- function(t, R, D, ka) {
    vapply(t, function(ti) {
      if (ti <= 0) return(0)
      if (ti <= D) {
        (R / ka) * (1 - exp(-ka * ti))
      } else {
        depot_D <- (R / ka) * (1 - exp(-ka * D))
        depot_D * exp(-ka * (ti - D))
      }
    }, numeric(1))
  }

  # -- Test 1 -----------------------------------------------------------------
  # state = 0 throughout  ->  rate(depot) = rate0 + 0 = rate0 exactly.
  # This is the baseline: no state dependence, behaviour identical to a fixed
  # rate infusion.
  test_that("state-dep rate with state=0 uses exactly rate0 (no state shift)", {
    mod <- rxode2({
      # Infusion rate is modelled as rate0 + state.
      # With state held at 0, this collapses to the plain fixed-rate case.
      rate(depot) <- rate0 + state
      d/dt(depot) <- -ka * depot
      d/dt(state) <- 0         # constant state (initial value is the only value)
    })

    amt   <- 100
    rate0 <- 10
    ka    <- 0.5
    R     <- rate0            # state=0 -> effective rate = rate0
    D     <- amt / R          # expected infusion duration
    times <- seq(0, 30, by = 0.5)

    et <- eventTable() |>
      add.dosing(dose = amt, rate = -1) |>  # rate=-1 -> use modeled rate
      add.sampling(times)

    s <- solve(mod, et,
               inits  = c(state = 0, depot = 0),
               params = c(rate0 = rate0, ka = ka))

    expected <- .depotAnalytic(times, R, D, ka)

    # Every sampled depot value must match the analytical trajectory computed
    # with rate = rate0.  A mismatch would mean the solver used a different rate.
    expect_equal(s$depot[s$time %in% times], expected,
                 tolerance = 1e-3,
                 label = "depot matches rate0-only analytical solution when state=0")
  })

  # -- Test 2 -----------------------------------------------------------------
  # state = S0 (constant, positive and negative) throughout.
  # The effective rate at infusion start is rate0 + S0; duration = amt/(rate0+S0).
  # The full depot trajectory must match the analytical formula for that rate.
  test_that("state-dep rate shifts infusion rate by state value: rate = rate0 + state_at_start", {
    mod <- rxode2({
      # Infusion rate includes a state-variable contribution.
      # With d/dt(state)=0 the state stays at its initial value for the whole
      # simulation, giving an analytically verifiable rate.
      rate(depot) <- rate0 + state
      d/dt(depot) <- -ka * depot
      d/dt(state) <- 0
    })

    amt   <- 100
    rate0 <- 10
    ka    <- 0.5
    times <- seq(0, 40, by = 0.5)

    for (S0 in c(2, 5, -3)) {
      R <- rate0 + S0         # effective rate at the moment the infusion starts
      D <- amt / R            # duration changes with state value

      et <- eventTable() |>
        add.dosing(dose = amt, rate = -1) |>
        add.sampling(times)

      s <- solve(mod, et,
                 inits  = c(state = S0, depot = 0),
                 params = c(rate0 = rate0, ka = ka))

      expected <- .depotAnalytic(times, R, D, ka)

      # The solver must use rate = rate0 + S0, not rate0 alone.
      # Verifying the full trajectory confirms both the rate magnitude and the
      # derived infusion duration (D = amt/R).
      expect_equal(s$depot[s$time %in% times], expected,
                   tolerance = 1e-3,
                   label = paste0("depot matches analytic solution for rate=",
                                  R, " (S0=", S0, ")"))

      # Confirm the trajectory with S0!=0 differs from the S0=0 baseline,
      # i.e., state truly shifts the rate and is not silently ignored.
      if (S0 != 0) {
        expected_s0 <- .depotAnalytic(times, rate0, amt / rate0, ka)
        expect_false(isTRUE(all.equal(expected, expected_s0)),
                     label = paste0("state S0=", S0,
                                    " produces a different trajectory than state=0"))
      }
    }
  })

  # -- Test 3 -----------------------------------------------------------------
  # The infusion end time must be exactly amt / (rate0 + state_at_start).
  # We probe this by checking depot at the expected stop time D and well after.
  # During the infusion depot is rising; after it decays.  Confirming depot(D)
  # matches the analytical peak locks in the exact stop time.
  test_that("state-dep rate: infusion ends at exactly amt/(rate0 + state_at_start)", {
    mod <- rxode2({
      rate(depot) <- rate0 + state
      d/dt(depot) <- -ka * depot
      d/dt(state) <- 0
    })

    amt   <- 100
    rate0 <- 5
    ka    <- 0.3

    for (S0 in c(0, 3, -1)) {
      R <- rate0 + S0
      D <- amt / R           # expected infusion end time

      # Sample finely around D so we capture exactly the stop-event time
      times <- sort(unique(c(seq(0, D + 8, by = 0.1), D)))

      et <- eventTable() |>
        add.dosing(dose = amt, rate = -1) |>
        add.sampling(times)

      s <- solve(mod, et,
                 inits  = c(state = S0, depot = 0),
                 params = c(rate0 = rate0, ka = ka))

      # depot(D) from the analytical formula
      depot_at_D_expected <- (R / ka) * (1 - exp(-ka * D))

      # observed depot at the exact stop-event sample time
      obs_at_D <- s$depot[which.min(abs(s$time - D))]

      # If rate was wrong the infusion would stop at a different time, yielding
      # a depot peak that is too large (rate too high, shorter duration) or too
      # small (rate too low, longer duration).
      expect_equal(obs_at_D, depot_at_D_expected,
                   tolerance = 1e-3,
                   label = paste0("depot at infusion stop D=", round(D, 4),
                                  " matches analytic (S0=", S0, ")"))
    }
  })

  # -- Test 4 -----------------------------------------------------------------
  # Multiple doses: each dose's rate is evaluated at the state present when
  # THAT dose's infusion begins.  Because state changes over time (exponential
  # decay here), consecutive doses get different effective rates and therefore
  # different durations.
  # We verify no NAs and that the two-dose trajectory is self-consistent (not
  # all identical) when the state is non-zero.
  test_that("state-dep rate: consecutive doses evaluated at their own start state", {
    mod <- rxode2({
      # rate depends on a decaying state variable so each dose sees a different
      # rate if the infusions are spaced far enough apart.
      rate(depot) <- rate0 + state
      d/dt(depot) <- -ka * depot
      d/dt(state) <- -kstate * state   # state decays between doses
    })

    et <- et() |>
      et(amt = 50, rate = -1, time = 0) |>    # dose 1 - state is large
      et(amt = 50, rate = -1, time = 15) |>   # dose 2 - state has decayed
      et(seq(0, 30, by = 0.5))

    s <- rxSolve(mod, et,
                 params = c(rate0 = 5, ka = 0.4, kstate = 0.2),
                 inits  = c(depot = 0, state = 4))

    # No numerical failures
    expect_false(any(is.na(s$depot)),
                 label = "no NAs in depot with two state-dep rate doses")

    # The trajectory around dose 2 must differ from what a state=0 run would
    # produce (state has not decayed to exactly 0 by t=15).
    s0 <- rxSolve(mod, et,
                  params = c(rate0 = 5, ka = 0.4, kstate = 0.2),
                  inits  = c(depot = 0, state = 0))

    expect_false(isTRUE(all.equal(s$depot, s0$depot)),
                 label = "nonzero initial state gives different depot trajectory than state=0")
  })
})