TSIR EPIDEMIOLOGICAL MODEL

Objective: To analyse the seasonal dynamics of measles using the Susceptible-Infectious-Recovered (SIR) epidemic model.

Clear workspace and load required packages.

# remove all objects from workspace
remove (list = objects())

#Load add-on packages 
library (deSolve)     #differential equation solver
## 
## Attaching package: 'deSolve'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot
# contains lsoda function - differential equation solver

Simulate sinusoidal dynamics of cosine function.

# timepoints
timepoints = seq (0, 2, by = 1/100)

# compute the cosine function
sinusoidal_dynamics_sin = cos (2 * pi * timepoints)

Plot sinusoidal dynamics of cosine function.

# plot cosine function
plot (timepoints, sinusoidal_dynamics_sin, type = "l", main = "Sinusoidal dynamics of cosine function", xlab = "Time", ylab = "cosine value")

Define the SIR model with seasonal dynamics

# finction to compute derivatives of the differential equations
seasonal_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 time-dependent beta
    beta = beta0 * (1 + beta1 * cos (2 * pi * current_timepoint))
    
    # compute derivatives
    dS = mu- (beta * S * I) - (mu * S)
    dI = (beta * S * I) - (gamma * I) - (mu * I)
    dR = (gamma * I) - (mu * R)
    
    # combine results
    results = c (dS, dI, dR)
    list (results)
    }
  )
}

Initialize parameters of measleas transmission dynamics.

# parameters
mu_value = 1/50                           # death rate -- unit: 1/year
infectious_period = 13                    # infectious period (days)
gamma_value = 365 / infectious_period     # gamma (recovery rate) -- unit: 1/year

# transmission rates
beta0_value = 1000       # basic transmission rate
beta1_value = 0.2        # amplitude of seasonality

# disease dynamics parameters
parameter_list = c (mu = mu_value, beta0 = beta0_value, beta1 = beta1_value, gamma = gamma_value)

# initial values for sub-populations
X = 999        # susceptible hosts
Y = 1          # infectious hosts
Z = 0          # 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)

Set simulation timeline.

# output timepoints
timepoints = seq (0, 20, by=1/100)     # timeline of 20 years

Simulate the seasonal measles epidemic.

# simulate the seasonal SIR epidemic
output = lsoda (initial_values, timepoints, seasonal_sir_model, parameter_list)

Plot the results.

# set graphics settings
prior_settings = par (mfrow = c (2, 2))           # (1 rows * 2 columns) graphs on 1 plot

# plot results of disease dynamics
plot (S ~ time, data = output, type='l', col = 'blue', xlim = range (10, 20), ylim = range (0, 0.05))       # susceptible hosts over time
plot (I ~ time, data = output, type='l', col = 'red', xlim = range (10, 20), ylim = range (0, 0.01))       # infectious hosts over time
plot (R ~ time, data = output, type='l', col = 'green', xlim = range (10, 20), ylim = range (0.95, 1))       # Susceptible hosts over time

# reset to prior graphics settings 
par (prior_settings)