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), # 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_men_acquirefrom_women = 50 # number of contacts per year
contact_rate_women_acquirefrom_men = 50 # number of contacts per year
transmission_probability_men_from_women = 0.05 # transmission probability
transmission_probability_women_from_men = 0.05 # transmission probability
infectious_period = 1            # infectious period (year)

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

beta_mw_value = contact_rate_men_acquirefrom_women * transmission_probability_men_from_women

beta_wm_value = contact_rate_women_acquirefrom_men * transmission_probability_women_from_men

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, 10, 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, Infectious and Recovered 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') 
text(30,0.2, "Sm", col = 'blue')

# 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) 
text(5,0.5, "Im", col = 'red')

# remain on same frame
par (new = TRUE)  

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

# remain on same frame
par (new = TRUE)  

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

# remain on same frame
par (new = TRUE)