Model Reference

Joel Hellewell

2024-07-10

Human states

The equations for the human compartments are in section 1.1.1 of Protocol S1 of the original Griffin paper.

In the odin model each of the human compartments (S, T, D, A, U, P) have three dimensions: age, biting heterogeneity category, and intervention category. The differential equations for the first age category are slightly different because they need to include people being born into them.

Susceptible

deriv(S[1, 1:nh, 1:num_int]) <- -FOI[i,j,k]*S[i,j,k] + rP*P[i,j,k] + rU*U[i,j,k] +
  cov[k]*eta*H*het_wt[j] - (eta+age_rate[i])*S[i,j,k]
deriv(S[2:na, 1:nh, 1:num_int]) <- -FOI[i,j,k]*S[i,j,k] + rP*P[i,j,k] + rU*U[i,j,k] -
  (eta+age_rate[i])*S[i,j,k] + age_rate[i-1]*S[i-1,j,k]

$\Lambda \mapsto$FOI

$r_P \mapsto$rP

$r_U \mapsto$rU

$\zeta \mapsto$het_wt are the weights from the biting heterogeneity groups, see Section 1.1.2 of Protocol S1.

The parameter eta ($\eta$) is the constant background death rate, cov is a vector of the proportion of the population in each intervention category (see interventions section later in this reference), age_rate is equal to the inverse of the size of the age group (see equilibrium solution), and H is the sum of all human compartments.

Treated

deriv(T[1, 1:nh, 1:num_int]) <- ft*clin_inc[i,j,k] - rT*T[i,j,k] -
  (eta+age_rate[i])*T[i,j,k]
deriv(T[2:na, 1:nh, 1:num_int]) <- ft*clin_inc[i,j,k] - rT*T[i,j,k] -
  (eta+age_rate[i])*T[i,j,k] + age_rate[i-1]*T[i-1,j,k]

$f_T \mapsto$ft (proportion of clinical disease that is treated)

$r_T \mapsto$rt

For clin_inc see section 1.7 of this guide.

(Clinical) Disease

deriv(D[1, 1:nh, 1:num_int]) <- (1-ft)*clin_inc[i,j,k] - rD*D[i,j,k] -
  (eta+age_rate[i])*D[i,j,k]
deriv(D[2:na, 1:nh, 1:num_int]) <- (1-ft)*clin_inc[i,j,k] - rD*D[i,j,k] -
  (eta+age_rate[i])*D[i,j,k] + age_rate[i-1]*D[i-1,j,k]

$r_D \mapsto$rD

Asymptomatic Disease

deriv(A[1, 1:nh, 1:num_int]) <- (1-phi[i,j,k])*FOI[i,j,k]*Y[i,j,k] - FOI[i,j,k]*A[i,j,k] +
  rD*D[i,j,k] - rA*A[i,j,k] - (eta+age_rate[i])*A[i,j,k]
deriv(A[2:na, 1:nh, 1:num_int]) <- (1-phi[i,j,k])*FOI[i,j,k]*Y[i,j,k] - FOI[i,j,k]*A[i,j,k] +
  rD*D[i,j,k] - rA*A[i,j,k] - (eta+age_rate[i])*A[i,j,k] + age_rate[i-1]*A[i-1,j,k]

$r_A \mapsto$rA

The parameter FOI is explained in Section 2.10

Sub-patent Disease

deriv(U[1, 1:nh, 1:num_int]) <- rA*A[i,j,k] - FOI[i,j,k]*U[i,j,k] - rU*U[i,j,k] -
  (eta+age_rate[i])*U[i,j,k]
deriv(U[2:na, 1:nh, 1:num_int]) <- rA*A[i,j,k] - FOI[i,j,k]*U[i,j,k] - rU*U[i,j,k] -
  (eta+age_rate[i])*U[i,j,k] + age_rate[i-1]*U[i-1,j,k]

$r_U \mapsto$rU

Prophylaxis

