Find timepoint when a state variable meets a condition

Question

I have set up a model that looks like this

library(deSolve)

model <- function(t, x, parms) {
  S <- x[1]
  E <- x[2]
  I <- x[3]
  R <- x[4]
  K <- x[5]
  V <- x[6]
  #
  with(as.list(parms), {
    Q <- ifelse(t < day0 | t > day0 + duration, 0, 0.04)
    dS <- -B * S * I - Q * S
    dE <- B * S * I - r * E
    dI <- r * E - g * I
    dR <- g * I 
    dK <- r * E
    dV <- Q * S
    res <- c(dS, dE, dI, dR, dK, dV)
    list(res)
  })
}

mtime = 120

step_size = 0.2

pp <- c(day0 = 30, duration = 5, B = 0.04, r = 1/7, g = 1/7)

init = c(S = 99, E = 1, I = 0, R = 0, K = 0, V = 0)

output <- as.data.frame(lsoda(y = init, # Initial conditions for population
                              times = seq(0, mtime, step_size), # Timepoints for evaluation
                              func = model, # Function to evaluate
                              parms = pp)
                        )

I want to replace Q <- ifelse(t < day0 | t > day0 + duration, 0, 0.04) with Q <- ifelse(t < I_change_time + day0 | t > I_change_time + day0 + duration, 0, 0.04), where I_change_time = first timepoint at which I >= some value.

My question is in lsoda, how do I "catch" the first time at which a state variable changed? In this case, I just want to "catch" the first time I >= some value, so that I can use it to implement a process between two timepoints.

If this is not the best way to go about it, could you please suggest an alternative solution? Thank you.

Answer 1

After thinking a little, I think it is easier to avoid the event function and to split the simulation in two segments:

1. The first simulation is run to the occurrence of root.
2. The final time and the final state value(s) are used as initial time and initial state for the second segment.
3. The parameter vector is changed.
4. A simulation for the second segment is run without root finding .
5. Both segments are glued together.

The following code shows a toy example with one state variable. It can of course be extended to an arbitrary number of states. As a positive side effect, it needs no "if" constructions.

library("deSolve")

pharmaco <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dy <- a - b * blood
    list(dy)
  })
}

parms1 <- c(
  a = 0,   # constant dosage
  b = 0.5  # decay
)

yini1  <- c(blood = 10)
times1 <- seq(from = 0, to = 10, by = 0.1)

## rootfun is triggered when return value is zero
rootfun <- function(t, y, parms) {
  y["blood"] - 1
}

## solve the first segment until a root occurs
## e.g. blood concentration reaches 1.0
## important: the solver must support root finding, e.g.
## lsoda, adams or bdf, but not vode, rk4 or ode45
out1 <- ode(func = pharmaco, times = times1, y = yini1,
            parms = parms1, rootfun = rootfun, method = "lsoda")

## outcomment to inspect results
# plot(out1)
# tail(out1)

## start 2nd segment from last time step found by root function
## and set initial values to final state at this time
times2 <- seq(from = out1[nrow(out1), 1], to=10, by = 0.1)
yini2  <- out1[nrow(out1), 2]

## change parameters, e.g. set constant dosage to 1
parms2 <- c(a = 1, b = 0.5)

## run 2nd segment without root finding
out2 <- ode(func = pharmaco, times = times2, y = yini2, parms = parms2)

## glue output matrices together and assign class "deSolve"
## to make it compatible to the deSolve plot function
out <- rbind(out1, out2)
class(out) <- "deSolve"

## plot output and indicate root with dashed red line
plot(out)
abline(h = 1, col="red", lty="dashed")
abline(v = times2[1], col="red", lty="dashed")