Apply event/root function to large set of equations R deSolve

Question

TLDR: My main challenge is just how do I write the root function that checks an arbitrary number of state variables, x, and then apply the event function such that all state variables n that have a value less than the threshold (n <= x) are acted upon by the event function?

I'm trying to use deSolve for a set of Lotka-Volterra equations, but with many state variables (i.e. not just a predator and prey but 20 interacting organisms).

I want to use a root function and event function to be constantly checking if any state variable values dip below a threshold value (e.g. 1.0) and if they do, use the event function to make that particular state variable 0. I've been messing around with a simple minimal example, but can't quite understand how to extend this to check all the state variables and just apply to the one(s) that is/are below the threshold.

The LV example from the deSolve package vignette

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
    
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
    
    return(list(c(dPrey, dPredator)))
  })
}

pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity

yini  <- c(Prey = 10, Predator = 2)
times <- seq(0, 50, by = 1)

I can apply my root and event functions to check for just the prey's values:

## event triggered if state variable less than 1
rootfun <- function (Time, State, Pars) {
  return(State[1] - 1)
}

## sets state variable = 1                                                  
eventfun <- function(Time, State, Pars) {
  return(c(State[1] <- 0, State[2]))
}
out   <- lsode(yini, times, LVmod, pars, 
               rootfunc = rootfun, 
               events = list(func = eventfun, root = TRUE))


## User specified plotting
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "time", ylab = "Conc",
        main = "Lotka-Volterra", lwd = 2)
legend("topright", c("prey", "predator"), col = 1:2, lty = 1:2)

And the result is this:

But now I want to extend this so that it checks all the state variables (in this case just the 2), but ideally in a way that is flexible to different numbers of state variables. I have tried messing around with doing this in some sort of loop but can't seem to figure it out. My main challenge is just how do I write the root function that checks an arbitrary number of state variables, x, and then apply the event function such that all state variables n that have a value less than the threshold (n <= x) are acted upon by the event function?

Perhaps useful information is at some point I would like to implement a separate (not root-based) event function to change a parameter at some pre-set times, so ideally the solution to this problem could interface with additional event function implementation.

Help much appreciated as always!!

Answer 1

One can use a vectorized version of the LV model and then write rootfun and eventfun also in vectorized style:

library(deSolve)

model <- function(t, y, parms) {
  with(parms, {
    dy <- r * y  + y * (A %*% y)
    list(dy)
  })
}

## int6eraction matrix
parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  A = matrix(c(
    0.0, 0.0, -0.2, 0.0, # prey 1
    0.0, 0.0, 0.0, -0.1, # prey 2
    0.2, 0.0, 0.0, 0.0,  # predator 1; eats prey 1
    0.0, 0.1, 0.0, 0.0), # predator 2; eats prey 2
    nrow = 4, ncol = 4, byrow = TRUE)
)

times = seq(0, 150, 1)
y0  = c(n1 = 1, n2 = 1, n3 = 2, n4 = 2)

out <- ode(y0, times, model, parms)
plot(out)

## defined as global variables for simplicity, can also be put into parms
threshold <- 0.2  # can be a vector of length(y0)
y_new     <- 1.0  # can be a vector of length(y0)

## uncomment the 'cat' lines to see what's going on
rootfun <- function (t, y, p) {
  #cat("root at t=", t, "\n")
  #cat("y old =", y, "\n")
  return(y - threshold)
}

eventfun <- function(t, y, p) {
  #cat("y old =", y, "\n")
  y <- ifelse(y <= threshold, y_new, y)
  #cat("y new =", y, "\n")
  return(y)
}

out <- ode(y0, times, model, parms, 
           events = list(func = eventfun, root = TRUE), rootfunc=rootfun)
plot(out)