Objective

To analyze the dynamics of the Susceptibles-Infectious-Susceptible (SIS) epidemic model as applied to a novel heterosexually transmitted infection.

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

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]        # susceptible 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, dSw, dIm, dIw)
           list (results)
         }
    )
}

Parameters

contact_rate_mw = 0.137                   # number of contacts per year in men
contact_rate_wm = 0.137                   # number of contacts per year in women
transmission_probability_mw = 0.05      # transmission probability in men
transmission_probability_wm = 0.1      # transmission probability in women
infectious_period = 365               # infectious period

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

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

Compute Ro - Reproductive number.

Ro_m = beta_value_mw / gamma_value
Ro_w = 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.

Xm = 10000        # susceptible men
Xw = 11000        # susceptible women
Ym = 1        # infectious men
Yw = 1        # infectious women

Compute total population.

N = Xm + Xw + Ym + Yw

Initial state values for the differential equations.

initial_values = c (Sm = Xm/N, Sw = Xw/N, Im = Ym/N, Iw = Yw/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 Susceptibles male sub-population.

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

Plot dynamics of Infectious male sub-population.

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

Plot dynamics of Susceptibles female sub-population.

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

Plot dynamics of Infectious female sub-population.

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

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

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

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

# remain on same frame
par (new = TRUE)    

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

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

Description and Results

The above graphs that susceptible men have been increasing over time while infectious men have been decreasing. The same applies to susceptible and infectious women, although their curves aren’t quite as sharp as the men’s over a 10 years period of time. The combined plots show that women are slightly more susceptible than men, and the same also appears to apply in the infectious category as well.

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 reproductive 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 proportion’s 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 male 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 male population using contraception \n", ylab = "Effective Reproductive Number for men (Re_m)")

Plot proportion of female 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)")

Description and Results

The below graphs show that as the proportion of men and women using contraception increases, the reproductive number of heterosexually transmitted infections deceased. If contraception use is 80%, then the minimum population proportion (or herd immunity) that needs to be maintained in order to prevent an outbreak among men is 96% and 98% for women.