deriv(P[1, 1:nh, 1:num_int]) <- rT*T[i,j,k] - rP*P[i,j,k] - (eta+age_rate[i])*P[i,j,k]
deriv(P[2:na, 1:nh, 1:num_int]) <- rT*T[i,j,k] - rP*P[i,j,k] - (eta+age_rate[i])*P[i,j,k] +
  age_rate[i-1]*P[i-1,j,k]

Clinical Incidence

# The number of individuals able to acquire clinical malaria
dim(Y) <- c(na,nh,num_int)
Y[1:na, 1:nh, 1:num_int] <- S[i,j,k]+A[i,j,k]+U[i,j,k]

# The number of new cases at this timestep
dim(clin_inc) <- c(na,nh,num_int)
clin_inc[1:na, 1:nh, 1:num_int] <- phi[i,j,k]*FOI[i,j,k]*Y[i,j,k]

An intermediate step used in the differential equations for the T, D, and A compartments. The parameter Y is just a sum of the S, A, and U compartments for each biting and age combination.

$\phi \Lambda(S(t) + A(t) + U(t)) \mapsto$clin_inc

The parameter $\phi$ (phi) is the probability of clinical disease upon receiving an infectious bite and is explained in Section 2.5.

Immunity states

The mathematical equations for (most of) the human immunity states are described in section 1.1.3 of Protocol S1 of the original Griffin paper.

There are several different types of immunity. In the original Griffin 2010 paper:

$I_B \mapsto$IB (immunity that prevents any infection after bite)

$I_C \mapsto$IC (immunity that reduces chance of clinical disease upon infection)

$I_CA \mapsto$ICA (clinical immunity due to exposure)

$I_CM \mapsto$ICM (clinical immunity due to maternal protection)

$I_A$ (immunity that reduces length of patent infection, no longer in model)

In the supplementary materials for Griffin Nat Comms 2014, $I_A$ is replaced with “detection immunity”, $I_D$, that reduces the probability of case detection and reduced infectiousness to mosquitoes (because of reduced parasite density).

x_I

x_I[] <- user() # intermediate variable for calculating immunity functions

This is a cryptically named parameter that is a constant defined in the equilibrium solution:

# From https://github.com/mrc-ide/dide-deterministic-malaria-model/blob/82069175e66a9ba199aca5fd97172d7b09a4ff1f/R/equilibrium-init-create.R#L68
 den <- 1/(1 + age_rate[1]/mpl$eta)
  for (i in 1:(na-1))
  {
    den[i+1] <- age_rate[i] * den[i]/(age_rate[i+1] + mpl$eta)  # proportion in each age_vector group
  }


# From https://github.com/mrc-ide/dide-deterministic-malaria-model/blob/82069175e66a9ba199aca5fd97172d7b09a4ff1f/R/equilibrium-init-create.R#L97
x_I <- den[1]/mpl$eta
  for (i in 2:na)
  {
    x_I[i] <- den[i]/(den[i - 1] * age_rate[i - 1])  #temporary variables
  }

The best way to understand it is as a way of controlling how aging and death change the levels of immunity. If your age groups are thinner (i.e. you have 2-4yo instead of 2-20yo) then immunity will be lost faster if transmission stops because the people with immunity from before transmission stopped will age out of that group faster. There is a different formula with helpful comments in the malariaEquilibrium package developed by Giovanni and Bob.

Maternally-acquired Immunity

init_ICM_pre[1:nh,1:num_int] <- PM*(ICA[age20l,i,j] + age_20_factor*(ICA[age20u,i,j]-ICA[age20l,i,j]))

deriv(ICM[1, 1:nh, 1:num_int]) <- -1/dCM*ICM[i,j,k] + (init_ICM_pre[j,k]-ICM[i,j,k])/x_I[i]
deriv(ICM[2:na, 1:nh, 1:num_int]) <- -1/dCM*ICM[i,j,k] - (ICM[i,j,k]-ICM[i-1,j,k])/x_I[i]

