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 = (birthrate) - (beta * S * I) - (deathrate*S)
           dI = ( beta * 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
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 = 98000       # 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, 900, 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 seasonality') 
text(130,0.3, "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(200,0.05, "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(150,0.65, "R", col = 'green')

Description: The SIR seasonal dynamic will change over the time period. At the beginning of period, susceptible is high, since susceptible expose to agents, number of susceptible decrease while number of infected increase. Number of recovered increase after host pass infectious period and got immunity. New born population will be filled in susceptible compartment over the time.