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.

sis_model = function (current_timepoint, state_values, parameters)
{
  # create state variables (local variables)
  Sm = state_values [1]        # susceptible men
  Im = state_values [2]        # infectious men
  Sw = state_values [3]        # susceptible women
  Iw = state_values [4]        # infectious women
  TPwm = state_values [5]        # estimate proportion women
  TPmw = state_values [6]        # estimate proportion men
  
  with ( 
    as.list (parameters),     # variable names within parameters can be used 
         {
           # compute derivatives
           dSm = (-beta_mw * Sm * Iw) + (gamma * Im)
           dIm = ( beta_mw * Sm * Iw) - (gamma * Im)
           dSw = (-beta_wm * Sw * Im) + (gamma * Iw)
           dIw = ( beta_wm * Sw * Im) - (gamma * Iw)
           dTPwm = (contact_rate_wm * transmission_probability_wm * disease_protection * compliant_population)
           dTPmw = (contact_rate_mw * transmission_probability_mw * disease_protection * compliant_population)
           
           # combine results
           results = c (dSm, dIm, dSw, dIw, dTPwm, dTPmw)
           list (results)
         }
    )
}

Parameters

contact_rate_mm = 0                     # number of contacts per day
contact_rate_mw = (50/365)              # number of contacts per day
contact_rate_ww = 0                     # number of contacts per day
contact_rate_wm = (50/365)              # number of contacts per day

transmission_probability_mm = 0       # transmission probability
transmission_probability_mw = 0.05       # transmission probability
transmission_probability_ww = 0       # transmission probability
transmission_probability_wm = 0.1       # transmission probability
disease_protection = .20              # protected population
compliance = .95                     # persons compliant with intervention
infectious_period = 365                # infectious period

Compute values of compliant_population (elimination proportion), beta_mw (tranmission rate), beta_wm (tranmission rate) and gamma (recovery rate)

beta_value_mw = contact_rate_mw * transmission_probability_mw
beta_value_wm = contact_rate_wm * transmission_probability_wm
gamma_value = 1 / infectious_period
compliant_population = 1 - compliance

Compute Ro - Reproductive number.

Rom = beta_value_mw / gamma_value
Row = beta_value_wm / gamma_value

Disease dynamics parameters.

parameter_list = c (compliant_population, beta_mw = beta_value_mw, beta_wm = beta_value_wm, gamma = gamma_value)

Initial values for sub-populations.

W = 10000     # susceptible hosts
X = 1           # infectious hosts
Y = 11000     # susceptible hosts
Z = 1           # infectious hosts
TPwm = 1
TPmw = 1

Compute total population.

N = W + X + Y + Z + TPmw + TPwm

Initial state values for the differential equations.

initial_values = c (Sm = W/N, Im = X/N, Sw = Y/N, Iw = Z/N, TPwm = 1, TPmw = 1)

Output timepoints.

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

Simulate the SIS epidemic.

output = lsoda (initial_values, timepoints, sis_model, parameter_list)

Plot dynamics of Susceptible Men sub-population.

plot (Sm ~ time, data = output, type='b', col = 'blue')

Plot dynamics of Infectious Men sub-population.

plot (Im ~ time, data = output, type='b', col = 'red')

Plot dynamics of Susceptible Women sub-population.

plot (Sw ~ time, data = output, type='b', col = 'orange')

Plot dynamics of Infectious Women sub-population.

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

Plot dynamics of Elimination Proportion Women.

plot (TPwm ~ time, data = output, type='b', col = 'purple')

Plot dynamics of Elimination Proportion Men.

plot (TPmw ~ time, data = output, type='b', col = 'black')

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

# susceptible men hosts over time
plot (Sm ~ time, data = output, type='b', ylim = c(0,1), col = 'blue', ylab = 'S, I, S', main = 'SIS epidemic') 

# remain on same frame
par (new = TRUE)    

# infectious men hosts over time
plot (Im ~ time, data = output, type='b', ylim = c(0,1), col = 'red', ylab = '', axes = FALSE) 

# remain on same frame
par (new = TRUE)  

# elimiation men hosts over time
plot (Im ~ time, data = output, type='b', ylim = c(0,1), col = 'black', ylab = '', axes = FALSE) 

# remain on same frame
par (new = TRUE)  

# elimiation men hosts over time
plot (Im ~ time, data = output, type='b', ylim = c(0,1), col = 'purple', ylab = '', axes = FALSE) 

# remain on same frame
par (new = TRUE)  

# susceptible women hosts over time
plot (Sw ~ time, data = output, type='b', ylim = c(0,1), col = 'orange', ylab = 'S, I, S', main = 'SIS epidemic') 

# remain on same frame
par (new = TRUE)    

# infectious women hosts over time
plot (Iw ~ time, data = output, type='b', ylim = c(0,1), col = 'green', ylab = '', axes = FALSE)