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, 10, 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.4     # 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, 10, by=1/100)   # timeline of 10 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 (0, 10), ylim = range (0, 0.05) )        # susceptible hosts over time
plot ( I ~ time, data = output, type='l', col = 'red', xlim = range (0, 10), 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 (0, 10), 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 (0, 10), 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 (0, 10), 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 (0, 10), ylim = range (0, 1))

This models aims to show and predict when an outbreak may occur and how intense the burden of disease of the outbreak will be. By year 10, the number of infectious and recovered hosts decline while the number of scusceptible hosts increase. There is no sign of an outbreak by year 10, but due to the number of scuscpetible hosts still remaining, it is possible another outbreak can occur.