The level of maternal immunity is assumed at birth to be a proportion ($P_{CM} \mapsto$PM) of the clinical immunity of a 20 year old ($I_{C}(20, t) \mapsto$ICA[age20l, i, j]). Maternal immunity decays at a constant rate ($d_M \mapsto$dCM)

Exposure-driven Immunity

deriv(ICA[1, 1:nh, 1:num_int]) <- FOI[i,j,k]/(FOI[i,j,k] * uCA + 1) - 1/dCA*ICA[i,j,k] -ICA[i,j,k]/x_I[i]
deriv(ICA[2:na, 1:nh, 1:num_int]) <- FOI[i,j,k]/(FOI[i,j,k] * uCA + 1) - 1/dCA*ICA[i,j,k] - (ICA[i,j,k]-ICA[i-1,j,k])/x_I[i]

In the original Griffin paper they consider four different immunity profile functions ($h(X)$). The odin model uses “Model 3” as described in Section 1.1.3 of Protocol S1 for exposure-driven immunity (take note that different immunity profiles are used for other types of immunity in the model!).

Therefore:

$h(\Lambda) = \frac{\Lambda}{\gamma\Lambda + 1} \mapsto$FOI[i,j,k]/(FOI[i,j,k] * uCA + 1)

where the value of $\gamma$ is uCA and controls the delay between infection and exposure-driven immunity increasing.

Clinical Immunity

IC[,,] <- ICM[i,j,k] + ICA[i,j,k]

Clinical immunity combines maternal and exposure-driven immunity.

Phi

Phi controls the probability of developing clinical disease for a given level of clinical immunity. The latest equation for $\phi$ is equation 7 in the supplementary materials for Griffin Nat Comms 2014.

phi[1:na,1:nh,1:num_int] <- phi0*((1-phi1)/(1+(IC[i,j,k]/IC0)^kC) + phi1)

$\phi_0 \mapsto$phi0

$\phi_1 \mapsto$phi1

$I_{C0} \mapsto$IC0

$\kappa_C \mapsto$kC

Infection-blocking Immunity

deriv(IB[1, 1:nh, 1:num_int]) <- EIR[i,j,k]/(EIR[i,j,k]* uB + 1) - IB[i,j,k]/dB - IB[i,j,k]/x_I[i]
deriv(IB[2:na, 1:nh, 1:num_int]) <- EIR[i,j,k]/(EIR[i,j,k]* uB + 1) - IB[i,j,k]/dB - (IB[i,j,k]-IB[i-1,j,k])/x_I[i]

The EIR at a given time for a given age group ($\epsilon(a, t)$) is in section 2.11.

Similar to exposure-driven immunity (Section 2.3) the increase in immunity uses a hill function: $h(\epsilon) = \frac{\epsilon}{\gamma\epsilon + 1}$

where here $\gamma$ is uB.

$d_B \mapsto$dB

b

b is the probability of any disease from an infectious bite, controlled by infection blocking immunity. It is a function that turns the values in IB into a probability. The latest equation for b is equation 6 in the supplementary materials for Griffin Nat Comms 2014.

b[1:na, 1:nh, 1:num_int] <- b0 * ((1-b1)/(1+(IB[i,j,k]/IB0)^kB)+b1)

$b_{0} \mapsto$b0

$b_{1} \mapsto$b1

$I_{B0} \mapsto$IB0

$\kappa_B \mapsto$kB

Detection Immunity

Detection immunity is introduced in equation 5 in the supplementary materials for Griffin Nat Comms 2014.

deriv(ID[1, 1:nh, 1:num_int]) <- FOI[i,j,k]/(FOI[i,j,k]*uD + 1) - ID[i,j,k]/dID - ID[i,j,k]/x_I[i]
deriv(ID[2:na, 1:nh, 1:num_int]) <- FOI[i,j,k]/(FOI[i,j,k]*uD + 1) - ID[i,j,k]/dID - (ID[i,j,k]-ID[i-1,j,k])/x_I[i]

Again a hill function is used (see Section 2.3) with $\gamma$ as uD.

