Remove all objects from workspace.

r 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.

sirc_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
  C = state_values [4]        # carriers
  
  with ( 
    as.list (parameters),     # variable names within parameters can be used 
         {
           # compute derivatives
           dS = (-beta * S * I) - (epsilon * beta * S * C)
           dI = (beta * S * I) + (epsilon * beta * S * C) - (gamma * I)
           dR = (1 - rho) * (gamma * I) + (tau * C) 
           dC = (rho * gamma * I) - (tau * C)
           
           # combine results
           results = c (dS, dI, dR, dC)
           list (results)
         }
    )
}

Parameters

contact_rate = 7                 # number of contacts per day
transmission_probability = 0.08      # transmission probability
infectious_period = 11             # infectious period
reduced_transmission_rate = 0.09     # chronic carriers compared to acute infections
acute_infections_proportion = 0.67     # acute infections that become carriers
length_of_time_in_carrier_state = 7     # length of time in carrier state

Compute values of beta (tranmission rate), epsilon (reduced transmission rate), tau (carrier state), rho (acute infections) and gamma (recovery rate)

beta_value = contact_rate * transmission_probability
gamma_value = 1 / infectious_period
epsilon_value = 0.09
tau_value = 1 / length_of_time_in_carrier_state
rho_value = 0.67

Compute Ro - Reproductive number.

Ro = beta_value / gamma_value

Disease dynamics parameters.

parameter_list = c (beta = beta_value, gamma = gamma_value, epsilon = epsilon_value, tau = tau_value, rho = rho_value)

Initial values for sub-populations.

V = 2          # carrier hosts
X = 17282      # susceptible hosts
Y = 19          # infectious hosts
Z = 768         # recovered hosts

Compute total population.

N = V + X + Y + Z

Initial state values for the differential equations.

initial_values = c (S = X/N, I = Y/N, C = V/N, R = Z/N)

Output timepoints.

timepoints = seq (0, 50, by=1)

Simulate the SIR/C epidemic.

output = lsoda (initial_values, timepoints, sirc_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 Carriers sub-population.

plot (C ~ time, data = output, type='b', col = 'pink')

Plot dynamics of Recovered sub-population.

plot (R ~ time, data = output, type='b', col = 'green')

Plot dynamics of Susceptibles, Infectious, Carriers 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, C', main = 'SIR/C epidemic') 

# 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) 

# remain on same frame
par (new = TRUE)  

# carrier hosts over time
plot (C ~ time, data = output, type='b', ylim = c(0,1), col = 'pink', ylab = '', axes = FALSE) 

# 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)

In the Susceptible Infectious Recovered/Carrier epidemiological model, a susceptible host becomes infectious, becomes a carrier, and/or recovers. They also may just become a carrier and then recover. In the above figure, the number of susceptible (blue), infectious (red), carrier (pink) and recovered hosts (green) are illustrated.