in R desolve SIR model with loop and if else

I have a simple SIR model, I am trying to implement a vaccination approach (V) where at first it is checked if the infected are higher than a threshold (100), and if there are still enough susceptibles left (50) it will vaccinate a certain number per timestep (50).

However what I am trying to do, is that once the condition is fulfilled it should vaccinate for 7 days (regardless of whether during these 7 days the infected are still higher than the threshold or not, for example if after day 4, I = 70 it should still continue, it should only stop if S < 50. After the 7 days are finished, it should check the condition again, and either start again for another 7 days or not.

What I have so far, I would be grateful if someone helped me implement that loop

sirV=function(time, y, params){

S    = y[1]
I    = y[2]
R    = y[3]
V   = y[4]

 with(as.list(params),{
 vac_helper = if (I > 100 & S > 50) {50}

 else {0}
 N = S+I+R+V
 dS = -S*beta*I/N  - vac_helper
 dI = S*beta*I/N - gamma*I
 dR = +gamma*I
 dV = vac_helper

 return(list(c(dS, dI, dR, dV)))
})
}


 myparameters = c(gamma=1/10,beta=0.2)
 times <- seq(0, 300)
 my_ode <- as.data.frame(ode( y=c(100000, 10, 0,0), times, sirV, myparameters))

Solution

Here is a suggestion, but I am not entirely sure it behaves the way you wanted to, so please check it first and come back to me if it doesn't. Note that end has to be in the global environment.

library(deSolve)

sirV = function(time, y, params){
  
  S    = y[1]
  I    = y[2]
  R    = y[3]
  V    = y[4]
  
  with(as.list(params),{
    
    # Has the previous vaccination ended, and do we still need to vaccinate?
    if(time > end & I > 100) {
      end <<- time + 7
    }
    
    # Can we vaccinate?
    vaccinate = ifelse(end > time & S > 50, 50, 0)
    
    N = S+I+R+V
    dS = -S*beta*I/N  - vaccinate
    dI = S*beta*I/N - gamma*I
    dR = gamma*I
    dV = vaccinate
    
    # Store some results in a global list so we can check what is happening under the hood
    temp = data.frame(time, end, S, I, R, V, dS, dI, dR, dV, vaccinate)
    catch_results[[length(catch_results)+1]] <<- temp

    return(list(c(dS, dI, dR, dV)))
  })
}

myparameters = c(gamma = 1/10, beta = 0.2)
times <- seq(from = 0, to = 300, by = 1)


end <- 0 # Reset end time
catch_results = list() # Catch results from inside the function
my_ode <- ode( y=c(100000, 10, 0, 0), times, sirV, myparameters)

plot(my_ode)

# Check this to see if we get the expected behavior, especially at around time = 189
results = dplyr::bind_rows(catch_results)

^{Created on 2021-08-22 by the reprex package (v0.3.0)}