$d_{ID} \mapsto$dID

Probability of Detection

fd[1:na] <- 1-(1-fD0)/(1+(age[i]/aD)^gammaD)
dim(p_det) <- c(na,nh,num_int)
p_det[,,] <- d1 + (1-d1)/(1 + fd[i]*(ID[i,j,k]/ID0)^kD)

The probability that an asymptomatic infection is detected by microscopy is $q \mapsto$p_det (equation 8 in supplementary materials for Griffin Nat Comms 2014).

The function $f_{D} \mapsto$fD is defined just below it.

$f_{D0} \mapsto$fD0

$a_D \mapsto$aD

$\gamma_D \mapsto$gammaD

$d_1 \mapsto$d1

$I_{D0} \mapsto$ID0

$\kappa_D \mapsto$kD

Force of Infection

# Force of infection, depends on level of infection blocking immunity
dim(FOI_lag) <- c(na,nh,num_int)
FOI_lag[1:na, 1:nh, 1:num_int] <- EIR[i,j,k] * (if(IB[i,j,k]==0) b0 else b[i,j,k])

# Current FOI depends on humans that have been through the latent period
dE <- user() # latent period of human infection.
dim(FOI) <- c(na,nh,num_int)
FOI[,,] <- delay(FOI_lag[i,j,k],dE)

EIR is defined in Section 2.11

EIR * b is the definition for the force of infection ($\Lambda(a, t)$) given in equation 2 of supplementary materials for Griffin Nat Comms 2014. The model lags the value of FOI_lag by the value dE (human latent period) and stores it in FOI.

EIR

EIR[,,] <- av_human[k] * rel_foi[j] * foi_age[i] * Iv/omega

Biting and age FOI heterogeneity: (Section 1.1.2 of Protocol S1)

$\zeta \mapsto$rel_foi

$\psi(a) \mapsto$foi_age

Mosquito model parameters: (Section 1.1.4 of Protocol S1)

$I^{v}_{M} \mapsto$Iv (infectious mosquito population, see Section 4.4)

$\alpha^{v} \mapsto$av_human (see Section 6.3)

$\omega \mapsto$omega

where $\omega = \int_{0}^{\infty} \eta(a)\psi(a) da$ is a normalising constant calculated in the equilibrium solution.

Seasonality function

theta2 <- if(ssa0 == 0 && ssa1  == 0 && ssa2  == 0 && ssb1  == 0 && ssb2  == 0 && ssb3  == 0 && theta_c  == 0)
  1 else max((ssa0+ssa1*cos(2*pi*t/365)+ssa2*cos(2*2*pi*t/365)+ssa3*cos(3*2*pi*t/365)+ssb1*sin(2*pi*t/365)+ssb2*sin(2*2*pi*t/365)+ ssb3*sin(3*2*pi*t/365) ) /theta_c,0.001)

The equation for the seasonality fourier function can be found in Section L of Winskill 2017 - S1 Appendix.

This function is a fourier series of the rainfall profile for a given admin 2 region. The seasonality function changes the carrying capacity of the mosquito population over time relative to a baseline carrying capacity (K0 - Section 5.2) that is calculated in the equilibrium solution and can be found in the code and in the supplementary material of White 2011 vector model (equation 7).

Mosquito States

Human infectiousness

How infectious humans with asymptomatic infections are to biting mosquitoes ($c_A$) can vary, whereas the same values for those with clinical ($c_D$) or sub-patent ($c_U$) infections are constant. Equations are in the section “Infectiouness to Mosquitoes” supplementary materials for Griffin Nat Comms 2014.

cA[,,] <- cU + (cD-cU)*p_det[i,j,k]^gamma1

$c_U \mapsto$cU

$c_D \mapsto$cD

$q \mapsto$p_det

$\gamma_I \mapsto$gamma1

Mosquito Force of Infection

