Introduction
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.
Diagram
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.
Parameters
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.
params = make_parameters_xde()
params$nVectors = 1
params$nHosts = 1
params$nHabitats = 5
params$nPatches = 3
params$nStrata = 4Aquatic Habitat Membership Matrix
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}\]
Egg Dispersal Matrix
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}\]
Aquatic Mosquito Parameters
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.
L0 <- rep(1, params$nHabitats)
psi <- rep(1/8, params$nHabitats)
phi <- rep(1/8, params$nHabitats)
theta <- c(1/10, 1/20, 1/40, 1/100, 1/10)
params = make_parameters_L_basic(params, psi = psi, phi = phi, theta = theta)
params = make_inits_L_basic(params, L0=L0)Adult Mosquito Parameters
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.
g <- 1/12
sigma <- 1/12/2
f <- 1/3
q <- 0.9
nu <- c(1/3,1/3,0)
eggsPerBatch <- 30
eip <- 12
calK <- t(matrix(
c(c(0, .6, .3),
c(.4, 0, .7),
c(.6, .4, 0)), 3, 3))
M0 <- rep(100, params$nPatches)
P0 <- rep(10, params$nPatches)
Y0 <- rep(1, params$nPatches)
Z0 <- rep(0, params$nPatches)
Omega <- make_Omega(g, sigma, calK, params$nPatches)
Upsilon <- expm::expm(-Omega*eip)
params = make_parameters_MYZ_RM(pars = params, g = g,
sigma = sigma, calK=calK,
eip=eip, f=f, q=q, nu=nu,
eggsPerBatch=eggsPerBatch,
solve_as = "ode")
params = make_inits_MYZ_RM_ode(params, M0=M0, P0=P0,
Y0=Y0, Z0=Z0)Mixing
We use the static demographic model, which assumes a constant population size (constant \(H\)).
H <- matrix(c(10,90, 100, 900), 4, 1)
params = setup_Hpar_static(params, 1, HPop=H) In this model, we define four population strata. We can describe their residency with a vector describing residence:
residence = c(2,2,3,3)
searchWtsH = c(1,1,1,1)
params = setup_BloodFeeding(params, 1, 1, residence=residence, searchWts=searchWtsH)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}\]
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}\]
TaR <- t(matrix(
c(c(0.01,0.01,0.001,0.001),
c(.95,.92,.04,.02),
c(.04,.02,.959,.929)), 4, 3
))
params$BFpar$TaR[[1]][[1]] <- TaRWe 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.
I0 <- as.vector(0.2*H)
r <- 1/200
b <- 0.55
c <- c(0.1, .02, .1, .02)
params = make_parameters_X_SIS(pars = params, b = b, c = c, r = r)
params = make_inits_X_SIS(params, S0 = H-I0, I0=I0)Simulation, the Long Way
Initial Conditions
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.
params = make_indices(params)Now we can set the initial conditions of the model.
y0 = get_inits(params)
params <- EggLaying(0, y0, params)Numerical Solution
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.
out <- deSolve::ode(y = y0, times = 0:365,
func = xDE_diffeqn, parms = params, method = "lsoda")
out1 <- outPlot Output
With a small amount of data wrangling made easier by the
data.table package, we can plot the output.
colnames(out)[params$ix$L[[1]]$L_ix+1] <- paste0('L_', 1:params$nHabitats)
colnames(out)[params$ix$MYZ[[1]]$M_ix+1] <- paste0('M_', 1:params$nPatches)
colnames(out)[params$ix$MYZ[[1]]$P_ix+1] <- paste0('P_', 1:params$nPatches)
colnames(out)[params$ix$MYZ[[1]]$Y_ix+1] <- paste0('Y_', 1:params$nPatches)
colnames(out)[params$ix$MYZ[[1]]$Z_ix+1] <- paste0('Z_', 1:params$nPatches)
colnames(out)[params$ix$X[[1]]$X_ix+1] <- paste0('X_', 1:params$nStrata)
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()
Using Setup
We create lists with all our parameters values:
MYZo = list(
solve_as = "ode",
g = 1/12, sigma = 1/12/2,
f = 1/3, q=0.9, nu=c(1/3,1/3,0),
eggsPerBatch = 30,
eip = 12,
M0 = 100, P0 = 10, Y0 = 1, Z0 = 0
)
EIPo <- list(eip=12)
xde_setup("mod534", "RM", "SIS", "basic",
nPatches = 3, HPop=c(10, 90, 100, 900),
membership=c(1,1,1,2,2), EIPopts = EIPo,
MYZopts=MYZo, calK=calK, Xopts=Xo,
residence=c(2,2,3,3), searchB=searchWtsH,
TimeSpent=TaR, searchQ = c(7,2,1,8,2), Lopts = Lo) -> mod534We solve and take the differences to check:
mod534 <- xde_solve(mod534,Tmax=365,dt=1)
mod534$outputs$orbits$deout -> out2Interestingly, the differences are small:
approx_equal(tail(out2, 1), tail(out1,1), tol = 1e-5)
#> time L1 L2 L3 L4 L5 MYZ1 MYZ2 MYZ3 MYZ4 MYZ5 MYZ6 MYZ7 MYZ8
#> [366,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> MYZ9 MYZ10 MYZ11 MYZ12 X1 X2 X3 X4 X5 X6 X7 X8
#> [366,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE