In collinmmccabe/enss: Effective Network Size Simulation

library(statnet)
library(igraph)

Null Model: Maximally-complete, unweighted networks

First, we generated sociomatrices with '1' signifying ties between a row index and a column index. Because we were focused on social networks, we replaced the diagonals of the matrices with '0' to show that there are no self-referencing, or looping ties.

  net <- matrix(data=1, nrow=6, ncol=6)
  diag(net) <- 0

net

We then generated a statnet network-type object to facilitate network-specific analyses.

  net <- network(net)

And we represented these as edgelists to make referencing random edges for disease spread more intuitive. We calculated the number of edges and vertices for each network as well, to reduce computation time during simulations.

  netud <- graph.edgelist(as.edgelist(net)[,], directed=FALSE)
  networkud <- cbind(unique(get.edgelist(netud)), weight=1,
                     n_nodes=c(vcount(netud), rep(NA, ecount(netud)/2-1)),
                     n_edges=c(ecount(netud)/2, rep(NA, ecount(netud)/2-1)))

networkud
max_n=50

We generated these graphs for network sizes from 6 (the minimum number of nodes required for modular networks) to 50.

Transmission Modes

Once we had our sample of null-model (maximally complete) networks, we simulated the spread of generic diseases. All diseases shared the same parameters for beta (the infection rate, or how likely a susceptible individual was to become infected upon interaction with an infected individual), gamma where applicable (the recovery rate, or the daily likelihood than an infected individual would recover to either a resistant state- for SIR, or a susceptible state- for SIS), per captia interaction rate per day (scaling with network size), and the number of days for which a disease would be simulated. For STD models, we also included infection rate modifiers, because transmission rate of an STD is dependent on the sex of the infected individual and the sex of the susceptible:

Infected -> | Male | Female -------------:|------|-------- Male | 3 | 0.5
Female | 1 | 0.1

We used iterative (looping), edge-selection-based algorithms for simulating the spread of disease through our populations. The algorithms for each of the disease modes are listed below.

SIS

SIS models assume that after a period of time, infected individuals will return to a susceptible state. So, as long as at least one individual in the group is infected at anyone time, the proportion of infected individuals to susceptible individuals will tend to oscillate around an equilibrium level. As this is the case, we are more concerned with the equilibrium ratio of infecteds to susceptibles in the population, and this is returned from the following function.

sim_SIS <- function(networkud, beta, gamma, intxn_per_day, days) {

  n = networkud[1,4]
  e = networkud[1,5]
  cdata <- networkud[,1:2]

  infection_status <- c(rep(1,n))
  index_infected <- sample(1:n, 1)
  infection_status[index_infected] = 2

  day_counter <- 0
  while(day_counter <= days) {

    int_counter <- 0
    while(int_counter <= intxn_per_day*n) {

      selected_edge <- sample(1:e,1)

      if (sum(infection_status[cdata[selected_edge,1:2]]) == 3) {
        if (beta >= runif(1,0,1)) { 
          infection_status[cdata[selected_edge,1:2]] = 2
        }
      }

      int_counter <- sum(int_counter,1)
    }

    for (j in which(infection_status %in% 2)) { 
      if (gamma >= runif(1,0,1)) {
        infection_status[j] = 1
      }
    }
    day_counter <- sum(day_counter,1)
    if (sum(infection_status%%2) == n) break
  }
  return(c(day_counter-1,sum(infection_status == 1),sum(infection_status == 2)))
}

SI

SI diseases assume no recovery, and so these will in theory spread to every susceptible individual in a population eventually, as long as individuals do not die and all individuals are connected in the network. In order to make sure that our diseases are reaching every individual, we reported the final infected ratio, assuming that all will reach 100%.