FOIvijk[1:na, 1:nh, 1:num_int] <- (cT*T[i,j,k] + cD*D[i,j,k] + cA[i,j,k]*A[i,j,k] + cU*U[i,j,k]) * rel_foi[j] * av_mosq[k]*foi_age[i]/omega
lag_FOIv=sum(FOIvijk)

FOIv <- delay(lag_FOIv, delayGam)

$\Lambda^{v}_{M} \mapsto$FOIvijk

An equation for mosqutio FOI is given in Section 1.1.4 of Protocol S1. However, at the end of Section 2.1 in Protocol S2 - Intervention model of the original Griffin 2010 paper there is a discussion about how using IRS changes the probability that a mosquito will bite a human to be different from the probability that a mosquito will bite and survive (since it may now land on a wall and die after feeding). This means that we calculate the human EIR using $\alpha^v \mapsto$av_human (because all infectious bites on humans are relevant even if the mosquito dies) and the mosquito FOI using $\alpha^v \mapsto$av_mosq (because we only care about mosquitoes that bite, become infected, and don’t die for FOI from the mosquito perspective).

See Section 6.3 for the code for av_mosq and av_human.

The FOI is lagged to account for the extrinsic incubation period of the parasite in the mosquito. tel ## Mosquito incidence

surv <- exp(-mu*delayMos)
ince <- FOIv * Sv
lag_incv <- ince * surv
incv <- delay(lag_incv, delayMos)

The probability that a mosquito survives from acquiring infection until sporozoite emergence is given in Section 1.1.4 of Protocol S1 as $P^{v}_{M} = exp(-\mu\tau_M)$ where:

$P^{v}_{M} \mapsto$surv

$\mu \mapsto$mu

$\tau_M \mapsto$delayMos

In the odin model the parameter incv calculates the incidence of infection in susceptible mosquitoes after accounting for the latent period and mosquito mortality during this time.

Mosquito compartments

# Number of mosquitoes born (depends on PL, number of larvae), or is constant outside of seasonality
betaa <- 0.5*PL/dPL
#betaa <- mv0 * mu0 * theta2

deriv(Sv) <- -ince - mu*Sv + betaa
deriv(Ev) <- ince - incv - mu*Ev
deriv(Iv) <- incv - mu*Iv

# Total mosquito population
mv = Sv+Ev+Iv

Mosquito equations are detailed in Section 1.1.4 of Protocol S1

Larval states

The Larval model was introduced in White P&V 2011.

Some of the notation in the larval model overlaps with the rest of the transmission model so the names are different in the odin code.

Egg oviposition rate

eov <- betaL / mu * (exp(mu/fv)-1)
beta_larval <- eov*mu*exp(-mu/fv)/(1-exp(-mu/fv)) # Number of eggs laid per day
b_lambda <- (gammaL*muLL/muEL-dEL/dLL+(gammaL-1)*muLL*dEL)
lambda <- -0.5*b_lambda + sqrt(0.25*b_lambda^2 + gammaL*beta_larval*muLL*dEL/(2*muEL*mu0*dLL*(1+dPL*muPL)))

$\omega \mapsto$lambda

($\omega$ is in equation 6 of section 3 of the supplementary material for White P&V 2011)

b_lambda is a smaller section of equation 6 for $\omega$ equal to

$-\frac{1}{2}(\gamma\frac{\mu^{0}_L}{\mu^{0}_E} - \frac{d_E}{d_L} + (\gamma - 1)\mu^{0}_L d_E)$

$\gamma \mapsto$gammaL

$\mu^{0}_L \mapsto$muLL

$d_L \mapsto$dLL

$d_E \mapsto$dEL

$\mu^{0}_E \mapsto$muEL

$\beta \mapsto$beta_larval (average eggs laid per day, determines new larvae production)

The parameter called eov requires some reverse engineering:

$eov = \frac{\beta_{L}(exp(\delta\mu_{M}) - 1)}{\mu_{M}}$

where $\mu_M$ is the mosquito mortality rate (mu) and $\delta = \frac{1}{f_v}$ is the days between oviposition (fv is the feeding rate in the model so 1/fv gives time between feeds and therefore oviposition after feeding).

