Remove all objects from workspace.

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

creat equations

sis_model = function (current_timepoint, state_values, parameters)
{  
  # create state variables (local variables)
  Sm = state_values [1]        # susceptibles men
  Im = state_values [2]        # infectious men
  Sw = state_values [3]        # susceptibles women
  Iw = state_values [4]        # infectious women
  
  with (as.list (parameters), 
         { # 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)
           
           # combine results
           results = c (dSm, dIm, dSw, dIw)
           list (results)
         } 
        )
}

Parameters

Assuming 80% protection by contraception and 0.15% both men and women are compliance

contact_rate_mw = 50/365 
contact_rate_wm = 50/365 
transmission_probability_mw = 0.05 * 0.2 *(1-0.15) # transmission prob*efficiency of protection* complian
transmission_probability_wm = 0.1 * 0.2 * (1-0.15)
infectious_period = 365

Compute values of beta (tranmission rate = contact rate*transmission prob) and gamma.

beta_mw_value = (50/365) * 0.05 * 0.2 * (1-0.15)

beta_wm_value = (50/365) * 0.1 * 0.2 * (1-0.15)

gamma_value = 1 /infectious_period

Compute Ro - Reproductive number.

Ro.m = beta_mw_value / gamma_value
Ro.W = beta_wm_value / gamma_value

Disease dynamics parameters.

parameter_list = c (beta_mw = beta_mw_value, beta_wm = beta_wm_value ,gamma = gamma_value)

Initial values for sub-populations.

A = 10000 # susceptible men
B = 11000 # susceptible women
C = 1  # infectious men
D = 1 # infectious women

Compute total population.

N = A+B+C+D

Initial state values for the differential equations.

initial_values = c (Sm = A/N, Im = C/N, Sw = B/N, Iw = D/N )

Output timepoints.

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

Simulate the SIS epidemic.

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

Plot dynamics of Susceptibles sub-population.

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

Plot dynamics of Infectious sub-population.

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

Plot dynamics of Susceptibles sub-population.

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

Plot dynamics of Infectious sub-population.

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

Plot dynamics of Susceptibles and infectious sub-populations in the same plot.

# susceptible men over time
plot (Sm ~ time, data = output, type='b', ylim = c(-0.2,0.8), col = 'blue', ylab = 'S, I, S', main = 'SIS Heterogeneous') 
text(1500,0.4, "Sm", col = 'blue')

# remain on same frame
par (new = TRUE)    


# susceptible women over time
plot (Sw ~ time, data = output, type='b', ylim = c(-0.2,0.8), col = 'purple', ylab = '', axes = FALSE) 
text(1500,0.6, "Sw", col = 'purple')

# remain on same frame
par (new = TRUE) 

# infectious men over time
plot (Im ~ time, data = output, type='b', ylim = c(-0.2,0.8), pch = 18,  col = 'red', ylab = '', axes = FALSE) 
text(1500,0.05, "Im", col = 'red')

# remain on same frame
par (new = TRUE) 

# infectious women over time
plot (Iw ~ time, data = output, type='b', ylim = c(-0.2,0.8), col = 'green', ylab = '', axes = FALSE) 
text(1500,-0.05, "Iw", col = 'green')

Description: If transmission probability of men acquired from women is 0.05, transmission probability of women acquired from men is 0.1 and contraception provide 80% protection. The minimum compliance population is only 0.15%.