sir_sim.R

library(ggplot2)
library(deSolve)
library(tidyr)
theme_set(theme_minimal())

# I: infected, S: susceptible, R: recovered
sir <- function(t, state, pars) {
  with(as.list(c(state, pars)), {
    N <- S + I + R
    dSdt <- -b * I/N * S
    dIdt <-  b * I/N * S - g * I
    dRdt <-  g * I
    return(list(c(dSdt, dIdt, dRdt)))
  })
}

pars <- c(b = 2, g = 0.5)
state <- c(S = 762, I = 1, R = 0)
times <- seq(1, 15, length = 100)
sol <-
  ode(
    y = state,
    times = times,
    func = sir,
    parms = pars,
    method = "ode45"
  )
sol <- as.data.frame(sol)
sol <- sol %>% 
  tidyr::pivot_longer(-time, names_to = "state", values_to = "value")

sol %>%
  ggplot(aes(time, value, color = state)) +
  geom_point(size = 0.5) + 
  geom_line() +
  guides(color = guide_legend(title = NULL)) +
  scale_color_discrete(labels = c("Infected","Recovered","Susceptible")) +
  xlab("Days") + ylab("Number of people") +
  ggtitle("Basic SIR model", 
          subtitle = "3-state ODE with beta = 2, gamma = 0.5, R0 = 4")