library(exDE)
library(MASS)
library(expm)
library(deSolve)
library(data.table)
library(ggplot2)This is the null model of aquatic mosquito dynamics; there are no
endogenous dynamics and the model is simply specified by
Lambda, the rate at which adult mosquitoes spontaneously
emerge from aquatic habitats.
Below we show an example, taken from the package tests, of using the trace-based aquatic model to keep an adult mosquito model at equilibrium when using unequal numbers of aquatic habitats per patch.
We first calculate equilibrium values following the Ross-Macdonald vignette. We use \(\mathcal{N}\) and \(\mathcal{U}\) to describe how aquatic
habitats are dispersed amongst patches, and how mosquitoes in each patch
disperse eggs to each habitat, respectively. Please note because we have
unequal numbers of aquatic habitats and patches, we use
MASS::ginv to compute the generalized inverse of \(\mathcal{N}\) to get \(\alpha\) required at equilibrium.
nPatches <- 3
nHabitats <- 4
f <- 0.3
q <- 0.9
g <- 1/20
sigma <- 1/10
tau <- 11
nu <- 1/2
eggsPerBatch <- 30
calN <- matrix(0, nPatches, nHabitats)
calN[1,1] <- 1
calN[2,2] <- 1
calN[3,3] <- 1
calN[3,4] <- 1
calU <- matrix(0, nHabitats, nPatches)
calU[1,1] <- 1
calU[2,2] <- 1
calU[3:4,3] <- 0.5
calK <- matrix(0, nPatches, nPatches)
calK[1, 2:3] <- c(0.2, 0.8)
calK[2, c(1,3)] <- c(0.5, 0.5)
calK[3, 1:2] <- c(0.7, 0.3)
calK <- t(calK)
Omega <- make_Omega(g, sigma, calK, nPatches)
Upsilon <- expm::expm(-Omega * tau)
kappa <- c(0.1, 0.075, 0.025)
Lambda <- c(5, 10, 8)
# equilibrium solutions
Omega_inv <- solve(Omega)
M_eq <- as.vector(Omega_inv %*% Lambda)
G_eq <- as.vector(solve(diag(nu+f, nPatches) + Omega) %*% diag(f, nPatches) %*% M_eq)
Y_eq <- as.vector(solve(diag(f*q*kappa) + Omega) %*% diag(f*q*kappa) %*% M_eq)
Z_eq <- as.vector(Omega_inv %*% Upsilon %*% diag(f*q*kappa) %*% (M_eq - Y_eq))
# the "Lambda" for the dLdt model
alpha <- as.vector(ginv(calN) %*% Lambda)Now we set up the model. Please see the Ross-Macdonald vignette for details on the
adult model. We use alpha as the Lambda
parameter of our forced emergence model, in
exDE::make_parameters_L_trace. Again, because we are not
running the full transmission model, we must use
exDE::MosquitoBehavior to pass the equilibrium values of
those bionomic parameters to exDE::xDE_diffeqn_mosy.
params <- list(
nPatches = nPatches,
nHabitats = nHabitats,
calU = calU,
calN = calN
)
params <- list2env(params)
# ODE
make_parameters_MYZ_RM_ode(pars = params, g = g, sigma = sigma, calK = calK, tau = tau, f = f, q = q, nu = nu, eggsPerBatch = eggsPerBatch, M0 = rep(0, nPatches), G0 = rep(0, nPatches), Y0 = rep(0, nPatches), Z0 = rep(0, nPatches))
make_parameters_L_trace(pars = params, Lambda = alpha)
make_indices(params)
y0 <- rep(0, max(params$Upsilon_ix))
y0[params$M_ix] <- M_eq
y0[params$G_ix] <- G_eq
y0[params$Y_ix] <- Y_eq
y0[params$Z_ix] <- Z_eq
y0[params$Upsilon_ix] <- as.vector(Upsilon)
MosyBehavior <- MosquitoBehavior(0, y, params)
out <- deSolve::ode(y = y0, times = 0:50, func = xDE_diffeqn_mosy, parms = params, method = 'lsoda', kappa = kappa, MosyBehavior = MosyBehavior)
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)
out <- out[, -c(params$Upsilon_ix+1)]
out <- as.data.table(out)
out <- melt(out, id.vars = 'time')
out[, c("Component", "Patch") := tstrsplit(variable, '_', fixed = TRUE)]
out[, variable := NULL]
ggplot(data = out, mapping = aes(x = time, y = value, color = Patch)) +
geom_line() +
facet_wrap(. ~ Component, scales = 'free') +
theme_bw()