I am trying to compare an SEIR model with no vaccination to one with vaccination, as you can see from my code:

Initial_infect = 1
R0 = 2.5
gamma = 1/7  # recovery rate
beta = R0 * gamma   # contact rate
sigma = 0.2  # infectious rate
N = 1e6


model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta * S * I / N
    dE <- beta * S * I / N - sigma * E
    dI <- sigma * E - gamma * I
    dR <- gamma * I
    return(list(c(dS, dE, dI, dR)))
  })
}


run_model <- function(vax, times) {
  initial_state <- c(
    S = (1 - vax) * N,
    E = 0,
    I = Initial_infect,
    R = vax * N
  )
  parms <- c(beta = beta, sigma = sigma, gamma = gamma, N = N)
  out <- ode(y = initial_state, times = times, func = model, parms = parms)
  return(out)
}

times <- seq(0, 200, by = 1)
results_no_vax <- run_model(0, times)

plot(S ~ times, data=results_no_vax, type="l", ylim=c(0,1e6), main = "SEIR model no vaccination", xlab="Time", ylab="Population")
lines(E~times, data=results_no_vax, col="red")
lines(I~times, data=results_no_vax, col="green")
lines(R~times, data=results_no_vax, col="blue")
legend("left", legend=c("Susceptible", "Exposed", "Infectous", "Reduced"), col = c("black", "red", "green", "blue"), lty=1)

vax = seq(0,1, by=0.1)
results_with_vax = lapply(vax, run_model, times)
plot(S ~ times, data=results_with_vax, type="l", ylim=c(0,1e6), main = "SEIR model with vaccination", xlab="Time", ylab="Population")
lines(E~times, data=results_with_vax, col="red")
lines(I~times, data=results_with_vax, col="green")
lines(R~times, data=results_with_vax, col="blue")
legend("left", legend=c("Susceptible", "Exposed", "Infectous", "Reduced"), col = c("black", "red", "green", "blue"), lty=1)

My problem being that when displaying the plot with the vaccination, it shows no difference to the one without the vaccination vector. So my question is if I did something wrong, or how am I able to pick a specific value from the vax vector (for example 0.5), to display on the plot without altering the run_model function and change the vax variable to be 0.5.

I so far have tried to have 2 different run_model functions (one for no vax, and one for vax), the no vax one had no vax variable, but the only thing that changed was shifting the S curve up, which doesn’t make much sense.