SIS model: Heterosexually transmitted infection

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

Compute Ro - Reproductive number.

Ro_mw = beta_value_mw / gamma_value
Ro_wm = beta_value_wm / gamma_value

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 sub-populations.

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

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

COMMENT: The graphs show that the number of susceptible men and women is decreasing over the 10-year period.

Plot dynamics of infectious sub-populations.

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

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

COMMENT: The graphs show that the number of infectious men and women is increasing over the 10-year period

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,1), col = 'blue', ylab = 'S, I, S', main = 'SIS epidemic') 

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

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

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 with interval 0.1 interval

fraction_contraception_m = seq (0, 1, 0.1)
fraction_contraception_w = seq (0, 1, 0.1)

Print fraction of populations using contraception

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

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 22.5 20 17.5 15 12.5 10 7.5 5 2.5 0
cat ("Effective reproductive number for women:", Re_w)
## Effective reproductive number for women: 50 45 40 35 30 25 20 15 10 5 0

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

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 = "sti" , 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 = "sti" , sub = subtitle, xlab = "Fraction of women population using contraception \n", ylab = "Effective Reproductive Number for women (Re_w)")

COMMENT: We can see from the graphs that as the proportion of population (men and women) using contraception increases, the reproductive number for STI decreases. The minimum proportions of population that should consistently use contraception to eliminate the outbreak (cfr. Herd immunity) are respectively 96% and 98% for men and women, considering that contraception provides 80% protection against the disease.