From equation 3 in the methods section of White P&V 2011 we know that:

$\beta_{max} = \frac{\epsilon_{max}\mu_M}{exp(\delta\mu_M - 1)}$

which can be re-arranged to give

$\epsilon_{max} = \frac{\beta_{max} exp(\delta\mu_M - 1)}{\mu_M}$

Therefore we can work out that $\epsilon_{max} \mapsto$eov and $\beta_{max} \mapsto$betaL.

Carrying capacity

K0 <- 2*mv0*dLL*mu0*(1+dPL*muPL)*gammaL*(lambda+1)/(lambda/(muLL*dEL)-1/(muLL*dLL)-1)

# Seasonal carrying capacity KL = base carrying capacity K0 * effect for time of year theta:
KL <- K0*theta2

The equation $K \mapsto$K0 is equation 7 in section 3 of the supplementary material for White P&V 2011.

KL is the larval carrying capacity which is scaled up according to seasonality by theta2 (see Section 3 above)

Birth and death rate

fv <- 1/( tau1/(1-zbar) + tau2 ) # mosquito feeding rate (zbar from intervention param.)
mu <- -fv*log(p1*p2) # mosquito death rate

See Section 4 of the supplementary material for White P&V 2011.

$f \mapsto$fv

$\mu_{M} \mapsto$mu

$\tau_1 \mapsto$tau1 (Time spent seeking blood meal during gonotrophic cycle)

$\tau_2 \mapsto$tau2 (Time spent resting during gonotrophic cycle)

$p_1 \mapsto$p10 Probability of surviving blood meal with no net coverage (the parameter p1 in the model is the probability of surviving a blood meal for a given intervention coverage - see Section 6.3)

$p_2 \mapsto$p2 Probability of surviving one resting cycle (does not change with intervention coverage)

$z \mapsto$zbar average probability of trying again after unsuccessful feeding attempt (see equation 11 in section 4 of the supplementary material for White P&V 2011).

Larval compartments

# (beta_larval (egg rate) * total mosquito) - den. dep. egg mortality - egg hatching
deriv(EL) <- beta_larval*mv-muEL*(1+(EL+LL)/KL)*EL - EL/dEL
# egg hatching - den. dep. mortality - maturing larvae
deriv(LL) <- EL/dEL - muLL*(1+gammaL*(EL + LL)/KL)*LL - LL/dLL
# pupae - mortality - fully developed pupae
deriv(PL) <- LL/dLL - muPL*PL - PL/dPL

Larval compartments are described in White P&V 2011.

Interventions

The intervention model details are described in Protocol S2 of the original Griffin 2010 paper.

Mortality and repellency

# Calculates decay for ITN/IRS
ITN_decay = if(t < ITN_IRS_on) 0 else exp(-((t-ITN_IRS_on)%%ITN_interval) * itn_loss)
IRS_decay = if(t < ITN_IRS_on) 0 else exp(-((t-ITN_IRS_on)%%IRS_interval) * irs_loss)

# The r,d and s values turn on after ITN_IRS_on and decay accordingly
d_ITN <- if(t < ITN_IRS_on) 0 else d_ITN0*ITN_decay
r_ITN <- if(t < ITN_IRS_on) 0 else r_ITN1 + (r_ITN0 - r_ITN1)*ITN_decay
s_ITN <- if(t < ITN_IRS_on) 1 else 1 - d_ITN - r_ITN

r_IRS <- if(t < ITN_IRS_on) 0 else r_IRS0*IRS_decay
d_IRS <- if(t < ITN_IRS_on) 0 else chi*d_IRS0*IRS_decay
s_IRS <- if(t < ITN_IRS_on) 1 else 1 - d_IRS

ITN_decay and IRS_decay control the decay in the effectiveness of ITNs and IRS in terms of mortality and repelling ability.

See Section 2.1.1.1 (pgs 8 & 9) in Protocol S2 for the equations for the parameters below:

$\gamma_N \mapsto$itn_loss

