In this vignette we demonstrate how to set up a model incorporating the Le Menach model of ITN based vector control (see this paper). We use the Ross-Macdonald model of adult mosquito dynamics, the SIS model of human population dynamics, and the basic competition model of aquatic mosquito dynamics to fill the dynamical components \(\mathcal{M}, \mathcal{X}, \mathcal{L}\).

## Parameters

This case study will use a simple model with 3 patches, 3 population strata, and 3 aquatic habitats. As usual, we set up parameter values and compute the values of state variables at equilibrium. As part of the equilibrium calculation we must compute \(\Upsilon(0) = \exp\left(-\int_{-\tau}^{0} \Omega(s) ds \right)\); the value of the integrated mosquito demography matrix at the initial time point. To simplify things, we simply assume that conditions were constant prior to \(t=0\) so that \(\Upsilon(0) = e^{-\Omega\tau}\).

```
pars = make_parameters_xde()
pars$nPatches <- 3
pars$nStrata <- 3
pars$nHabitats <- 3
# human parameters
b <- 0.55
c <- 0.15
r <- 1/200
wf <- rep(1, pars$nStrata)
pfpr <- runif(n = pars$nStrata, min = 0.25, max = 0.35)
H <- rpois(n = pars$nStrata, lambda = 1000)
X <- rbinom(n = pars$nStrata, size = H, prob = pfpr)
membershipH = c(1:3)
searchWtsH = rep(1,3)
Psi <- matrix(
data = c(
0.9, 0.05, 0.05,
0.05, 0.9, 0.05,
0.05, 0.05, 0.9
), nrow = pars$nStrata, ncol = pars$nPatches, byrow = T
)
Psi <- t(Psi)
# adult mosquito parameters
f <- 0.3
q <- 0.9
g <- 1/10
sigma <- 1/100
nu <- 1/2
eggsPerBatch <- 30
tau <- 11
# mosquito movement calK
calK <- matrix(0, pars$nPatches, pars$nPatches)
calK[upper.tri(calK)] <- 1/(pars$nPatches-1)
calK[lower.tri(calK)] <- 1/(pars$nPatches-1)
calK <- calK/rowSums(calK)
calK <- t(calK)
Omega <- make_Omega(g = g, sigma = sigma, K = calK, nPatches = pars$nPatches)
Omega_inv <- solve(Omega)
Upsilon <- expm::expm(-Omega * tau)
Upsilon_inv <- expm::expm(Omega * tau)
```

## Equilibrium

Now we compute the equilibrium conditions for the adult mosquitoes and human populations, such that the PfPR in the human population would be maintained at the input levels if conditions were unchanging.

```
# derived EIR to sustain equilibrium pfpr
EIR <- diag(1/b, pars$nStrata) %*% ((r*X) / (H - X))
# ambient pop
W <- Psi %*% H
# biting distribution matrix
beta <- diag(wf) %*% t(Psi) %*% diag(1/as.vector(W), pars$nPatches)
# kappa
kappa <- t(beta) %*% (X*c)
# equilibrium solutions
Z <- diag(1/(f*q), pars$nPatches) %*% ginv(beta) %*% EIR
MY <- diag(1/as.vector(f*q*kappa), pars$nPatches) %*% Upsilon_inv %*% Omega %*% Z
Y <- Omega_inv %*% (diag(as.vector(f*q*kappa), pars$nPatches) %*% MY)
M <- MY + Y
G <- solve(diag(nu+f, pars$nPatches) + Omega) %*% diag(f, pars$nPatches) %*% M
Lambda <- Omega %*% M
```

Given the equilibrium value of \(\Lambda\) required to sustain mosquito populations at such a level sufficient to maintain transmission at that level of PfPR, as well as values for \(\psi\) (maturation rate) and \(\phi\) (density independent mortality), we compute equilibrium values of \(L\) (aquatic mosquito density) and \(\theta\) (density dependent mortality).

