```
library(exDE)
library(expm)
library(data.table)
library(ggplot2)
```

We show how to setup, solve, and analyze models of mosquito-borne
pathogen transmission dynamics and control using modular software. This
vignette is designed to explain modular notation by constructing a model
with five aquatic habitats (\(l=5\)),
three patches (\(p=3\)), and four human
population strata (\(n=4\)). We call it
`5-3-4`

.

The model `5-3-4`

is designed to illustrate some important
features of the framework and notation. We assume that:

the first three habitats are found in patch 1; the last two are in patch 2; patch 3 has no habitats.

patch 1 has no residents; patches 2 and 3 are occupied, each with two different population strata;

Transmission among patches is modeled using the concept of

*time spent*, which is similar to the*visitation rates*that have been used in other models. While the strata have a residency (i.e; a patch they spend most of their time in), each stratum allocates their time across all the habitats.

We already know three important parameters, \(l\), \(p\)
and \(n\) because they are determined
in the early stages of model building. The `exDE`

package
expects all parameters to be contained in a list object, containing
`nHabitats`

, `nPatches`

, and `nStrata`

which correspond to `l`

, `p`

and
`n`

.

```
<- list(
params nHabitats = 5,
nPatches = 3,
nStrata = 4
)<- list2env(params) params
```

The aquatic habitat membership matrix \(\mathcal{N}\) is a \(p \times l\) matrix mapping aquatic habitats to the patches which contain them. It should be attached directly to the parameters list.

\[\begin{equation} {\cal N} = \left[ \begin{array}{ccccc} 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 0 & 0\\ \end{array} \right] \end{equation}\]

```
<- matrix(
calN data = c(1,1,1,0,0,
0,0,0,1,1,
0,0,0,0,0),
nrow = params$nPatches, ncol = params$nHabitats, byrow = TRUE
)$calN <- calN params
```

The egg dispersal matrix \(\mathcal{U}\) is a \(l \times p\) matrix describing how eggs laid by adult mosquitoes in a patch are allocated among the aquatic habitats in that patch. It is also attached directly to the parameters list.

\[\begin{equation} {\cal U} = \left[ \begin{array}{ccccc} .7 & 0 & 0\\ .2 & 0 & 0\\ .1 & 0 & 0\\ 0 & .8 & 0\\ 0 & .2 & 0\\ \end{array} \right] \end{equation}\]

```
<- matrix(c(.7, .2, .1, .8, .2), 5, 1)
xi $calU <- t(calN %*% diag(as.vector(xi))) params
```

For this simulation, we use the basic competition model of larval
dynamics (see more here). It requires
specification of three parameters, \(\psi\) (maturation rates), \(\phi\) (density-independent mortality
rates), and \(\theta\)
(density-dependent mortality terms), and initial conditions. The
function `exDE::make_parameters_L_basic`

does basic checking
of the input parameters and returns a list with the correct class for
method dispatch. The returned list is attached to the main parameter
list with name `Lpar`

.

```
<- rep(1, params$nHabitats)
L0 <- rep(1/8, params$nHabitats)
psi <- rep(1/8, params$nHabitats)
phi <- c(1/10, 1/20, 1/40, 1/100, 1/10)
theta
make_parameters_L_basic(pars = params, psi = psi, phi = phi, theta = theta, L0 = L0)
```

We use the ODE version of the generalized Ross-Macdonald model (see
more here). Part of the specification of
parameters includes the construction of the mosquito dispersal matrix
\(\mathcal{K}\), and the mosquito
demography matrix \(\Omega\). Like for
the aquatic parameters, we use
`exDE::make_parameters_MYZ_RM_ode`

to check parameter types
and return a list with the correct class for method dispatch. We attach
the returned list to the main parameter list with name
`MYZpar`

.

