Ross-Macdonald transmission model

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:

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:

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:

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:

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

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

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

Finally plug this into the first equation to solve λ:

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 pEIP instead of e−gEIP.

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

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

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.

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.

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