Example
Let’s do an example. We’ll need to set up both the aquatic and mosquito components. The example will run for three years. We will make \(\lambda\) a vector giving emergence for each patch, for each day of the year, peaking in the middle of the year with 100 mosquitoes emerging per day, per patch.
We also assume that daily survival p is sinusoidal, with
a period of 365 days.
We let the EIP be a constant 5 days. We will also set \(\kappa = 0.05\) for each day.
tmax <- 365 * 3
p <- 3
lambda <- dnorm(x = 1:365, mean = 180, sd = 90)
lambda <- lambda * (100/max(lambda))
lambda <- t(replicate(p, lambda))
psurv <- (sin((1:365)/365*2*pi) + 1.01)/2 * 0.9
EIP <- 5
f <- 0.3
q <- 1
psi <- matrix(
c(
0.9, 0.05, 0.05,
0.05, 0.9, 0.05,
0.05, 0.05, 0.9
), nrow = 3, ncol = 3,
byrow = TRUE
)
par(mar = c(5,5,2,5))
plot(lambda[1, ], type = "l", col = "red", xlab = "Day", ylab = "Lambda (red)")
par(new = TRUE)
plot(psurv, type = "l", axes = F, xlab = NA, ylab = NA, col = "blue")
axis(side = 4)
mtext(side = 4, line = 3, 'Survival Probability (blue)')
Deterministic simulation
Now we can set up the models. Let’s first do a deterministic
simulation. We use make_MicroMoB() to set up the base model
object and setup_mosquito_RM() to setup the RM mosquito
model for the mosquito component. setup_aqua_trace() sets
up the aquatic (immature) component as a simple trace (forced)
model.
M <- c(100, 100, 100)
Y <- c(0, 0, 0)
Z <- c(0, 0, 0)
mod <- make_MicroMoB(tmax = tmax, p = 3)
setup_mosquito_RM(mod, stochastic = FALSE, f = f, q = q, eip = EIP, p = psurv, psi = psi, M = M, Y = Y, Z = Z)
setup_aqua_trace(model = mod, lambda = lambda, stochastic = FALSE)
det_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), patch = 1:3, value = NaN)
det_out <- det_out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(det_out, day)
# run it
while(get_tnow(mod) <= tmax) {
mod$mosquito$kappa <- rep(0.05, 3)
step_aqua(model = mod)
step_mosquitoes(model = mod)
det_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
det_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
det_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
mod$global$tnow <- mod$global$tnow + 1L
}
ggplot(det_out) +
geom_line(aes(x = day, y = value, color = state)) +
facet_grid(patch ~ .)
Stochastic simulation
Now we run 10 stochastic simulations, using the same parameters.
sto_out <- mclapply(X = 1:10, FUN = function(runid) {
mod <- make_MicroMoB(tmax = tmax, p = 3)
setup_mosquito_RM(mod, stochastic = TRUE, f = f, q = q, eip = EIP, p = psurv, psi = psi, M = M, Y = Y, Z = Z)
setup_aqua_trace(model = mod, lambda = lambda, stochastic = TRUE)
out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), patch = 1:3, value = NaN)
out <- out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(out, day)
while(get_tnow(mod) <= tmax) {
mod$mosquito$kappa <- rep(0.05, 3)
step_aqua(model = mod)
step_mosquitoes(model = mod)
out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
mod$global$tnow <- mod$global$tnow + 1L
}
out[, 'run' := as.integer(runid)]
return(out)
})
sto_out <- data.table::rbindlist(sto_out)
ggplot(sto_out) +
geom_line(aes(x = day, y = value, color = state, group = interaction(run, state)), alpha = 0.35) +
facet_grid(patch ~ .)