SI model overview

1. SI model function

This is our ODE model function.

• Line 1: tells R we are creating a function called SI_model and that this function takes 3 inputs: variables called times, yinit, and pars:

• Line 2: pars & ynit: names and values for all of our parameters and initial model conditions, respectively.
• Line 3 & 4: differential equations for an SI model. In R you can use “<-” and “=” essentially interchangeably.
• Line 5: tells the function to return the values for the derivative of the number of susceptibles and infecteds and then close out the function. The output derivatives MUST be in the same order as the initial conditions (here, that means susceptibles must come before infecteds for both).

SI_model <- function(times,yinit,pars){ # input
  with(as.list(c(yinit,pars)), {
    dS <- b - beta*I*S - mu_s*S
    dI <- beta*I*S - mu_i*I
    return(list(c(dS, dI)))})       # output 
}

2. Setting Parameters.

We are defining the values of the parameters used in our model.

• beta : transmission parameter
  
• mu_s: mortality rate for susceptibles.
  
• mu_i: mortality rate for infecteds.
• b: birth rate (assumed not to depend on population size).
• pars: last line combines all of these values into a vector named pars that retains both the names and the values of each parameter

This way we can refer to all of the parameters at once by referring to the variable pars. If you type pars into the command line after running all of the code, it will output the values you have assigned to each parameter.

beta = 8e-6 # beta is transmission parameter
mu_s = 1/30 # mu_s is the mortality rate for susceptibles
mu_i = 1/10 # mu_i is the mortality rate for infecteds
b = 1000 # b is the birth rate (assumed not to depend on population size)
pars <- cbind(b,beta,mu_s,mu_i)

3. Setting Time Frame

Here, we are setting the time frame of interest. I’ve decided to let my start time be 0 and my end time be 200 years. I then ask R for a vector ranging from year 0 to year 200 in increments of 1 and name this times.

start_time = 0 # start date
end_time = 200 # end date 
times <- seq(start_time, end_time, by = 1) # gives a sequence from start to end

4. Setting Initial Conditions

Here I am telling R now many susceptibles (30000) and how many infecteds (10) are in my population at time 0 when I want the model to start.

# initial starting values
yinit <- c(S = 30000, I = 10)

5. Running The Model

This line of code is where R is doing all the work!

To breakdown the code: We are using the function ode to numerically integrate our model. We are telling R the names we have given to our initial conditions (y=yinit), the times we’re interested in (times=times), the parameters (parms=pars), and the function that contains our ODE system (func=SI_model). We’re also telling it to use a Runge-Kutta integration method (method=’rk4’). Finally, as.data.frame is simply a built-in R function that tells R how to format the results.

results <- as.data.frame(ode(y=yinit,times=times,func=SI_model,parms=pars,method='rk4'))

6. Plotting Results

This final bit of code simply deals with plotting the results in a way that looks nice.

colors <- c('Susceptible'='black','Infected'='red', 'Total'='blue')
SIplot <- ggplot(data=results) + 
  geom_line(aes(x=time,y=S,col='Susceptible')) + 
  geom_line(aes(x=time,y=I,col='Infected')) + 
  geom_line(aes(x=time,y=(I+S),col='Total')) +
  scale_colour_manual(name="Legend",values=colors) +
  xlab('Year') + ylab('Number of People')+
  theme_minimal()

SIplot