```
# aquatic habitat membership matrix (assume 1 habitat per patch)
calN <- matrix(0, pars$nPatches, pars$nHabitats)
diag(calN) <- 1
# egg dispersal matrix (assume 1 habitat per patch)
calU <- matrix(0, pars$nHabitats, pars$nPatches)
diag(calU) <- 1
alpha <- as.vector(solve(calN) %*% Lambda)
psi <- 1/10
phi <- 1/12
eta <- as.vector(calU %*% (G * nu * eggsPerBatch))
L <- alpha/psi
theta <- (eta - psi*L - phi*L)/(L^2)
```

## Setup

Now that all state variables have been solved at equilibrium, we can
set up the parameters and components of the system. We also need to set
up a time-varying function to compute the coverage of ITNs at any time
point. We use a sine curve with a period of 365 days which goes from 0
to 1 over that period, with the phase shifted so that at 0 it returns 0.
The function also returns 0 for any \(t<0\). This must be a function that
takes a single argument `t`

(time) and returns a scalar
value.

We use the null model of human demographic dynamics, which assumes \(H\) is constant for all time.

```
pars$calU <- calU
pars$calN <- calN
# ITN coverage
ITN_cov <- function(t) {ifelse(t < 0, 0, (sin(2*pi*(t-365/4) / 365) + 1) / 2)}
# set parameters
pars = make_parameters_exogenous_null(pars)
pars = make_parameters_vc_lemenach(pars, phi = ITN_cov)
pars = make_parameters_MYZ_RM(pars = pars, g = g,
sigma = sigma, calK = calK,
f = f, q = q, nu = nu,
eggsPerBatch = eggsPerBatch,
tau = tau)
pars = make_inits_MYZ_RM(pars = pars, M0 = M,
Y0 = Y, Z0 = Z,Upsilon0=Upsilon)
pars = make_parameters_L_basic(pars = pars, psi = psi,
phi = phi, theta = theta)
pars = make_inits_L_basic(pars = pars, L0 = L)
pars = make_parameters_X_SIS(pars = pars, b = b, c = c, r = r)
pars = make_inits_X_SIS(pars = pars, X0=X)
pars = make_parameters_demography_null(pars = pars, H=H, membershipH=membershipH,
searchWtsH=searchWtsH, TimeSpent=Psi)
pars = make_parameters_mi_null(pars)
pars = make_indices(pars = pars)
```

## Solve

After setting parameters we must set the initial conditions of the model. Then we’re ready to solve the model! We integrate over a 5 year interval, to show the effects of the “seasonal” ITN coverage.

```
# ICs
y0 <- get_inits(pars)
```

```
# solve the model
out <- dede(y = y0, times = 0:(365*5), func = xDE_diffeqn, parms = pars)
```

We can plot the output to study the effects of the seasonal ITN coverage on the state variables. We see that total and parous mosquito densities \(M\) and \(G\) now vary in a sinusoidal manner as coverage changes, as does \(L\), larval density, as the number of ovipositing mosquitoes changes. Because in addition to increasing mosquito mortality, ITNs also decrease the feeding rate and proportion of bloodmeals taken on humans due to the repellency effect, we see the densities of infected and infectious mosquitoes (\(Y\) and \(Z\)) steadily decrease over time as there are no longer sufficient bloodmeals being taken to sustain prevalence. A similar pattern appears in human prevalence \(X\), with slower dynamics due to the slow recovery rate of untreated infection.

```
colnames(out)[pars$L_ix+1] <- paste0('L_', 1:pars$nHabitats)
colnames(out)[pars$M_ix+1] <- paste0('M_', 1:pars$nPatches)
colnames(out)[pars$G_ix+1] <- paste0('G_', 1:pars$nPatches)
colnames(out)[pars$Y_ix+1] <- paste0('Y_', 1:pars$nPatches)
colnames(out)[pars$Z_ix+1] <- paste0('Z_', 1:pars$nPatches)
colnames(out)[pars$X_ix+1] <- paste0('X_', 1:pars$nStrata)
out <- out[, -c(pars$Upsilon_ix+1)]
out <- as.data.table(out)
out <- melt(out, id.vars = 'time')
out[, c("Component", "Stratification") := tstrsplit(variable, '_', fixed = TRUE)]
out[, variable := NULL]
ggplot(data = out, mapping = aes(x = time, y = value, color = Stratification)) +
geom_line() +
facet_wrap(. ~ Component, scales = 'free') +
theme_bw()
```