---
title: "Queueing model for human superinfection"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Queueing model for human superinfection}
  %\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)
```

One of the earliest elaborations on Ross's early SIS style models of malaria infection
in human populations was a queueing model which was developed when it became apparent
that _superinfection_, when an individual host is simultaneously infected with multiple
distinct parasite broods, was an important concept in malaria epidemiology.

The basic model was presented as an $(M/M/\infty)$ queueing process, with state space
$X_{0}, X_{1}, \dots$, where the subscripts denote the number of parasite broods
an individual is infected with (such that $0$ is an uninfected person). The number of broods
infecting a person is known as the multiplicity of infection (MOI).

The force of infection is denoted $h$ and the recovery rate from compartment $m$ is $\rho_{m}$.
Then the deterministic dynamics can be expressed as:

\begin{equation}
\dot{X}_{0} = -h X_{0} + \rho_{1} X_{1} \\
\dot{X}_{m} = -(h + \rho_{m})X_{m} + h X_{m-1} + \rho_{m+1} X_{m+1}
\end{equation}

The prevalence of disease (also known as the parasite rate in malaria) is $X = 1 - \frac{X_{0}}{H}$,
where $H$ is the total human population.

In our formulation, we let the recovery rate have the following form:

\begin{equation}
\rho_{m} = r m^{\sigma}
\end{equation}

When $\sigma = 1$, parasite broods clear independently. By setting $\sigma > 1$ clearance
rates become faster and competition is simulated; $\sigma < 1$ means slower clearance
due to facilitation between parasites. When broods clear independently the
distribution of MOI in the population is Poisson with mean $h/r$.

In **Micro-MoB** the `MOI` vector is allowed to grow as needed to accommodate arbitrarily
large values of MOI, but may become computationally expensive in such cases.

## Simulation

Let's check that we recover approximately the correct distribution over MOI.
First we set up and run a simulation for 1000 days. We use `make_MicroMoB()`
to set up the base model object and `setup_humans_MOI()` with our chosen parameters
to set up the multiplicity of infection human model.

```{R}
h <- 0.025
r <- 1/200
b <- 0.55

EIR <- -log(1 - h) / b

n <- 1
tmax <- 1e3
MOI_init <- matrix(data = c(1e5, rep(0, 1e2)), nrow = 101, ncol = n)

mod <- make_MicroMoB(tmax = tmax, p = 1)
setup_humans_MOI(model = mod, stochastic = TRUE, theta = matrix(1, nrow = n, ncol = 1), H = colSums(MOI_init), MOI = MOI_init, r = r, b = b)

MOI_out <- matrix(data = NaN, nrow = nrow(MOI_init), ncol = tmax + 1)
MOI_out[, 1L] <- MOI_init

while (mod$global$tnow <= mod$global$tmax) {
  mod$human$EIR <- EIR
  step_humans(model = mod)
  MOI_out[, mod$global$tnow + 1L] <- mod$human$MOI
  mod$global$tnow <- mod$global$tnow + 1L
}
```

With these parameter values we expect a mean MOI of about 5. We see that the simulation
converges to this equilibrium distribution after some time.

```{R}
weighted.mean(x = 0:100, w = MOI_out[, tmax])
```

Now we can plot the compartment sizes.

```{R}
# plot output
MOI_out_dt <- as.data.table(t(MOI_out))
MOI_out_dt <- suppressWarnings(melt(MOI_out_dt))
setnames(MOI_out_dt, new = c("MOI", "Count"))
levels(MOI_out_dt$MOI) <- 0:100
MOI_out_dt[, "Day" := 0:tmax, by = MOI]

ggplot(data = MOI_out_dt) +
  geom_line(aes(x = Day, y = Count, group = MOI, color = as.numeric(MOI)))
```

Another check for Poisson-ness of the MOI distribution is to plot the compartment sizes at the last time point.
The theoretical distribution is plotted as a blue line.

```{R}
MOI_final <- MOI_out_dt[Day == tmax, ]
MOI_final[, "prop" := Count / sum(Count)]
MOI_final[, "MOI" := 0:100]
MOI_final[, "theoretical" := dpois(x = MOI, lambda = h/r)]

ggplot(MOI_final, aes(MOI, prop)) +
    geom_bar(stat = 'identity') +
    geom_line(aes(x = MOI, y = theoretical), color = "blue")
```