sim_SI <- function(networkud, beta, intxn_per_day, days) {

  n = networkud[1,4]
  e = networkud[1,5]
  cdata <- networkud[,1:2]

  infection_status <- c(rep(1,n))
  index_infected <- sample(1:n, 1)
  infection_status[index_infected] = 2

  day_counter <- 0
  while(day_counter <= days) {

    int_counter <- 0
    while(int_counter <= intxn_per_day*n) {

      selected_edge <- sample(1:e,1)

      if (sum(infection_status[cdata[selected_edge,1:2]]) == 3) {
        if (beta >= runif(1,0,1)) { 
          infection_status[cdata[selected_edge,1:2]] = 2
        }
      }

      int_counter <- sum(int_counter,1)
    }

    day_counter <- sum(day_counter,1)
    if (sum(infection_status%%2) == 0) break
  }
  return(c(day_counter-1,sum(infection_status == 1),sum(infection_status == 2)))
}

STD

STDs are really just special cases of SI diseases, where the transmission rate varies based on the sex of interacting individuals. The beta modifiers for each type of interaction (male-to-female, MM; female-to-male, FM; male-to-male, MM; and female-to-female, FF) are therefore listed as additional inputs for the STD simulation function.

sim_STD <- function(networkd, beta, intxn_per_day, days, MM, MF, FM, FF) {

  n = networkd[1,4]
  e = networkd[1,5]
  cdata <- networkd[,1:2]
  sexes <- networkd[1:n,6]

  infection_status <- c(rep(1,n))
  index_infected <- sample(1:n, 1)
  infection_status[index_infected] = 2

  day_counter <- 0
  while(day_counter <= days) {

    int_counter <- 0
    while(int_counter <= intxn_per_day*n) {

      selected_edge <- sample(1:e,1)

      if (sum(infection_status[cdata[selected_edge,1:2]]) == 3) {
        sex_ind = 0; beta_mod = 0
        if (infection_status[cdata[selected_edge,1]] == 2) {
          sex_ind <- (sexes[cdata[selected_edge,1]] * 2) - sexes[cdata[selected_edge,2]]+1
        } else {
          sex_ind <- (sexes[cdata[selected_edge,2]] * 2) - sexes[cdata[selected_edge,1]]+1
        }

        switch(sex_ind, {beta_mod <- beta * MF}, {beta_mod <- beta * MM},
               {beta_mod <- beta * FF}, {beta_mod <- beta * FM})

        if (beta_mod >= runif(1,0,1)) { 
          infection_status[cdata[selected_edge,1:2]] = 2
        }
      }

      int_counter <- sum(int_counter,1)
    }

    day_counter <- sum(day_counter,1)
    if (sum(infection_status%%2) == 0) break
  }
  return(c(day_counter-1,sum(infection_status == 1),sum(infection_status == 2)))
}

SIR

Finally, the most complex model we tested, SIR, assumes that after a period of infection, individuals recover and become immune to future infection. For this reason, we recorded the maximum number of individuals infected at any point in the simulation to gauge the peak prevalence of the disease.

sim_SIR <- function(networkud, beta, gamma, intxn_per_day, days) {

  n = networkud[1,4]
  e = networkud[1,5]
  cdata <- networkud[,1:2]

  infection_status <- c(rep(1,n))
  index_infected <- sample(1:n, 1)
  infection_status[index_infected] = 2

  max_infected <- 1

  day_counter <- 0
  while(day_counter <= days) {

    int_counter <- 0
    while(int_counter <= intxn_per_day*n) {

      selected_edge <- sample(1:e,1)

      if (sum(infection_status[cdata[selected_edge,1:2]]) == 3) {
        if (beta >= runif(1,0,1)) { 
          infection_status[cdata[selected_edge,1:2]] = 2
        }
      }

      int_counter <- sum(int_counter,1)
    }

    for (j in which(infection_status %in% 2)) { 
      if (gamma >= runif(1,0,1)) {
        infection_status[j] = 3
      }
    }

    curr_infected <- sum(infection_status == 2)
    if (curr_infected > max_infected) {
      max_infected <- curr_infected
    }

    day_counter <- sum(day_counter,1)

    if (sum(infection_status%%2) == n) break
  }
  return(c(day_counter-1, sum(infection_status == 1), sum(infection_status == 2),
           sum(infection_status == 3), max_infected))
}