knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(MicroMoB) library(ggplot2) library(data.table) library(parallel)
The simple behavioral state mosquito model has two behavioral states which mosquitoes can exist in: blood feeding ($B$) and oviposition ($Q$). When mosquitoes are in $B$ they will attempt to blood feed until they are successful, at which point they transition to $Q$ and attempt to oviposit an egg batch. Upon emergence, mosquitoes are primed for blood feeding and are in $B$. They transition between these states until they die, which occurs according to the state dependent probabilities $p_{B}$ and $p_{Q}$ (these may also vary by location and time). The model does not consider male mosquitoes.
The model also considers infection. Uninfected (susceptible) mosquitoes $M$ may become infected if they are in $B$, successfully take a blood meal, and are infected (with probability $\kappa$). They then transition to the infected class $Y$, in behavioral state $Q$. The extrinsic incubation period (EIP) may vary with time, and they advance until they become infectious (if they survive), where they remain until death. Both dynamics operate simultaneously.
The deterministic behavioral state model has the following form:
\begin{equation}
\left[
\begin{array}{cc}
B_{t+1} \
Q_{t+1} \
\end{array}
\right]
=
\left[
\begin{array}{ccc}
(1 - \psi_b) \Psi_{b b} & \psi_q \Psi_{q b} \
\psi_b \Psi_{b q} & (1 - \psi_q) \Psi_{q q} \
\end{array}
\right]
\left[
\begin{array}{cc}
p_b B_{t} \
p_q Q_{t} \
\end{array}
\right]
+
\left[
\begin{array}{c}
\Lambda_{t} \
0 \
\end{array}
\right]
\end{equation}
The state is a column vector $\left[\begin{array}{cc} B \ Q \ \end{array}\right]$. We assume that there are $p$ locations where mosquitoes go to seek blood hosts, so that the first $p$ elements correspond to the number of mosquitoes in the $B$ state at those places. There are $l$ locations where mosquitoes go to oviposit (aquatic habitats), so the last $l$ elements in the vector are mosquitoes in the $Q$ state. There is no requirement that the set of points where mosquitoes blood feed and oviposit be distinct, although they may be.
The infection states are similar to the Ross-Macdonald model, see vignette("RM_mosquito")
for more details.
The parameters in the state updating equation are:
The stochastic model has similar updating dynamics to the deterministic implementation, except that all survival and success probabilities are used in binomial draws and movement is drawn from a multinomial distribution.
We assume that $p = l = 1$ and that the total mosquito density $M = B + Q$ is known, and that we want to solve for the emergence rate $\Lambda$ such that the system is at equilibrium. Rewriting the equations when we substitute $Q = M - B$ and $B = M - Q$ we solve the state variables as:
\begin{equation} Q = \frac{Mp_{B}\Psi_{B}}{p_{B}\Psi_{B} - p_{Q}(1-\Psi_{Q}) + 1} \ B = \frac{M-Mp_{Q}(1-\Psi_{Q})}{p_{B}\Psi_{B} - p_{Q}(1-\Psi_{Q}) + 1} \end{equation}
Then the first equation can simply be rearranged to yield:
\begin{equation} \Lambda = B - p_{B}(1-\Psi_{B})B - p_{Q}\Psi_{Q}Q \end{equation}
And now the model with 1 point of each type can be set up at equilibrium. We will
use the Beverton-Holt model of aquatic ecology demonstrated in vignette("BH_aqua")
,
which will be parameterized to provide the correct equilibrium $\Lambda$.
p <- l <- 1 tmax <- 1e2 M <- 120 pB <- 0.8 pQ <- 0.95 PsiB <- 0.5 PsiQ <- 0.85 B <- (M - (M*pQ*(1-PsiQ))) / ((pB*PsiB) - (pQ*(1-PsiQ)) + 1) Q <- (M*pB*PsiB) / ((pB*PsiB) - (pQ*(1-PsiQ)) + 1) lambda <- B - (pB*(1-PsiB)*B) - (pQ*PsiQ*Q) nu <- 25 eggs <- nu * PsiQ * Q # static pars molt <- 0.1 surv <- 0.9 # solve L L <- lambda * ((1/molt) - 1) + eggs K <- - (lambda * L) / (lambda - L*molt*surv)
Let's set up the model. We use make_MicroMoB()
to set up the base model object, and
setup_aqua_BH()
for the Beverton-Holt aquatic model with our chosen parameters.
setup_mosquito_BQ()
will set up a behavioral state model of adult mosquito dynamics.
We run a deterministic simulation and store output in a matrix. Note that we calculate
the f
and q
parameters to achieve the correct PsiB
probability; normally these
would be updated dynamically during the bloodmeal but we are running a mosquito-only
simulation so we set these deterministically.
# deterministic run mod <- make_MicroMoB(tmax = tmax, p = p, l = l) setup_aqua_BH(model = mod, stochastic = FALSE, molt = molt, surv = surv, K = K, L = L) setup_mosquito_BQ(model = mod, stochastic = FALSE, eip = 5, pB = pB, pQ = pQ, psiQ = PsiQ, Psi_bb = matrix(1), Psi_bq = matrix(1), Psi_qb = matrix(1), Psi_qq = matrix(1), nu = nu, M = c(B, Q), Y = matrix(0, nrow = 2, ncol = 6)) out_det <- data.table::CJ(day = 1:tmax, state = c('L', 'A', 'B', 'Q'), value = NaN) out_det <- out_det[c('L', 'A', 'B', 'Q'), on="state"] data.table::setkey(out_det, day) mod$mosquito$q <- 0.3 mod$mosquito$f <- log(1 - PsiB) / -0.3 while (get_tnow(mod) <= tmax) { step_aqua(model = mod) step_mosquitoes(model = mod) out_det[day == get_tnow(mod) & state == 'L', value := mod$aqua$L] out_det[day == get_tnow(mod) & state == 'A', value := mod$aqua$A] out_det[day == get_tnow(mod) & state == 'B', value := mod$mosquito$M[1]] out_det[day == get_tnow(mod) & state == 'Q', value := mod$mosquito$M[2]] mod$global$tnow <- mod$global$tnow + 1L }
Now we run the same model, but using the option stochastic = TRUE
for our dynamics,
and draw 10 trajectories.
# stochastic runs out_sto <- mclapply(X = 1:10, FUN = function(runid) { mod <- make_MicroMoB(tmax = tmax, p = p, l = l) setup_aqua_BH(model = mod, stochastic = TRUE, molt = molt, surv = surv, K = K, L = L) setup_mosquito_BQ(model = mod, stochastic = TRUE, eip = 5, pB = pB, pQ = pQ, psiQ = PsiQ, Psi_bb = matrix(1), Psi_bq = matrix(1), Psi_qb = matrix(1), Psi_qq = matrix(1), nu = nu, M = c(B, Q), Y = matrix(0, nrow = 2, ncol = 6)) out <- data.table::CJ(day = 1:tmax, state = c('L', 'A', 'B', 'Q'), value = NaN) out <- out[c('L', 'A', 'B', 'Q'), on="state"] data.table::setkey(out, day) mod$mosquito$q <- 0.3 mod$mosquito$f <- log(1 - PsiB) / -0.3 while (get_tnow(mod) <= tmax) { step_aqua(model = mod) step_mosquitoes(model = mod) out[day == get_tnow(mod) & state == 'L', value := mod$aqua$L] out[day == get_tnow(mod) & state == 'A', value := mod$aqua$A] out[day == get_tnow(mod) & state == 'B', value := mod$mosquito$M[1]] out[day == get_tnow(mod) & state == 'Q', value := mod$mosquito$M[2]] mod$global$tnow <- mod$global$tnow + 1L } out[, 'run' := as.integer(runid)] return(out) })
Now we process the output and plot the results. Deterministic solutions are solid lines and each stochastic trajectory is a faint line.
out_sto <- data.table::rbindlist(out_sto) ggplot(data = out_sto) + geom_line(aes(x = day, y = value, color = state, group = run), alpha = 0.35) + geom_line(data = out_det, aes(x = day, y = value, color = state)) + facet_wrap(. ~ state, scales = "free")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.