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

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.

```
<- 0.025
h <- 1/200
r <- 0.55
b
<- -log(1 - h) / b
EIR
<- 1
n <- 1e3
tmax <- matrix(data = c(1e5, rep(0, 1e2)), nrow = 101, ncol = n)
MOI_init
<- make_MicroMoB(tmax = tmax, p = 1)
mod 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)
<- data.table::CJ(day = 1:tmax, MOI = 0:(nrow(MOI_init)-1), value = NaN)
human_out ::setkey(human_out, day)
data.table
while (get_tnow(mod) <= mod$global$tmax) {
$human$EIR <- EIR
modstep_humans(model = mod)
==get_tnow(mod), value := as.vector(mod$human$MOI)]
human_out[day$global$tnow <- mod$global$tnow + 1L
mod }
```

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.

```
weighted.mean(x = 0:100, w = human_out[day == tmax, value])
#> [1] 5.05693
```

Now we can plot the compartment sizes.

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

```
<- human_out[day == tmax, ]
human_final 'theoretical' := dpois(x = MOI, lambda = h/r)]
human_final[, 'empirical' := value / sum(value)]
human_final[,
ggplot(human_final, aes(MOI, empirical)) +
geom_bar(stat = 'identity') +
geom_line(aes(x = MOI, y = theoretical), color = "blue")
```