# This is the R version of program 2.5 from page 39 of 
# "Modeling Infectious Disease in humans and animals" 
# by Keeling & Rohani.
 
# It is the simple SIS epidemic without births or deaths.

# Set up equations as a function
SIS = function(time,state,pars){
  with(as.list(c(state,pars)),{
    dSdt =  -beta*S*I + gamma*I
    dIdt = beta*S*I - gamma*I
    return(list(c(dSdt,dIdt)))
  })
}

# Set parameters (all rates are per day)
beta=1.4247         # Transmission rate 
gamma=0.14286       # Recovery rate
pars = c(beta,gamma)

# Initial conditions
I0 = 1e-6       # initial proportion of population that are infectious
S0 = 1-I0       # initial proportion of population that are susceptible
y0 = c(S=S0,I=I0)

# Time range
MaxTime=70
step=0.1
times = seq(0,MaxTime,step)

# Check all parameters are valid
if(S0<=0){
  print(paste("Initial level of susceptibles is less than or equal to zero, S0 =",S0))  
}
if(I0<=0){
  print(paste("Initial level of infecteds is less than or equal to zero, I0 =",I0))  
}
if(beta<=0){
  print(paste("Transmission rate is less than or equal to zero, beta =",beta))  
}
if(gamma<=0){
  print(paste("Recovery rate is less than or equal to zero, gamma =",gamma))  
}
if(MaxTime<=0){
  print(paste("Maximum run time is less than or equal to zero, MaxTime =",MaxTime))  
}
if(beta<gamma){
  print(paste("Basic reproductive ratio is less than 1, R =",beta/gamma))
}

# Solve differential equations
library(deSolve)
out = as.data.frame(ode(y0,times,SIS, pars))

# Plot results using ggplot2
library(ggplot2)

p1=ggplot() +
  geom_line(data=out,aes(y=S,x= time),color="green3",linewidth=0.6,linetype=1) +
  xlab("Time") + ylab("Susceptibles") +
  scale_y_continuous(limits=c(0,1),expand=expansion(mult=c(0,0.05))) +
  scale_x_continuous(limits=c(0,70),expand=expansion(mult=c(0.0,0.0))) +
  theme_classic() +
  theme(legend.position = "none") +
  theme(panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 0.6))

p2=ggplot() +
  geom_line(data=out,aes(x=time,y=I),color="red",linewidth=0.6,linetype=1) +
  xlab("Time") + ylab("Infectious") +
  scale_y_continuous(limits=c(0,1),expand=expansion(mult=c(0,0.05))) +
  scale_x_continuous(expand=expansion(mult=c(0,0))) +
  theme_classic() +
  theme(legend.position="none") +
  theme(panel.border = element_rect(color = "black", 
                                    fill = NA, linewidth = 0.6))

library(patchwork)
p1/p2