seiv_params <- list(
    ave_lifespan = 5*365, # days
    ave_rabies_duration = 10, # days
    beta = (2 / 7) * 0.5, # potential infectious interactions per dog per week
    N = 100000,
    ave_incubation_period = 6, # days
    ave_vaccination_period = 365, # days
    ave_vaccination_duration = 3 * 365 # days
)

seiv_pre_vax <- within(seiv_params, {
    ave_vaccination_period <- Inf
})

dogs_zero <- c(S = seiv_params$N - 1, E = 0, I = 1, V = 0)

seiv <- function(t, y, parms){
  with(c(as.list(y), parms), {
    death <- y * (1/ave_lifespan)
    disease_death <- I * (1/ave_rabies_duration)
    birth <- sum(death) + disease_death

    infection <- beta * S * I / N
    incubation <- 1 / (ave_incubation_period) * E
    vaccination <- 1 / (ave_vaccination_period) * S
    immunity_loss <- 1 / (ave_vaccination_duration) * V

    dSdt <- (birth + immunity_loss) - (infection + vaccination)
    dEdt <- (infection) - (incubation)
    dIdt <- (incubation) - (disease_death)
    dVdt <- (vaccination) - (immunity_loss)
    # Note: Population size is constant, so don't need to specify dRdt
    return(list(c(dSdt, dEdt, dIdt, dVdt) - death))
  })
}

library(deSolve)                # Load libary to be used for numerical integration

time.out <- seq(0, 5*365, 0.1)     # INCLUDE THE APPROPRIATE VALUE TO GET A SEQUENCE
                               # FROM 0 to 365 BY STEPS OF 0.1


ts_seiv_pv <- data.frame(lsoda(
  y = dogs_zero,               # Initial conditions for population
  times = time.out,             # Timepoints for evaluation
  func = seiv,                   # Function to evaluate
  parms = seiv_pre_vax                # Vector of parameters
))

with(
  ts_seiv_pv,
  plot(time, I, type = "l", xlab = "Time in days", ylab = "Number infected", main = "Canine Rabies in Tanzania", xlim = c(0, 5*365), bty = "n")
)

dogs_pv_zero <- (tail(ts_seiv_pv, 1) |> unlist())[-1]

ts_seiv <- tail(data.frame(lsoda(
        y = dogs_pv_zero,               # Initial conditions for population
        times = time.out,             # Timepoints for evaluation
        func = seiv,                   # Function to evaluate
        parms = seiv_params                # Vector of parameters
    )), 100)
ts_seiv$frac_vac <- 1

for (frac_vac in seq(.90, .10, by = -.1)) {
    local_params <- within(seiv_params, ave_vaccination_period <- ave_vaccination_period / frac_vac)
    dt <- tail(data.frame(lsoda(
        y = dogs_pv_zero,               # Initial conditions for population
        times = time.out,             # Timepoints for evaluation
        func = seiv,                   # Function to evaluate
        parms = local_params                # Vector of parameters
    )), 100)
    dt$frac_vac <- frac_vac
    ts_seiv <- rbind(ts_seiv, dt)
}



plot(ts_seiv$time,               # Time on the x axis
     ts_seiv$I,                  # Number infected (I) on the y axis
     xlab = "Time in days",     # Label the x axis
     ylab = "Number infected",  # Label the y axis
     main = "Canine Rabies in Tanzania",    # Plot title
     xlim = c(0, 5*365),           #
     type = "l",                # Use a line plot
     bty = "n")                 # Remove the box around the plot