---
title: "Ross-Macdonald transmission model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ross-Macdonald transmission 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)
```

In MicroMoB, the goal is to get various models fill in a set of component interfaces;
if you want to write a new aquatic model, you fill in the interface and then plug
into the existing framework, and let the simulation handle the rest! Here we
present the most simple case of a Ross-Macdonald mosquito model linked to a
trace derived (exogeneously forced) model of emergence, and a SIS (Susceptible-Infectious-Susceptible)
human model. This is in fact quite similar to the "classic" Ross-Macdonald
(though best written by Aron and May) model.

## Endemic equilibrium solution

First we will write out the closest approximating ODE model and derive an endemic
equilibrium to check that our discrete time models are approximately correct.

Let the human component be described by:

\begin{equation}
\dot{S} = -b fqZ (1/N) S + rI \\
\dot{I} = b fqZ (1/N) S - rI
\end{equation}

here, $fqZ$ is the EIR, and $N$ is the total human population size. If we know
the values of $S$ and $I$ we can solve for the density of infectious mosquitoes
required to give those values at endemic equilibrium:

\begin{equation}
  Z = \frac{rIN}{bfqS}
\end{equation}

Now let us look at the system of ODEs describing the mosquitoes. Actually they are
delay differential equations (DDE) to incorporate the delay between infection and
infectiousness:

\begin{equation}
  \dot{M} = \lambda - gM \\
  \dot{Y} = fq\kappa (M-Y) - gY \\
  \dot{Z} = fq\kappa_{t-EIP} (M_{t-EIP} - Y_{t-EIP}) e^{-gEIP} - gZ
\end{equation}

We are only interested in equilibrium solutions, where all derivatives are zero
and the state variables are the same for all time, so we can ignore the delays:

\begin{equation}
  0 = \lambda - gM \\
  0 = fq\kappa (M-Y) - gY \\
  0 = fq\kappa (M - Y) e^{-gEIP} - gZ
\end{equation}

Plug our solution for $Z$ into the last equation to get an expression for $M-Y$ in
terms of known quantities:

\begin{equation}
M-Y = \frac{gZ}{fq\kappa e^{-gEIP}}
\end{equation}

Now plug $M-Y$ into the second equation to get an expression for $Y$, after some
simplification we get:

\begin{equation}
  Y = \frac{Z}{e^{-gEIP}}
\end{equation}

To get $M$, add $(M-Y) + Y$ and simplify to get:

\begin{equation}
M = \frac{Z(g + fq\kappa)}{fq\kappa e^{-gEIP}}
\end{equation}

Finally plug this into the first equation to solve $\lambda$:

\begin{equation}
  \lambda = g \left( \frac{Z(g + fq\kappa)}{fq\kappa e^{-gEIP}} \right)
\end{equation}

We can set up some model parameters and calculate numeric solutions for those values. Note that the actual equilibrium calculations we use in the code are different by a factor of $p$ in the terms of $M$ because of how events are ordered in the discrete time model, and that we use $p^{EIP}$ instead of $e^{-gEIP}$.

```{R}
# mosquito parameters
f <- 0.3
q <- 1
eip <- 14

lifespan <- 20
g <- 1/lifespan
p <- 1- g

# human parameters
b <- 0.55
c <- 0.15
r <- 1/200

S <- 1e3
I <- 300
N <- S + I

# transmission parameters
kappa <- (I/N)*c

