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
           newbeta = beta*(1+ beta1*cos(omega*current_timepoint))
           dS = (birthrate) - (newbeta * S * I) - (deathrate*S)
           dI = ( newbeta * S * I) - (gamma * I) - (deathrate*I)
           dR = (gamma * I) - (deathrate*R)
           
           # combine results
           results = c (dS, dI, dR)
           list (results)
         } 
        )
}

Parameters

contact_rate = 15            # number of contacts per day
transmission_probability = 0.03    # transmission probability
infectious_period = 5              # infectious period
birthrate = 1/250                   # birthdate
deathrate = 1/250                   # deathrate

Compute values of beta,gammae,deathrate,birthrate.

beta_value = contact_rate * transmission_probability
beta1 = amplitude_of_seasonality = 0.5
period_of_seasonality = 365*10             # day
pi = 22/7
omega = 2*pi / period_of_seasonality 
gamma_value = 1 / infectious_period
birthrate = 1/250                  
deathrate = 1/250

Compute Ro - Reproductive number.

Ro = beta_value / gamma_value

Disease dynamics parameters.

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

Initial values for sub-populations.

X = 9800       # susceptible hosts
Y = 500           # infectious hosts
Z = 1500           # 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, 365*20, 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 = 'Amplitude=0.5, period of seasonality=365*10') 
text(800,0.9, "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(2000,0.1, "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(1500,0.20, "R", col = 'green')