---
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)

human_out <- data.table::CJ(day = 1:tmax, MOI = 0:(nrow(MOI_init)-1), value = NaN)
data.table::setkey(human_out, day)

while (get_tnow(mod) <= mod$global$tmax) {
  mod$human$EIR <- EIR
  step_humans(model = mod)
  human_out[day==get_tnow(mod), value := as.vector(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 = human_out[day == tmax, value])
```

Now we can plot the compartment sizes.

```{R}
ggplot(data = human_out) +
  geom_line(aes(x = day, y = value, group = MOI, color = 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}
human_final <- human_out[day == tmax, ]
human_final[, 'theoretical' := dpois(x = MOI, lambda = h/r)]
human_final[, 'empirical' := value / sum(value)]

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