---
title: "Simple behavioral state mosquito model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple behavioral state mosquito model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
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.

## Deterministic model

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:

  * $\psi_b$: probability of successful blood feeding (vector of length $p$); this parameter is computed from $f, q$ (themselves calculated during the bloodmeal algorithm) as $1-e^{-fq}$.
  * $\psi_q$: probability of successful oviposition (vector of length $l$).
  * $\Psi_{b b}$: transition probability matrix for movement among blood feeding haunts. It has dimension $p\times p$, and has _columns_ that sum to 1 (note state vectors are on the right).
  * $\Psi_{q b}$: transition probability matrix for movement from aquatic habitats to blood feeding haunts. It has dimension $p\times l$.
  * $\Psi_{b q}$: transition probability matrix for movement from blood feeding haunts to aquatic habitats. It has dimension $l\times p$.
  * $\Psi_{q q}$: transition probability matrix for movement among aquatic habitats. It has dimension $l\times l$. 
  * $p_{B}$: daily survival probability for blood feeding mosquitoes.
  * $p_{Q}$: daily survival probability for ovipositing mosquitoes.
  
## Stochastic model

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.

## Simulation

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$.

```{R}
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.

```{R}
# 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.

```{R}
# 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.

```{R}
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")
```