$r_{N0} \mapsto$r_ITN0

$r_{NM} \mapsto$r_ITN1

$r_N \mapsto$r_ITN

$d_{N0} \mapsto$d_ITN0

$d_N \mapsto$d_ITN

$s_N \mapsto$s_ITN

$r_{S0} \mapsto$r_IRS0

$\gamma_S \mapsto$irs_loss

$r_S \mapsto$r_IRS

$\chi \mapsto$chi

$d_{S0} \mapsto$d_IRS0

$d_S \mapsto$d_IRS

$s_S \mapsto$s_IRS

Probabilities of biting and survival

See Table S2.2 in Protocol S2 of the original Griffin 2010 paper.

# probability that mosquito bites and survives for each intervention category
dim(w_) <- 4
w_[1] <- 1
w_[2] <- 1 - bites_Bed + bites_Bed*s_ITN
w_[3] <- 1 - bites_Indoors + bites_Indoors*(1-r_IRS)*s_IRS
w_[4] <- 1 - bites_Indoors + bites_Bed*(1-r_IRS)*s_ITN*s_IRS + (bites_Indoors - bites_Bed)*(1-r_IRS)*s_IRS
w[] <- w_[i]
dim(w) <- num_int

# probability that mosq feeds during a single attempt for each int. cat.
dim(yy_) <- 4
yy_[1] <- 1
yy_[2] <- w_[2]
yy_[3] <- 1 - bites_Indoors + bites_Indoors*(1-r_IRS)
yy_[4] <- 1 - bites_Indoors + bites_Bed*(1-r_IRS)*s_ITN + (bites_Indoors - bites_Bed)*(1-r_IRS)
yy[] <- yy_[i]
dim(yy) <- num_int

# probability that mosquito is repelled during a single attempt for each int. cat.
dim(z_) <- 4
z_[1] <- 0
z_[2] <- bites_Bed*r_ITN
z_[3] <- bites_Indoors*r_IRS
z_[4] <- bites_Bed*(r_IRS+ (1-r_IRS)*r_ITN) + (bites_Indoors - bites_Bed)*r_IRS
z[] <- z_[i]
dim(z) <- num_int

$\Phi_I \mapsto$bites_Indoors

$\Phi_B \mapsto$bites_Bed

Average biting rates

See pages 5-6 in Protocol S2 of the original Griffin 2010 paper.

# Calculating Z (zbar) and W (wbar) - see Supplementary materials 2 for details
dim(zhi) <- num_int
dim(whi) <- num_int
zhi[1:num_int] <- cov[i]*z[i]
whi[1:num_int] <- cov[i]*w[i]
zh <- if(t < ITN_IRS_on) 0 else sum(zhi)
wh <- if(t < ITN_IRS_on) 1 else sum(whi)
# Z (zbar) - average probability of mosquito trying again during single feeding attempt
zbar <- Q0*zh
# W (wbar) - average probability of mosquito successfully feeding during single attempt
wbar <- 1 - Q0 + Q0*wh

# p1 is the updated p10 given that interventions are now in place:
p1 <- wbar*p10/(1-zbar*p10)
Q <- 1-(1-Q0)/wbar # updated anthropophagy given interventions
av <- fv*Q # biting rate on humans
dim(av_mosq) <- num_int
av_mosq[1:num_int] <- av*w[i]/wh # rate at which mosquitoes bite each int. cat.
dim(av_human) <- num_int
av_human[1:num_int] <- av*yy[i]/wh # biting rate on humans in each int. cat.

$Z \mapsto$zbar

$W \mapsto$wbar

$Q^v \mapsto$Q

$p_1 \mapsto$p1

The average biting rate ($\alpha_v \mapsto$av_mosq or av_human) is calculated using w or yy depending on whether it matters that the mosquito may die after feeding. For the FOI on mosquitoes, we only want to count infections in mosquitoes that don’t then die due to IRS. For the EIR for humans, infectious bites count even if mosquitoes die due to IRS afterwards.