R: Varying parameters based on states in ODE

Question

Please pardon me as this is my first time typing a post on SO. If you think I need to format my changes do let me know :)

I am currently trying to run an SIR model that introduces vaccination, and this vaccination will stop after hitting a threshold (e.g., 50% of the population), while the infections are ongoing.

In my compartmental model, I have the following states:

1. S
2. I
3. R
4. P
5. Sn
6. In
7. Rn

People from S will be vaccinated at a rate α, and the vaccination is either a success or a failure (based on some vaccine efficacy rate).

I understand that there are solutions on SO, but most of what I could find is stopping the parameter at some time, hence time-dependent parameter. However, I am interested in state-dependent parameter.

My current solution is to use the vaccination rate, α, is a state. And once the threshold has been reached (supposed to be captured by the root function), the event function will set α to 0.

I have tried using rootfunc and eventfunc, but I am honestly not sure if I am doing this correctly and hence the error.

Here, I have created the root and event function.

root_func <- function(t, states, parameters) {
    with(as.list(c(states, parameters)),
         {
             return(P + Rn - 0.5 * (S + I + R + P + Sn + In + Rn))  
         })
    
}

event_func <- function(t, states, parameters) {
    with(as.list(c(states, parameters)),
         {
             alpha <- 0
             return(list(c(S, I, R, P, Sn, In, Rn, alpha)))
         })
}

Next, I have created my SIR model.

model_1_take <- function(t, states, parameters) {
    with(
        as.list(c(states, parameters)),
        {
            N <- S + I + R + 
                P + 
                Sn + In + Rn
            
            alpha <- ifelse(t >= 100, 0, alpha)
            
            dS <- (mu * N) - 
                (beta * S * (I + In) / N) - 
                (alpha * epsilon * S) - 
                (alpha * (1 - epsilon) * S) - 
                (nu * S)
            dI <- (beta * S * (I + In) / N) - 
                (gamma * I) - 
                (nu * I)
            dR <- (gamma * I) - 
                (alpha * R) - 
                (nu * R)
            dP <- (alpha * epsilon * S) - 
                (nu * P)
            dSn <- (alpha * (1 - epsilon) * S) - 
                (beta * Sn * (I + In) / N) - 
                (nu * Sn)
            dIn <- (beta * Sn * (I + In) / N) - 
                (gamma * In) - 
                (nu * In)
            dRn <- (alpha * R) +
                (gamma * In) -
                (nu * Rn)
            
            return(list(c(dS, dI, dR, dP, dSn, dIn, dRn)))
        }
    )
}

Here, I am initializing my states. Notice that α is in the states, though it should technically be a parameter.

initN <- 5e6
initRn <- 0
initR <- 0
initIn <- 0
initI <- 1
initP <- 0
initSn <- 0
initS <- initN - initRn - initR - initI - initIn - initP - initSn

states_1_take <- c(
  S = initS, 
  I = initI, 
  R = initR,
  P = initP,
  Sn = initSn,
  In = initIn,
  Rn = initRn,
  alpha = 0.02
)

Lastly, I initialize my parameters.

R_0 <- 2
gamma <- 1 / 5
epsilon <- 0.75
rho <- 0.2
mu <- 1 / (365 * 70)
nu <- 1 / (365 * 70)
beta <- R_0 * (gamma + nu)

parameters_1_take <- c(
    beta = beta, 
    gamma = gamma,
    epsilon = epsilon,
    rho = rho, 
    mu = mu,
    nu = nu
)

When I want to run the model, I use the following code.

out <- ode(
        y = states_1_take,
        times = times, 
        func = model_1_take,
        parms = parameters_1_take,
        method = "lsodes",
    )

The error I have gotten is:

Error in checkEventFunc(Eventfunc, times, y, rho) : The number of values returned by events$func() (1) must equal the length of the initial conditions vector (8)

Would deeply appreciate any help!

Answer 1

To debug the problem, just call the model function with the default solver and without events to reproduce the error, e.g.:

times <- 0:100
out <- ode(
        y = states_1_take,
        times = times, 
        func = model_1_take,
        parms = parameters_1_take,
        method = "lsoda",
    )

It says that the model returns 7 derivatives while the initial vector contains 8 states.

One can also call the model function alone:

model_1_take(0, states_1_take, parameters_1_take)

The reason for the error is, that alpha is obviously no state variable. One way to solve this is adding alpha as a state that has normally a derivative dAlpha = 0 in the model function:

model_1_take <- function(t, states, parameters) {
  # ...
  dAlpha <- 0
            
  return(list(c(dS, dI, dR, dP, dSn, dIn, dRn, dAlpha)))
  # ...
}

It is then only changed in the event function.

Now we can run the model first without the events and then with events. Here we get another error because the event function should only return only the state vector, so it should be:

event_func <- function(t, states, parameters) {
  with(as.list(c(states, parameters)), {
    alpha <- 0
    c(S, I, R, P, Sn, In, Rn, alpha)
  })
}

times <- 0:100
out <- ode(
         y = states_1_take,
         times = times, 
         func = model_1_take,
         parms = parameters_1_take,
         method = "lsoda", 
         events = list(func=event_func, root = TRUE), 
         rootfun = root_func)

plot(out)

Two final notes:

• One can use the default lsoda solver here, that automatically switches between a stiff and a nonstiff solver has also root finding and event capabilities. The lsodes solver is only needed if the structure of a sparse Jacobian needs to be specified.
• The state variables have different orders of magnitude. This can sometimes lead to numerical problems of stability and accuracy. Here it can be wise to rescale the states to a "common range" between say 0.001 to 1000 or to assign individual absolute tolerances to the states: atol = c(.....)