https://github.com/morgankain/RRV_HostVectorCompetence
Tip revision: be7e87c3c4c8af0420a8dd42cdcff5586fdbad90 authored by Morgan Kain on 25 May 2021, 16:23:11 UTC
Merge pull request #1 from morgankain/add-license-1
Merge pull request #1 from morgankain/add-license-1
Tip revision: be7e87c
16_R0_calc.R
#####
## Calculate R0 and all metrics of host competence using model components
#####
## Grab the parameters and the output from the model components
num_hosts <- dim(h_to_m_trans_all_samps_adj_to_com)[2]
num_mosq <- dim(h_to_m_trans_all_samps_adj_to_com)[3]
prop_hosts <- host_prop_for_R0_all_samps
inf_days <- seq(1, length(inf_days), 1)
prop_mosq <- mosq_prop_for_R0$prop
daily_bites <- rep(bite_rate, num_mosq)
## Bump up Aedes aegypti to take into account the multiple feedings within a gonotrophic cycle (if Ae aegypti is present in the dataset)
if (length(grep("aegypti", dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]])) > 0) {
daily_bites[grep("aegypti", dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]])] <-
daily_bites[grep("aegypti", dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]])] * 2
}
mos_surv <- mosq_surv_for_R0_all_samps_adj_to_com
mosq_bite_pref <- mosq_bite_pref_all_samps_adj_to_com
R0.out <- data.frame(
model = "Community_with_uncertainty"
, R0_primary = numeric(n_samps)
, foi_on_A = numeric(n_samps)
, foi_on_H = numeric(n_samps)
, rel_foi_on_H = numeric(n_samps)
, foi_from_A = numeric(n_samps)
, foi_from_H = numeric(n_samps)
, rel_foi_from_H = numeric(n_samps))
## Setup mosquito biting preference with and without including host abundance.
## Note: mosquito biting _preference_ is modeled as their intrinsic propensity to bite a given species, which is given as their proportional
## increase or decrease relative to 1 which designates no preference and an equivalence to biting that species randomly (in proportion to the
## relative abundance of that species in the community). These intrinsic preferences are then multiplied by the actual observed abundance and
## scaled to a proportion. Note that this can be calculated assuming that all hosts are equally abundant, but this isn't very ecologically interesting
## because it isn't a real scenario
mosq_bite <- array(dim = dim(mosq_bite_pref))
for (k in 1:n_samps) {
mosq_bite[,,k] <- mosq_bite_pref[,,k] *
{
if (use.host_abundance) {
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(prop_hosts[, k], num_mosq), byrow = TRUE)
} else {
matrix(nrow = num_mosq, ncol = num_hosts, data = 1/num_hosts, byrow = TRUE)
}
}
mosq_bite[,,k] <- sweep(mosq_bite[,,k], 1, rowSums(mosq_bite[,,k]), "/")
}
## Total number of new hosts infected by a mosquito over its lifespan...
mosq_trans_piece <- colSums(
## Set mosquito survival to vary by species ore set all species to the mean
if (use.mosq_survival) {
mosq_surv_for_R0_all_samps_adj_to_com * m_to_h_trans_all_samps_adj_to_com
} else {
array(dim = dim(mosq_surv_for_R0_all_samps_adj_to_com), data = rowMeans(mosq_surv_for_R0_all_samps_adj_to_com[,,1])) * m_to_h_trans_all_samps_adj_to_com
}
)
## ... which also includes a mosquitoes daily bite rate...
mosq_trans_piece <- mosq_trans_piece * matrix(ncol = n_samps, nrow = num_mosq, data = rep(daily_bites, n_samps))
## ... and the status of the hosts that the mosquitoes are biting:
## WAIFW_left is I mosquitoes infecting S hosts. Can consider mosquito biting preference, or some subset of that such as host abundance, or none of the above
WAIFW_left <- array(data = NA, dim = c(num_hosts, num_mosq, n_samps))
for (k in 1:n_samps) {
## Proportion of these hosts susceptible. From the disease's perspective a proportion of bites are "wasted" (don't lead to a new infection) in the
## proportion of hosts that are seropositive already
if (use.mosq_bite_pref) {
WAIFW_left[,,k] <- t(
{
if (use.host_seroprev) {
## Here use: mosquito biting preference (already scaled or not by host abundance) scaled by host seropositivity
sweep(mosq_bite[,,k], 2, (1 - host_sero$prop_positive), "*") *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
} else {
## Or just mosquito biting preference
mosq_bite[,,k] *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
}
}
)
} else {
WAIFW_left[,,k] <- t(
{
if (use.host_abundance) {
if (use.host_seroprev) {
## Here use host abundance, but assume mosquito biting is random in the community (bites only determined by host abundance)
sweep(t(matrix(nrow = num_hosts, ncol = num_mosq, data = prop_hosts[, k])), 2, (1 - host_sero$prop_positive), "*") *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
} else {
## Same thing, just don't weight by seropositivity
t(matrix(nrow = num_hosts, ncol = num_mosq, data = prop_hosts[, k])) *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
}
} else {
if (use.host_seroprev) {
sweep(t(matrix(nrow = num_hosts, ncol = num_mosq, data = 1/num_hosts)), 2, (1 - host_sero$prop_positive), "*") *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
} else {
t(matrix(nrow = num_hosts, ncol = num_mosq, data = 1/num_hosts)) *
matrix(nrow = num_mosq, ncol = num_hosts, data = rep(mosq_trans_piece[, k], num_hosts))
}
}
}
)
}
}
## WAIFW_right is I hosts infecting S mosquitoes. Can consider just titer and infection probability or scale by the abundance of hosts and mosquitoes
## In brief to explain the logic here: since we have number of mosquitoes per host in the community, one infected host among many of its kind is assumed to be bit
## at random by X mosquitoes per each of this host type (scaled by the mosquitoes biting preference etc.). This mosquito/host ratio trick and assumption of homogeneous
## mixing of the mosquito population and no preferential biting on I or S hosts translates one infected host to the total number of infected mosquitoes per day or
## over that infected hosts infectious period
WAIFW_right_s1 <- apply(h_to_m_trans_all_samps_adj_to_com, 3:4, FUN = colSums) ## This can be interpreted as one mosquito of each type biting each host once per day, how many total mosquitoes of each type would be infected
WAIFW_right_s2 <- array(dim = dim(WAIFW_right_s1))
WAIFW_right <- array(data = NA, dim = c(num_mosq, num_hosts, n_samps))
#####
## !! For full details on the calculation before see the notes page: --- transmission_notes.R ---
## The calculation below is as in the notes page, order of operations are just a bit different
#####
for (k in 1:n_samps) {
## This step calculates the raw capability of _AN INFECTED_ host of _species X_ to infect each mosquito which is given by:
## The total number of mosquitoes this infected host infects over its infectious period, which is determined by:
## its titer on each day, and the number and identity of mosquitoes biting it each day, which in turn is determined by:
## that mosquitoes capability of picking up infection, the abundance of each mosquito relative to this infected host, and the biting preference of these mosquitoes
## NOTE: If mosquito biting preference is NOT considered, the relative abundance of the species of the infected host is irrelevant for this
## arm of the transmission cycle because N mosquitoes are simply biting randomly. However, if mosquito biting preference is considered (for which
## host abundance does play a role), the single infected host of species X may be bit more or less frequently than random
temp_WAIFW <- {
if (use.mosq_bite_pref) {
WAIFW_right_s1[,,k] *
## mosquito bites on each host per day = the proportion of each of their bites on each host type * their daily biting rate
t(mosq_bite[,,k]) *
matrix(nrow = num_hosts, ncol = num_mosq, data = rep(daily_bites, num_hosts)) *
## Mosquito to host ratio (the number of S mosquitoes biting each infected host at the bite rate).
## Because mosquito _relative abundance_ is used (below), this m_to_h_ratio captures the sum of all mosquito individuals in the community
## relative to the sum of all host individuals in the community
m_to_h_ratio
} else {
WAIFW_right_s1[,,k] *
matrix(nrow = num_hosts, ncol = num_mosq, data = 1/num_hosts) *
matrix(nrow = num_hosts, ncol = num_mosq, data = rep(daily_bites, num_hosts)) *
m_to_h_ratio
}
}
WAIFW_right_s2[,,k] <- {
## mosquito _relative abundance_
if (use.mosq_abundance) {
temp_WAIFW * matrix(nrow = num_hosts, ncol = num_mosq, data = rep(prop_mosq, num_hosts), byrow = TRUE)
} else {
temp_WAIFW * matrix(nrow = num_hosts, ncol = num_mosq, data = 1/num_mosq)
}
}
## If scaling by the proportion of all hosts that show titer, multiply by that ratio to get the average hosts' response (i.e. some hosts are going to lead
## to no infection). By including this here we are assuming essentially that a host of each type becomes _exposed_ but there is some probability that that
## host is actually _infectious_
if (use.cond_titer) {
WAIFW_right_s2[,,k] <- sweep(WAIFW_right_s2[,,k], 1, prop_inf_for_R0$num_inf, "*")
}
WAIFW_right[,,k] <- t(WAIFW_right_s2[,,k])
}
## Rename WAIFW_right as host competence (host-mosquito, the first form of host competence) for use later
host_competence <- WAIFW_right
dimnames(host_competence) <- list(
dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]]
, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]]
, NULL)
## Calculate R0
physiol_mat <- array(dim = c(num_hosts, num_hosts, n_samps))
temp_mmat.f <- array(dim = c(num_hosts, num_hosts, n_samps))
## For host-to-host transmission which we are defining here as the second measure of host competence
FOI_on_t <- matrix(nrow = n_samps, ncol = num_hosts, data = 0)
FOI_from_t <- matrix(nrow = n_samps, ncol = num_hosts, data = 0)
physiol_mat_mm <- array(dim = c(num_mosq, num_mosq, n_samps))
temp_mmat.f_mm <- array(dim = c(num_mosq, num_mosq, n_samps))
FOI_on_t_mm <- matrix(nrow = n_samps, ncol = num_mosq, data = 0)
FOI_from_t_mm <- matrix(nrow = n_samps, ncol = num_mosq, data = 0)
## Method for calculating R0 sticks these matrices in the diagonal of a larger matrix then takes the eigen
## For details on the matrix algebra component here see:
## source("17_matrix_algebra_exploration.R")
for (k in 1:n_samps) {
temp_mat <- matrix(data = 0, nrow = num_hosts + num_mosq, ncol = num_hosts + num_mosq)
temp_mat[(num_hosts+1):(num_hosts + num_mosq), 1:num_hosts] <- WAIFW_right[,,k]
temp_mat[1:(num_hosts), (1+num_hosts):(num_hosts+num_mosq)] <- WAIFW_left[,,k]
R0 <- Re(eigen(temp_mat)$values)
R0 <- R0[which(R0 > 0)]
R0.out[k, 2] <- R0[1]
## And for FOI, the pairwise host-host or mosquito-mosquito transmission capability:
temp_mmat <- WAIFW_left[,,k] %*% WAIFW_right[,,k]
temp_mmat.f[,,k] <- temp_mmat
physiol_mat[,,k] <- temp_mmat
## New infected hosts of each type after a generation of host-mosquito-host. Will be all the same when there is no preference weighting
FOI_on_t[k, ] <- rowSums(temp_mmat)
FOI_from_t[k, ] <- colSums(temp_mmat)
## Full transmission pairwise identity matrix
temp_mmat_mm <- WAIFW_right[,,k] %*% WAIFW_left[,,k]
temp_mmat.f_mm[,,k] <- temp_mmat_mm
physiol_mat_mm[,,k] <- temp_mmat_mm
## FOI on and from using definition 1 from above
FOI_on_t_mm[k, ] <- rowSums(temp_mmat_mm)
FOI_from_t_mm[k, ] <- colSums(temp_mmat_mm)
## The largest FOI value ON and FOI ON humans
R0.out[k, c(3, 4)] <- FOI_on_t[k, ][c(which(FOI_on_t == max(FOI_on_t[k, ]))[1], which(host_prop$host_species == "human"))]
## FOI ON humans relative to the average FOI ON
R0.out[k, 5] <- FOI_on_t[k, ][which(host_prop$host_species == "human")] / mean(FOI_on_t[k, ])
## The largest FOI value FROM and FOI FROM humans
R0.out[k, c(6, 7)] <- FOI_from_t[k, ][c(which(FOI_from_t[k, ] == max(FOI_from_t[k, ]))[1], which(host_prop$host_species == "human"))]
## FOI FROM humans relative to the average FOI FROM
R0.out[k, 8] <- FOI_from_t[k, ][which(host_prop$host_species == "human")] / mean(FOI_from_t[k, ])
}
####
## Calculate host and mosquito competence for just the one-step physiology (hosts ability to infect mosquitoes
## and mosquitoes ability to pick up infection from hosts)
####
source("18_host_mosq_comp_remaining.R")
####
## Summary statistics on FOI_on and FOI_from (which is in fact host-to-host transmission competence)
####
source("19_host_competence_cleanup.R")
source("20_mosq_competence_cleanup.R")
## And for all calculations stick together previous and current estimates
if (complexity_counter == 1) {
####
## Compile summary from hosts
####
## Host competence as a host's ability to infect mosquitoes
host_competence.r.gg.f <- host_competence.r.gg
host_competence.gg.f <- host_competence.gg
## Host competence as a host's ability to complete the life cycle of host-host
host_competence.r.gg.hh.f <- host_competence.r.gg.hh
host_competence.gg.hh.f <- host_competence.gg.hh
## FOI ON (host)
FOI_on_t.gg.f <- FOI_on_t.gg
## Host to host transmission matrix
physiol_mat.gg.f <- physiol_mat.gg
####
## Compile summary from mosquitoes
####
## Mosquito competence as a mosquitoes ability to pick up infection
mosq_competence.r.gg.f <- mosq_competence.r.gg
mosq_competence.gg.f <- mosq_competence.gg
## Mosquito competence as a mosquitoes ability to complete a mosquito-mosquito life cycle
mosq_competence.r.gg.mm.f <- mosq_competence.r.gg.mm
mosq_competence.gg.mm.f <- mosq_competence.gg.mm
## FOI ON (mosq)
FOI_on_t.gg.f_mm <- FOI_on_t.gg_mm
## mosquito to mosquito transmission matrix
physiol_mat.gg.f.mm <- physiol_mat.gg_mm
} else {
host_competence.r.gg.f <- rbind(host_competence.r.gg.f, host_competence.r.gg)
host_competence.gg.f <- rbind(host_competence.gg.f, host_competence.gg)
host_competence.r.gg.hh.f <- rbind(host_competence.r.gg.hh.f, host_competence.r.gg.hh)
host_competence.gg.hh.f <- rbind(host_competence.gg.hh.f, host_competence.gg.hh)
mosq_competence.r.gg.f <- rbind(mosq_competence.r.gg.f, mosq_competence.r.gg.mm)
mosq_competence.gg.f <- rbind(mosq_competence.gg.f, mosq_competence.gg.mm)
mosq_competence.r.gg.mm.f <- rbind(mosq_competence.r.gg.mm.f, mosq_competence.r.gg.mm)
mosq_competence.gg.mm.f <- rbind(mosq_competence.gg.mm.f, mosq_competence.gg.mm)
FOI_on_t.gg.f_mm <- rbind(FOI_on_t.gg.f_mm, FOI_on_t.gg_mm)
FOI_on_t.gg.f <- rbind(FOI_on_t.gg.f, FOI_on_t.gg)
physiol_mat.gg.f <- rbind(physiol_mat.gg.f, physiol_mat.gg)
physiol_mat.gg.f.mm <- rbind(physiol_mat.gg.f.mm, physiol_mat.gg_mm)
}
print(
paste(round(complexity_counter / nrow(model.runs), 2)*100, "%", "Complete", sep = " ")
)