This program simulates the temporal seasonal dynamics of measles outbreaks.

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.

# function 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 measles 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')        # recovered hosts over time
plot (R ~ time, data = output, type='l', col = 'green', xlim = range (10, 20), ylim = range (0.95, 1))     

# reset to prior graphics settings
par (prior_settings)

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

# susceptible hosts over time
plot (S ~ time, data = output, type='l', col = 'blue', xlim = range (10, 20), ylim = range (0, 1), ylab = 'S, I, R', main = 'SIR epidemic' ) 

# remain on same frame
par (new = TRUE)    

# infectious hosts over time
plot ( I ~ time, data = output, type='l', col = 'red', xlim = range (10, 20), ylim = range (0, 1) ) 

# remain on same frame
par (new = TRUE)  

# recovered hosts over time
plot (R ~ time, data = output, type='l', col = 'green', xlim = range (10, 20), ylim = range (0, 1))