```
<- rep(1/12, params$nPatches)
g <- rep(1/12/2, params$nPatches)
sigma <- t(matrix(
calK c(c(0, .6, .3),
c(.4, 0, .7),
c(.6, .4, 0)), 3, 3))
<- rep(1/3, params$nPatches)
f <- rep(0.9, params$nPatches)
q <- c(1/3,1/3,0)
nu <- 30
eggsPerBatch <- 12
tau <- rep(100, params$nPatches)
M0 <- rep(10, params$nPatches)
G0 <- rep(1, params$nPatches)
Y0 <- rep(0, params$nPatches)
Z0
<- make_Omega(g, sigma, calK, params$nPatches)
Omega <- expm::expm(-Omega*tau)
Upsilon
make_parameters_MYZ_RM_ode(pars = params, g = g, sigma = sigma, calK = calK, tau = tau, f = f, q = q, nu = nu, eggsPerBatch = eggsPerBatch, M0 = M0, G0 = G0, Y0 = Y0, Z0 = Z0)
```

Although not directly used in this example, we create the residency membership matrix \(\mathcal{J}\), a \(p \times n\) matrix indicating which patch each human population strata resides in.

\[\begin{equation} {\cal J} = \left[ \begin{array}{cccc} 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ \end{array} \right] \end{equation}\]

```
<- t(matrix(
calJ c(c(0,0,0,0),
c(1,1,0,0),
c(0,0,1,1)), 4, 3
))
```

We then create the time at risk matrix \(\Psi\), a \(p \times n\) matrix describing how each human strata spends their time across patches.

\[\begin{equation} \Psi = \left[ \begin{array}{cccc} 0.01 & .01 & .001 & .001 \\ 0.95 & .92 & .04 & .02 \\ 0.04 & .02 & .959 & .929 \\ \end{array} \right] \end{equation}\]

```
<- t(matrix(
Psi c(c(0.01,0.01,0.001,0.001),
c(.95,.92,.04,.02),
c(.04,.02,.959,.929)), 4, 3
))
```

We use the basic SIS (Susceptible-Infected-Susceptible) model for the
human component (see more here). We set it
up like the rest of the components, using
`exDE::make_parameters_X_SIS`

to make the correct return
type, which is attached to the parameter list with name
`Xpar`

.

```
<- matrix(c(10,90, 100, 900), 4, 1)
H <- as.vector(0.2*H)
X0 <- rep(1/200, params$nStrata)
r <- rep(0.55, params$nStrata)
b <- c(0.1, .02, .1, .02)
c
make_parameters_X_SIS(pars = params, b = b, c = c, r = r, Psi = Psi, X0 = X0, H = as.vector(H))
```

We need to specify models for exogenous forcing and vector control. We simply use the null models for both (no presence of vector control and no exogenous forcing).

```
make_parameters_exogenous_null(params)
make_parameters_vc_null(params)
```

After the parameters for `5-3-4`

have been specified, we
can generate the indices for the model and attach them to the parameter
list.

`make_indices(params)`

Now we can set the initial conditions of the model.

```
<- rep(NaN, max(params$X_ix))
y $L_ix] <- as.vector(L0)
y[params$M_ix] <- as.vector(M0)
y[params$G_ix] <- as.vector(G0)
y[params$Y_ix] <- as.vector(Y0)
y[params$Z_ix] <- as.vector(Z0)
y[params$Upsilon_ix] <- as.vector(Upsilon)
y[params$X_ix] <- as.vector(X0) y[params
```

Now we can pass the vector of initial conditions, `y`

, our
parameter list `params`

, and the function
`exDE::xDE_diffeqn`

to the differential equation solvers in
`deSolve::ode`

to generate a numerical trajectory. The
classes of `Xpar`

, `MYZpar`

, and `Lpar`

in `params`

will ensure that the right methods are invoked
(dispatched) to solve your model.

`<- deSolve::ode(y = y, times = 0:365, func = xDE_diffeqn, parms = params, method = "lsoda") out `

With a small amount of data wrangling made easier by the
`data.table`

package, we can plot the output.

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