inst/doc/BH_aqua.R

In MicroMoB: Discrete Time Simulation of Mosquito-Borne Pathogen Transmission

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(MicroMoB)
library(ggplot2)
library(data.table)
library(parallel)

## -----------------------------------------------------------------------------
patches <- 1
tmax <- 1e2

M <- 120
p <- 0.9
lambda <- M*(1-p)

nu <- 25
f <- 0.3
eggs <- nu * f * M

# static pars
molt <-  0.1
surv <- 0.9

# solve L
L <- lambda * ((1/molt) - 1) + eggs
K <- - (lambda * L) / (lambda - L*molt*surv)

## -----------------------------------------------------------------------------
# deterministic run
mod <- make_MicroMoB(tmax = tmax, p = patches)
setup_aqua_BH(model = mod, stochastic = FALSE, molt = molt, surv = surv, K = K, L = L)
setup_mosquito_RM(model = mod, stochastic = FALSE, f = f, q = 0.9, eip = 10, p = p, psi = diag(1), nu = nu, M = M, Y = 0, Z = 0)

out_det <- data.table::CJ(day = 1:tmax, state = c('L', 'A', 'M'), value = NaN)
out_det <- out_det[c('L', 'A', 'M'), on="state"]
data.table::setkey(out_det, day)

while (get_tnow(mod) <= tmax) {
  
  step_aqua(model = mod)
  step_mosquitoes(model = mod)
  
  out_det[day == get_tnow(mod) & state == 'L', value := mod$aqua$L]
  out_det[day == get_tnow(mod) & state == 'A', value := mod$aqua$A]
  out_det[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]

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

## -----------------------------------------------------------------------------
# stochastic runs
out_sto <- mclapply(X = 1:10, FUN = function(runid) {
  
  mod <- make_MicroMoB(tmax = tmax, p = patches)
  setup_aqua_BH(model = mod, stochastic = TRUE, molt = molt, surv = surv, K = K, L = L)
  setup_mosquito_RM(model = mod, stochastic = TRUE, f = f, q = 0.9, eip = 10, p = p, psi = diag(1), nu = nu, M = M, Y = 0, Z = 0)
  
  out <- data.table::CJ(day = 1:tmax, state = c('L', 'A', 'M'), value = NaN)
  out <- out[c('L', 'A', 'M'), on="state"]
  data.table::setkey(out, day)
  
  while (get_tnow(mod) <= tmax) {
    
    step_aqua(model = mod)
    step_mosquitoes(model = mod)
    
    out[day == get_tnow(mod) & state == 'L', value := mod$aqua$L]
    out[day == get_tnow(mod) & state == 'A', value := mod$aqua$A]
    out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
  
    mod$global$tnow <- mod$global$tnow + 1L
  }
  
  out[, 'run' := as.integer(runid)] 
  return(out)
})

out_sto <- data.table::rbindlist(out_sto)

## -----------------------------------------------------------------------------
ggplot(out_sto) +
    geom_line(aes(x = day, y = value, color = state, group = run), alpha = 0.35) +
    geom_line(data = out_det, aes(x = day, y = value, color = state)) +
    facet_wrap(. ~ state, scales = "free")