In R, SIR model,not constant parameters

Question

On this page a SIR model in R is shown, https://rstudio-pubs-static.s3.amazonaws.com/382648_93783f69a2fd4df98ade8751c21abbad.html, the solution of it and the optimization of the $\beta$ and $\gamma$ parameter is also executed. (see below)

In this code both $\beta$ and $\gamma$ are assumed to be constant over the whole time. What I want is to to have a time varying beta, it does not need to change each day, we have fourteen days of data, it would suffice if it would change after seven days, i.e we have $\beta_1$ for days[0:6] and $\beta_2$ for days[7:13] and then do the optimization algorithm like below for both, i.e. in the end I want to receive a vector for the optimal values of (\beta_1, \beta_2, \gamma) whereas gamma stayed constant the whole time. Would it be possible with a modification of the code given? If yes could someone help how to modify it to receive the desired output.

   day  cases
   0    1
   1    6
   2    26
   3    73
   4    222
   5    293
   6    258
   7    236
   8    191
   9    124
   10   69
  11    26
  12    11
  13    4

 #here beta is assumed to be constant
sir_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
dS <- -beta * I * S
dI <-  beta * I * S - gamma * I
dR <-  gamma * I
return(list(c(dS, dI, dR)))
})
}

parameters_values <- c(
beta  = 0.004, # infectious contact rate (/person/day)
gamma = 0.5    # recovery rate (/day)
)

initial_values <- c(
S = 999,  # number of susceptibles at time = 0
I =   1,  # number of infectious at time = 0
R =   0   # number of recovered (and immune) at time = 0
)

time_values <- seq(0, 10) # days


sir_values_1 <- ode(
y = initial_values,
times = time_values,
func = sir_equations,
parms = parameters_values 
)

sir_values_1

sir_values_1 <- as.data.frame(sir_values_1)
sir_values_1

sir_1 <- function(beta, gamma, S0, I0, R0, times) {
require(deSolve) # for the "ode" function

# the differential equations:
sir_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
  dS <- -beta * I * S
  dI <-  beta * I * S - gamma * I
  dR <-  gamma * I
  return(list(c(dS, dI, dR)))
  })
  }

# the parameters values:
parameters_values <- c(beta  = beta, gamma = gamma)

# the initial values of variables:
initial_values <- c(S = S0, I = I0, R = R0)

# solving
out <- ode(initial_values, times, sir_equations, parameters_values)

# returning the output:
as.data.frame(out)
}
sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = seq(0, 10))

flu <- read.table("https://uc8f29367cc06ca2f989ead2cd8e.dl.dropboxusercontent.com/cd/0/inline/BNzBF_deK5fmfGXWCB9a5YO95JkiLNFRc2Jq1w-qGNqQMXxnpn-yL-cAVoE1JQG7D4Od_SkG8YVKesqBr7wMoQHHSTNbHU_hhyahK7up0EDEft-u7Vf4xZJvu4cTNuUjXFb-QaHlOfBPnFhKspeb7RbO/file", header = TRUE)

predictions <- sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = flu$day)
predictions

model_fit <- function(beta, gamma, data, N = 763, ...) {
I0 <- data$cases[1] # initial number of infected (from data)
times <- data$day   # time points (from data)
# model's predictions:
predictions <- sir_1(beta = beta, gamma = gamma,   # parameters
                   S0 = N - I0, I0 = I0, R0 = 0, # variables' intial values
                   times = times)                # time points
# plotting the observed prevalences:
 with(data, plot(day, cases, ...))
# adding the model-predicted prevalence:
with(predictions, lines(time, I, col = "red"))
}

predictions <- sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = flu$day)
predictions

ss <- function(beta, gamma, data = flu, N = 763) {
I0 <- data$cases[1]
times <- data$day
 predictions <- sir_1(beta = beta, gamma = gamma,   # parameters
                   S0 = N - I0, I0 = I0, R0 = 0, # variables' intial values
                   times = times)                # time points
sum((predictions$I[-1] - data$cases[-1])^2)
}
ss(beta = 0.004, gamma = 0.5)

beta_val <- seq(from = 0.0016, to = 0.004, le = 100)
ss_val <- sapply(beta_val, ss, gamma = 0.5)

min_ss_val <- min(ss_val)
min_ss_val

beta_hat <- beta_val[ss_val == min_ss_val]
beta_hat

 plot(beta_val, ss_val, type = "l", lwd = 2,
 xlab = expression(paste("infectious contact rate ", beta)),
 ylab = "sum of squares")
# adding the minimal value of the sum of squares:
abline(h = min_ss_val, lty = 2, col = "grey")
# adding the estimate of beta:
abline(v = beta_hat, lty = 2, col = "grey")

ss(beta = 0.004, gamma = 0.5)

ss2 <- function(x) {
ss(beta = x[1], gamma = x[2])
}
ss2(c(0.004, 0.5))



starting_param_val <- c(0.004, 0.5)
ss_optim <- optim(starting_param_val, ss2)

Answer 1

This is certainly possible. All you need is an if statement in your gradient function:

beta <- if (time<6) beta1 else beta2

or

beta <- ifelse(time<6, beta1, beta2))

and make sure your parameter vector includes both beta1 and beta2.