Making an SEIR model parameter to depend on incidence of a compartment at each time step using ode function

Below is the code in my attempt to make parameter r depend on the incidence of I in a SEIR model. The model runs if (c(0, diff(I)))/10000 is replaced with any value less than one. However, I would like the parameter r to depend on the value of I minus the value I at the previous time step at each time step.

library(deSolve)

seirs <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    dS= -lambda*S + rho*R
    dE= lambda*S - gamma_h*E
    dI=gamma_h*E - r*I
    dR=r*I - rho*R
    dCInc = lambda*S
    output <- c(dS, dE, dI, dR, dCInc)
    list(output)
  })
}

start<-c(S=5000, E=3000, I=2000,R=0, CInc = 0 )

parms <- c(lambda = 0.3, gamma_h = 1/20, r= (c(0, diff(I)))/10000, rho = 1/(1*365))
times <- seq(0, 10000, 1)
result <-ode(times=times, y=start, func=seirs,parms=parms)

A parameter is usually a fixed value. A time-dependent value can be implemented by several methods:

1. Integrate the dependency into the model function, e.g. by adding an additional differential equation for the parameter,
2. Use a forcing function and interpolate from an external "signal",
3. Use events, to modify the state variables. This method can be combined with method #1.

The methods #2 and #3 are described at the ?forcings and ?events help pages of the deSolve package.

To make a variable (or parameter) dependent on the difference of a state variable of the current and the previous time step, it first needs to define, what is meant by a "time step" in a differential-equation system. Differential equations are continuous by definition, so in theory the "time step" is infinitesimally small. In a numerical simulation, it has a fixed small value, that depends on the numerical solver. It is in most cases automatically adjusted and can vary.

So if "time step" means "very small", then the difference between I and the previous I, divided by the time step is just the derivative dI, see method #1 above. This is in my opinion the preferred approach

It can also be, that "time step" should mean a discrete time step, e.g. 1 day, assuming times is in days. In this case, one can use a delay differential equation (DDE). In deSolve, it can be implemented using the lagvalue-function and the dede-solver.

The following code contains an example for such a DDE. To get a plausible picture, the parameter values, the state variables and the times vector were changed. I also changed the difference from I - I_tau to I_tau - I to get positive values and added r as "auxiliary variables" to the output.

Please take this as a technical answer, because the motivation of the original question is still unclear to me, i.e. which biological process is should be modeled with such an approach. Possibly, method #1 from above may be a better option.

library(deSolve)

seirs <- function(t, x, parms)  {

  with(as.list(c(parms, x)), {
    if (t < tau)
      I_tau <- I # default value can also be set manually
    else
      I_tau <- lagvalue(t - tau, 1)

    r      <- k_r * (I_tau - I)
      
    dS     <- -lambda*S + rho*R
    dE     <- lambda*S - gamma_h*E
    dI     <- gamma_h*E - r*I
    dR     <- r*I - rho*R
    dCInc  <- lambda*S
    
    list(c(dS, dE, dI, dR, dCInc), r = r)
  })
}

start <- c(S=5000, E=3000, I=2000, R=0, CInc = 0)

parms <- c(lambda = 0.3, gamma_h = 1/20, k_r= 0.01, rho = 0.01, tau = 1)

times <- seq(0, 10, 1)
result <-dede(times=times, y=start, func=seirs, parms=parms)

plot(result)