# equilibrium solutions
Z <- (r*I*N) / (b*f*q*S)
Y <- Z / (p^eip)
M <- (Z*(g + (f*q*p*kappa))) / (f*q*p*kappa*(p^eip))
lambda <- g*M
```

## MicroMoB simulation

Now let's set up a simulation with only a single patch and only a single human stratum.

```{R}
patches <- 1
nstrata <- 1
tmax <- 365 * 2
theta <- diag(nstrata)
psi <- diag(patches)
```

We first run a deterministic simulation. The first thing we need to do is create
the base model object using `make_MicroMoB()` and initialize all the components with specific models, which we do
below. `setup_humans_SIS()` sets up the human component to use a SIS model of infection.

We do not use `setup_alternative_trace()` or `setup_visitor_trace()` to set up the other blood hosts
and visitor components in this simple model, so we use `compute_bloodmeal_simple()` to compute
$\kappa$ and $EIR$ instead of the more complex `compute_bloodmeal()`.

`setup_mosquito_RM()` sets up the adult mosquito component using the Ross-Macdonald model, and
`setup_aqua_trace()` sets up the aquatic mosquito component using a trace.

```{R}
mod <- make_MicroMoB(tmax = tmax, p = patches)
setup_humans_SIS(mod, stochastic = FALSE, theta = theta, H = N, X = I, b = b, c = c, r = r)
setup_aqua_trace(mod, stochastic = FALSE, lambda = lambda)
setup_mosquito_RM(mod, stochastic = FALSE, f = f, q = q, eip = eip, p = p, psi = psi, M = M, Y = Y, Z = Z)
```

Now we can make matrices to store output and run the model. Note that we call `compute_bloodmeal()`
_before_ any other updating functions are called. This is crucial for properly updating
models in **Micro-MoB**.

```{R}
mosy_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), value = NaN, species = "mosquito")
mosy_out <- mosy_out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(mosy_out, day)

human_out <- data.table::CJ(day = 1:tmax, state = c('S', 'I'), value = NaN, species = "human")
human_out <- human_out[c('S', 'I'), on="state"]
data.table::setkey(human_out, day)

while (get_tnow(mod) <= tmax) {
  
  compute_bloodmeal_simple(model = mod)
  step_aqua(model = mod)
  step_mosquitoes(model = mod)
  step_humans(model = mod)
  
  mosy_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
  mosy_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
  mosy_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
  
  human_out[day == get_tnow(mod) & state == 'S', value := mod$human$H - mod$human$X]
  human_out[day == get_tnow(mod) & state == 'I', value := mod$human$X]

  mod$global$tnow <- mod$global$tnow + 1L
}

det_out <- rbind(mosy_out, human_out)
```

Now we draw 10 trajectories from the stochastic simulation, and plot output. We plot
the cloud of stochastic trajectories as faint lines and the deterministic solution as
solid lines.

```{R}
sto_out <- mclapply(X = 1:10, FUN = function(runid) {
  
  mod <- make_MicroMoB(tmax = tmax, p = patches)
  setup_humans_SIS(mod, stochastic = TRUE, theta = theta, H = N, X = I, b = b, c = c, r = r)
  setup_aqua_trace(mod, stochastic = TRUE, lambda = lambda)
  setup_mosquito_RM(mod, stochastic = TRUE, f = f, q = q, eip = eip, p = p, psi = psi, M = M, Y = Y, Z = Z)
  
  mosy_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), value = NaN, species = "mosquito")
  mosy_out <- mosy_out[c('M', 'Y', 'Z'), on="state"]
  data.table::setkey(mosy_out, day)
  
  human_out <- data.table::CJ(day = 1:tmax, state = c('S', 'I'), value = NaN, species = "human")
  human_out <- human_out[c('S', 'I'), on="state"]
  data.table::setkey(human_out, day)
    
  while (get_tnow(mod) <= tmax) {
    
    compute_bloodmeal_simple(model = mod)
    step_aqua(model = mod)
    step_mosquitoes(model = mod)
    step_humans(model = mod)
    
    mosy_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
    mosy_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
    mosy_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
    
    human_out[day == get_tnow(mod) & state == 'S', value := mod$human$H - mod$human$X]
    human_out[day == get_tnow(mod) & state == 'I', value := mod$human$X]
  
    mod$global$tnow <- mod$global$tnow + 1L
  }
  
  out <- rbind(mosy_out, human_out)
  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 = run), alpha = 0.3) +
    geom_line(data = det_out, aes(x = day, y = value, color = state)) +
    facet_wrap(species ~ state, scales = "free")
```