####################################
# Parameters
####################################
force.feml.monog <- F # Should we enforce a rule of female monogamy?
force.male.monog <- F # Should we enforce a rule of male monogamy?
n.femls <- 5000 # Number of females
n.males <- 5000 # Number of females
expected.meandeg <- 0.8 # Mean degree (= average # of relationships a person is at any moment in time
avg.duration <- 10 # Mean duration of a relationship (in months)
numtimesteps <- 1200 # Length of simulations (in months)
init.number.of.infected.femls <- 100 # How many females are infect at the outset of the simulation?
init.number.of.infected.males <- 100 # How many males are infect at the outset of the simulation?
beta.by.time.since.inf <- c(rep(0.2055,3), # Probability of transmission per month for active relationship
rep(0.0088,100),rep(0.0614,9),rep(0,10)) # From Hollingsworth et al. 2008
####################################
# Basic calculations and book-keeping
####################################
n.pop <- n.femls + n.males # Total pop size
expected.edges <- round(expected.meandeg*(n.pop)/2) # Expected # of relationships ("edges") in the population at any time
prob.dissolution <- 1/avg.duration # Daily prob of dissolution for an existing relationship
time.of.aids.death <- length(beta.by.time.since.inf) # Length of time from HIV infection until death (in months)
cum.num.feml.aids.deaths <- 0 # Cumulative # of female AIDS deaths
cum.num.male.aids.deaths <- 0 # Cumulative # of female AIDS deaths
feml.prev <- vector() # Sets up a vector to store female prevalence at each time step
male.prev <- vector() # Sets up a vector to store male prevalence at each time step
####################################
# Creating the female data frame
####################################
femls <- data.frame(row.names=1:n.femls) # Creates a data frame to store info about the females. One row per female.
femls$hiv.status <- rep(0,n.femls) # Create a variable in the female data frame called "hiv.status". Give all females the value 0 (for the moment).
femls$inf.time <- NA # Create a variable in the female data frame called "inf.time" (infection time). Give all females an NA (for the moment).
if (init.number.of.infected.femls>0) { # If there are any initially infected females,
init.inf.f <- sample(1:n.femls, init.number.of.infected.femls) # Sample from vector of female ids with size of initially infected females
femls$hiv.status[init.inf.f] <- 1 # Give them infection
init.inf.time.f <- sample(0:(-time.of.aids.death+2),
init.number.of.infected.femls, replace=TRUE) # Sample times backwards from present to length of infection
femls$inf.time[init.inf.f] <- init.inf.time.f # Assign.
}
####################################
# Male data frame
####################################
# All parallel to female data frame above
males <- data.frame(row.names=1:n.males)
males$hiv.status <- 0
males$inf.time <- NA
if (init.number.of.infected.males>0) {
init.inf.m <- sample(1:n.males, init.number.of.infected.males) # Sample from vector of female ids with size of initially infected females
males$hiv.status[init.inf.m] <- 1 # Give them infection
init.inf.time.m <- sample(0:(-time.of.aids.death+2),
init.number.of.infected.males, replace=TRUE) # Sample times backwards from present to length of infection
males$inf.time[init.inf.m] <- init.inf.time.m # Assign.
}
####################################
# Initial contact network
####################################
edgelist <- data.frame(row.names=1:expected.edges) # Create a data frame that will be used to hold the IDs of the relationship pairs
if (force.feml.monog==FALSE) { # If we are *not* enforcing female monogamy,
edgelist$f <- sample(1:n.femls,expected.edges,replace=T) # sample from the female IDs with replacement, and assign them to a column in the data frame called "f"
} else { # otherwise,
edgelist$f <- sample(1:n.femls,expected.edges,replace=F) # sample from the female IDs without replacement, and assign them to a column in the data frame called "f"
}
if (force.male.monog==FALSE) { # Same for males
edgelist$m <- sample(1:n.males,expected.edges,replace=T)
} else {
edgelist$m <- sample(1:n.males,expected.edges,replace=F)
}
table(tabulate(edgelist$f)) # These rather odd lines provide a check to see a table of relationships per person in the population. When monogamy is enforced for
table(tabulate(edgelist$m)) # either sex, the corresponding table should be limited to 0's and 1's; when it os not enforced, the table is not limited.
#################################### # Everything until now was set-up; now this is the core of the simulation.
# Time loop # We will simulate forward in time, considering transmission, vital dynamics, and relational dynamics at each time point.
####################################
for (time in 1:numtimesteps) {
# Transmissions
hiv.status.feml.partners <- femls$hiv.status[edgelist$f] # Create a vector containing the HIV status for the female partner in each relationship
hiv.status.male.partners <- males$hiv.status[edgelist$m] # Create a vector containing the HIV status for the male partner in each relationship
sdpf <- which(hiv.status.male.partners==0 & # Make a vector of IDs for the relationships that are serodiscordant with positive female (SDPF)
hiv.status.feml.partners == 1)
f.in.sdpf <- edgelist$f[sdpf] # Make a vector of IDs for the females in the SDPF relationships
inf.time.sdpf <- time - femls$inf.time[f.in.sdpf] # Make a vector that identifies how long before present these females were infected
prob.trans.sdpf <- beta.by.time.since.inf[inf.time.sdpf] # Determine probability of transmission in the SDPM relationships, based on time since the male was infected
trans.sdpf <- rbinom(length(sdpf),1,prob.trans.sdpf) # Flip a weighted coin for each SDPM relationship to determine if transmission occurs
newly.inf.males <- edgelist$m[sdpf[trans.sdpf==1]] # Double-indexing (an R specialty!) This lets you get the identities of the newly infected males in one step.
males$hiv.status[newly.inf.males] <- 1 # Assign the newly infected males an hiv.status of 1
males$inf.time[newly.inf.males] <- time # Assign the newly infected males an infection time of time
sdpm <- which(hiv.status.feml.partners==0 & # All parallel for serodiscordant with positive male (SDPM)
hiv.status.male.partners == 1)
m.in.sdpm <- edgelist$m[sdpm]
inf.time.sdpm <- time - males$inf.time[m.in.sdpm]
prob.trans.sdpm <- beta.by.time.since.inf[inf.time.sdpm]
trans.sdpm <- rbinom(length(sdpm),1,prob.trans.sdpm)
newly.inf.femls <- edgelist$f[sdpm[trans.sdpm==1]]
femls$hiv.status[newly.inf.femls] <- 1
femls$inf.time[newly.inf.femls] <- time
# Deaths to AIDS
femls.dying.of.AIDS <- which(time-femls$inf.time==time.of.aids.death) # Which females have been HIV+ long enough to die of AIDS?
males.dying.of.AIDS <- which(time-males$inf.time==time.of.aids.death) # Which females have been HIV+ long enough to die of AIDS?
cum.num.feml.aids.deaths <- # Increase cumulative # of female AIDS deaths by the newly dying females
cum.num.feml.aids.deaths + length(femls.dying.of.AIDS)
cum.num.male.aids.deaths <- # Increase cumulative # of male AIDS deaths by the newly dying males
cum.num.male.aids.deaths + length(males.dying.of.AIDS)
# End of ties because of death
edges.with.feml.dying.of.aids <- # Determine the IDs of those relationships involving dying women
which(edgelist$f %in% femls.dying.of.AIDS)
edges.with.male.dying.of.aids <-
which(edgelist$m %in% males.dying.of.AIDS) # Determine the IDs of those relationships involving dying women
edges.with.either.dying.of.aids <-
c(edges.with.feml.dying.of.aids,edges.with.male.dying.of.aids) # Combine the two
if (length(edges.with.either.dying.of.aids>0)){ # Remove those edges from the edgelist
edgelist <- edgelist[-edges.with.either.dying.of.aids,]
}
# End of other ties randomly
edges.coinflip <- rbinom(dim(edgelist)[1],1,prob.dissolution) # Flip a weighted coin for each remaning relationship
edges.to.break <- which(edges.coinflip==1) # Make vector of those edges to break
if (length(edges.to.break>0)) { # If there are any edges to break,
edgelist <- edgelist[-edges.to.break,] # then break them.
}
# Add new edges
num.ties.to.add <- expected.edges - dim(edgelist)[1] # Determine how many edges to add. (We assume same as # broken in this simple model.)
if (num.ties.to.add>0) { # If there are edges to add,
if (force.feml.monog==F) { # and if we are *not* enforcing female monogamy,
f <- sample(1:n.femls,num.ties.to.add,replace=T) # sample from the female IDs with replacement,
} else { # else if we are enforcing female monogamy,
femls.with.no.ties <- setdiff(1:n.femls,edgelist$f) # determine which women do not currently have a relationship, and
f <- sample(femls.with.no.ties,num.ties.to.add,replace=F) # sample from them without replacement
}
if (force.male.monog==F) { # and if we are *not* enforcing male monogamy,
m <- sample(1:n.males,num.ties.to.add,replace=T) # sample from the male IDs with replacement,
} else { # else if we are enforcing male monogamy,
males.with.no.ties <- setdiff(1:n.males,edgelist$m) # determine which men do not currently have a relationship, and
m <- sample(males.with.no.ties,num.ties.to.add,replace=F) # sample from them without replacement
}
new.edges <- data.frame(f,m) # Combine the new males and females into pairs in a data frame
edgelist <- data.frame(rbind(edgelist,new.edges)) # "Bind" that data frame onto the end of our existing esgelist
row.names(edgelist) <- 1:dim(edgelist)[1] # Rename the rows of the updated edgelist so they are consecutive numbers again.
}
# Insert new births
new.feml.ids <- femls.dying.of.AIDS # Make a vector of the IDs of the women who just died. We shall "replace" them with a new arrival, by
femls$hiv.status[new.feml.ids] <- 0 # setting the HIV status of that node back to 0, and
femls$inf.time[new.feml.ids] <- NA # setting the infection time for that node back to NA.
new.male.ids <- males.dying.of.AIDS # All parallel with males.
males$hiv.status[new.male.ids] <- 0
males$inf.time[new.male.ids] <- NA
# Track prevalence
feml.prev[time] <- mean(femls$hiv.status) # Calculate current female prevalence and append it onto the female prevalence vector.
male.prev[time] <- mean(males$hiv.status) # Calculate current male prevalence and append it onto the male prevalence vector.
}
windows(20,10) # These two lines set up a plotting window
par(mfrow=c(1,2)) # with two panels.
plot(feml.prev,main= "Female prevalence over time ",ylim=c(0,0.5)) # Plot female prevalence over time
plot(male.prev,main= "Male prevalence over time ",ylim=c(0,0.5)) # Plot male prevalence over time