Remove all objects from workspace.

remove (list = objects() )

Load add-on packages - deSolve - contains lsoda function - differential equation solver.

library (deSolve)
## 
## Attaching package: 'deSolve'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot

Function to compute derivatives of the differential equations.

sir_model = function (current_timepoint, state_values, parameters)
{  
  # create state variables (local variables)
  S = state_values [1]        # susceptibles
  I = state_values [2]        # infectious
  R = state_values [3]        # recovered
  
  with (
    as.list (parameters), # variable names within parameters can be used 
         { # compute derivatives
           dS = (-beta * S * I)
           dI = ( beta * S * I) - (gamma * I)
           dR = (gamma * I)
           
           # combine results
           results = c (dS, dI, dR)
           list (results)
         } 
        )
}

Parameters

contact_rate = 100                 # number of contacts per day
transmission_probability = 0.5    # transmission probability
infectious_period = 20            # infectious period

Compute values of beta (tranmission rate) and gamma (recovery rate).

beta_value = contact_rate * transmission_probability
gamma_value = 1 / infectious_period

Compute Ro - Reproductive number.

Ro = beta_value / gamma_value

Disease dynamics parameters.

parameter_list = c (beta = beta_value, gamma = gamma_value)

Initial values for sub-populations.

X = 980         # susceptible hosts
Y = 10          # infectious hosts
Z = 10           # recovered hosts

Compute total population.

N = X + Y + Z

Initial state values for the differential equations.

initial_values = c (S = X/N, I = Y/N, R = Z/N)

Output timepoints.

timepoints = seq (0, 30, by=1)

Simulate the SIR epidemic.

output = lsoda (initial_values, timepoints, sir_model, parameter_list)

Plot dynamics of Susceptibles sub-population.

plot (S ~ time, data = output, type='b', col = 'blue')

Plot dynamics of Infectious sub-population.

plot (I ~ time, data = output, type='b', col = 'red')

Plot dynamics of Recovered sub-population.

plot (R ~ time, data = output, type='b', col = 'green')

Plot dynamics of Susceptibles, Infectious and Recovered sub-populations in the same plot.

# susceptible hosts over time
plot (S ~ time, data = output, type='b', ylim = c(0,1), col = 'blue', ylab = 'S, I, R', main = 'SIR epidemic') 
text(15,0.1, "S", col = 'blue')

# remain on same frame
par (new = TRUE)    

# infectious hosts over time
plot (I ~ time, data = output, type='b', ylim = c(0,1), col = 'red', ylab = '', axes = FALSE) 
text(10,0.7, "I", col = 'red')

# remain on same frame
par (new = TRUE)  

# recovered hosts over time
plot (R ~ time, data = output, type='b', ylim = c(0,1), col = 'green', ylab = '', axes = FALSE)   
text(20,0.6, "R", col = 'green')