remove (list = objects())
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
  
  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)
    
    # combine results
    results = c (dSm, dIm, dSw, dIw)
    list (results)
    }
  )
}

Parameters

contact_rate = 50                      # number of contact per year
transmission_probability_mw = 0.05     # transmission probability men from women
transmission_probability_wm = 0.1      # transmission probability women from men
infectious_period = 1                  # infectious period

Compute values of beta (transmission rate) and gamma (recovery rate).

beta_value_mw = contact_rate * transmission_probability_mw
beta_value_wm = contact_rate * transmission_probability_wm
gamma_value = 1 / infectious_period

Estimate Ro using contact rate, transmission probability and infectious period.

contact_rate = 50                      
transmission_probability_mw = 0.05     
transmission_probability_wm = 0.1      
infectious_period = 1
R0m = contact_rate * transmission_probability_mw * infectious_period
R0w = contact_rate * transmission_probability_wm * infectious_period
c ("R0m = ", R0m)
## [1] "R0m = " "2.5"
c ("R0w = ", R0w)
## [1] "R0w = " "5"

Disease dynamics parameters

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

Initial values for sub-populations

W = 10000      # susceptible men
X = 1          # infectious men
Y = 11000      # susceptible women
Z = 1          # infectious women

Compute total population

N = W + X + Y + Z

Initial state values for the differential equations.

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

Output timepoints.

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

Simulate the SIS epidemic

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

Plot dynamics of Susceptible and Infectious sub-populations.

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

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

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

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

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

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

# remain on same frame
par (new = TRUE) 

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

# remain on same frame
par (new = TRUE) 
## Warning in par(new = TRUE): calling par(new=TRUE) with no plot

Compute new transmission probability considering contraception use provides 80% protection

new_transmission_probability_mw = contact_rate * (transmission_probability_mw * 0.2)
new_transmission_probability_wm = contact_rate * (transmission_probability_wm * 0.2)
c ("new_transmission_probability_mw = ", new_transmission_probability_mw)
## [1] "new_transmission_probability_mw = "
## [2] "0.5"
c ("new_transmission_probability_wm = ", new_transmission_probability_wm)
## [1] "new_transmission_probability_wm = "
## [2] "1"

Estimate Ro using contact rate, new transmission probability (using contraception) and infectious period

R0_m = contact_rate * new_transmission_probability_mw * infectious_period
R0_w = contact_rate * new_transmission_probability_wm * infectious_period
c ("R0_m = ", R0_m)
## [1] "R0_m = " "25"
c ("R0_w = ", R0_w)
## [1] "R0_w = " "50"

Generate sequence of numbers of fraction of population using contraception from 0 to 1

fraction_contraception_m = seq (0, 1, 0.2)
fraction_contraception_w = seq (0, 1, 0.15)

Print fraction of populations using contraception

cat ("Fraction of men population using contraception:", fraction_contraception_m)
## Fraction of men population using contraception: 0 0.2 0.4 0.6 0.8 1
cat ("Fraction of women population using contraception:", fraction_contraception_w)
## Fraction of women population using contraception: 0 0.15 0.3 0.45 0.6 0.75 0.9

Compute effective reproductive numbers

Re_m = R0_m * (1 - fraction_contraception_m)
Re_w = R0_w * (1 - fraction_contraception_w)

Print effective reproductuve numbers

cat ("Effective reproductive number for men:", Re_m)
## Effective reproductive number for men: 25 20 15 10 5 0
cat ("Effective reproductive number for women:", Re_w)
## Effective reproductive number for women: 50 42.5 35 27.5 20 12.5 5

Compute minimum proportion of populations that should consistently use contraception to eliminate the outbreak (Herd immunity threshold)

population_proportion_threshold_m = 1 - (1/R0_m)
population_proportion_threshold_w = 1 - (1/R0_w)

Print population proportions threshold

cat ("population proportion threshold men= ", population_proportion_threshold_m)
## population proportion threshold men=  0.96
cat ("population proportion threshold women= ", population_proportion_threshold_w)
## population proportion threshold women=  0.98

Plot proportion of men population using contraception (versus) effective reproductive number

subtitle = paste ("R0_m = ", R0_m, ", population proportion threshold men = ", round (population_proportion_threshold_m, digits = 4), "; men", sep="")
plot (fraction_contraception_m, Re_m, main = "Contraception Intervention" , sub = subtitle, xlab = "Fraction of men population using contraception \n", ylab = "Effective Reproductive Number for men (Re_m)")

Plot proportion of women population using contraception (versus) effective reproductive number

subtitle = paste ("R0_w = ", R0_w, ", population proportion threshold women = ", round (population_proportion_threshold_w, digits = 4), "; women", sep="")
plot (fraction_contraception_w, Re_w, main = "Contraceptioin Intervention" , sub = subtitle, xlab = "Fraction of women population using contraception \n", ylab = "Effective Reproductive Number for women (Re_w)")

This model is used to analyze, simulate, and estimate the minimum proportion of both men and women populations that are needed to elminate an STI outbreak given that the contraception intervention gives 80% protection against the disease.The resulting graphs reveal that as the proportion of both the mean and women populations using contraception increases, the reproductive number forthe disease decreases. At minimum, 96% of men and 98% of women need to consistently used contraception in order to eliminate the outbreak.