Plot vector field of SIR model

SIRmod.R

library(deSolve)
library(ggplot2)
library(tidyr)
library(dplyr)
library(foreach)
SIRmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS <- - beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I 
    return(list(c(dS,  dI, dR)))
  })
}

pars  <- c(beta=2, gamma=1)
yini  <- c(S=0.5, I=0.01, R=0)
times <- seq(0, 50, by = 0.1)
out   <- ode(yini, times, SIRmod, pars)

df_out <- pivot_longer(data.frame(out), S:R) %>% 
  mutate(name=factor(name, levels = c("S","I","R")))

ggplot(df_out, aes(x=time, y=value, colour=name, linetype = name))+
  geom_line()+
  theme_bw()

sol_ode_from_states <- function(df_state, times, func, parms){
  res_df <- foreach(i=1:nrow(df_state), .packages='deSolve', .combine=rbind) %do% {
    ode_out <- ode(y=c(unlist(df_state[i,])), times=times, func=func, parms=parms)
    data.frame(group=i, ode_out)
  }
  return(res_df)
}

df_state <- expand.grid(S=seq(.01, 1, 0.05), I=seq(.01, 1, 0.05)) %>% 
  mutate(R=0)
times <- seq(0, 0.1, by = 0.01)
res <- sol_ode_from_states(df_state=df_state, times=times, func=SIRmod, parms=pars)

ggplot(data=res, aes(x=S, y=I, group=group, color=time)) +
  geom_path(arrow = arrow(length = unit(3, "pt"))) +
  scale_colour_gradient2(low="grey80", high="firebrick")+
  geom_vline(xintercept = pars["gamma"]/pars["beta"], linetype=2)+
  theme_bw()
ggsave("SIRmod.png", width = 5, height = 5)