Basic Mosquito Infection Dynamics in Discrete Time

In a Single Patch

We define several variables describing adult female mosquitoes. Let $M_t$ denote the total density of uninfected mosquitoes at time $t.$ We assume that mosquitoes emerge from aquatic habitats at the rate ${\mathrm{\Lambda }}_t$ and that a proportion of mosquitoes, $p$, survives each day. The total density of mosquitoes is described by an equation:

$$M_{t+1}={\mathrm{\Lambda }}_t+p{M}_t$$

We let $f$ denote the proportion of mosquitoes that blood feeds each day, and let $P_t$ denote the density of parous mosquitoes:

$$P_{t+1}=f(M_t-P_t)+p{P}_t$$

To model the infection process we need two additional parameters:

• $q$ - the proportion of blood meals that are taken on humans

• ${\kappa }_t$ - fraction of human blood meals that infect the mosquito

Let $U$ denote the density of uninfected mosquitoes:

$$U_{t+1}={\mathrm{\Lambda }}_t+p{e}^{-fq{\kappa }_t}{U}_t$$ Let $Y_{i,t}$ denote a the density of cohorts of infected mosquitoes infected $i$ days ago:

$$Y_{1,t+1}=p(1-e^{-fq{\kappa }_t})U_t$$ We let $\tau$ denote the oldest cohort of infected mosquitoes that we would wish to track. For mosquitoes infected more than 1 day ago, we let $G_i$ be the fraction that becomes infectious the next day. If we modeled the EIP with a fixed delay, then $G_i=1$ for $i<\tau$, and $G_{\tau }=1$. Mosquitoes would become infectious on day $\tau +1.$

$$Y_{i+1,t+1}=p(1-G_i)Y_{i,t}$$

and

$$Y_{\tau ,t+1}=p(1-G_{\tau -1})Y_{\tau -1,t}+p(1-G_{\tau })Y_{\tau ,t}$$

$Z_t$ denote the density of infectious mosquitoes at time $t$

$$Z_{t+1}=\sum _{i}p{G}_i{Y}_i+p{Z}_t$$

As we have defined our variables, one of them is not necessary, since

$$U=M-\sum _i{Y}_i-Z$$

The Multi-Patch Model

In a model with $n$ patches, the model is defined as before, but now we let $M_t,$ $P_t,$ $U_t$, and $Z_t$ be vectors of length $p,$ and we let $G$ be a vector of length $\tau ,$ (where $G_1=0$). We let $Y_t$ denote a $p\times \tau$ matrix, where each row is a cohort of mosquitoes, and where column represents the mosquitoes in a location. We let $Y_{i,t}$ denote the $i^{th}$ cohort:

Now, we also need to define mosquito movement. We let $\sigma$ denote a vector describing the proportion of mosquitoes that emigrates from each patch, and let $\mathcal{K}$ denote a matrix describing the fraction of emigrating mosquitoes that end up in every other patch, where $$\text{diag}\mathcal{K}=0.$$ and where the columns sum up to a number that is less than or equal to one, where the gap represents mortality occurring during dispersal.

We let $\mathrm{\Omega }$ denote the demographic matrix where:

$$\mathrm{\Omega }=\text{diag}\left(p(1-\sigma )\right)+\text{diag}\left(p\sigma \right)\cdot \mathcal{K}$$ Now the dynamics are very similar:

$$\begin{array}{rl}{M}_{t+1}& ={\mathrm{\Lambda }}_t+\mathrm{\Omega }\cdot M_t\\ P_{t+1}& =f(M_t-P_t)+\mathrm{\Omega }\cdot P_t\\ U_{t+1}& ={\mathrm{\Lambda }}_t+\mathrm{\Omega }\cdot \left(e^{-fq{\kappa }_t}{U}_t\right)\\ Y_{1,t+1}& =\mathrm{\Omega }\cdot \left((1-e^{-fq{\kappa }_t})U_t\right)\\ Y_{i+1,t+1}& =\mathrm{\Omega }\cdot Y_{i,t}(1-G_i)\\ Y_{\tau ,t+1}& =\mathrm{\Omega }\cdot \left(Y_{\tau -1,t}(1-G_{\tau -1})\right)+\mathrm{\Omega }\cdot \left(Y_{\tau ,t}(1-G_{\tau })\right)\\ Z_{t+1}& =\mathrm{\Omega }\cdot (Y_t\cdot \text{diag}(G))+\mathrm{\Omega }\cdot Z_t\end{array}$$

Implementation Notes:

In the implementation, we could choose models where the EIP is distributed over some period but where the fraction maturing after time $\tau$ is continuous. We have chosen a rotating index to track age of infection cohorts in the matrix $Y;$ using modulo arithmetic, each columns tracks a single cohort for $\tau$ days. On the $\tau +1$ day, the column is added to the one before, and then we compute $Y_{1,t}.$ The computation thus does this:

Noting that in the modulo arithmetic, $\tau +1=1$, we get:

$$\begin{array}{rl}{Y}_{\tau ,t+1}& =Y_{\tau ,t+1}+Y_{\tau +1,t}\\ Y_{1,t+1}& =\mathrm{\Omega }\cdot \left((1-e^{-fq{\kappa }_t})U_t\right)\end{array}$$

Verification

We can verify our model by solving for steady states.

$$\overline{M}=\mathrm{\Lambda }\cdot (1-\mathrm{\Omega }{)}^{-1}$$

$$\overline{P}=f\overline{M}\cdot (1-\mathrm{\Omega }+\text{diag}f{)}^{-1}$$

$$\overline{U}=\mathrm{\Lambda }\cdot (1-\mathrm{\Omega }\cdot \text{diag}{e}^{-fq\kappa }{)}^{-1}$$

$${\overline{Y}}_1=\mathrm{\Omega }\cdot (1-e^{-fq\kappa })\cdot \overline{U}$$ For $1<i<\tau ,$ we get a recursive relationship:

$${\overline{Y}}_{i+1}=\mathrm{\Omega }\cdot {\overline{Y}}_i(1-G_i)$$ and for $Y_{\tau },$ we get:

$${\overline{Y}}_{\tau }=\mathrm{\Omega }\cdot \left(\text{diag}(1-G)\right){\overline{Y}}_{\tau -1}\cdot {(1-\mathrm{\Omega }(1-G_{\tau }))}^{-1}$$ and

$$\overline{Z}=\mathrm{\Omega }\cdot G\cdot Y{(1-\mathrm{\Omega })}^{-1}$$

Demo

rm1 <- xds_setup(MYZname = "RM", Xname = "trivial", xds = 'dts')

Om <- with(rm1$MYZpar[[1]], compute_Omega_dts(pp, ssigma, mmu, calK))

dts_solve(rm1, 200) -> rm1

#dts_steady(rm1) -> rm1

with(rm1$outputs$orbits$MYZ[[1]], {
  plot(time, M, type = "l") 
  lines(time, U, col = "darkblue") 
  lines(time, Y, col = "purple") 
  lines(time, Z, col = "darkred") 
})

xds_plot_YZ(rm1)

compute_MYZ_equil = function(pars, Lambda, kappa, i=1){
  with(pars$MYZpar[[i]],{

    Mbar <-  Lambda/(1-p) 
    Pbar <- f*Mbar/(1 - p + f)
    Ubar <- Lambda/(1-p*exp(-f*q*kappa))
    
    Y1 <- Omega*(1-exp(-f*q*kappa))*Ubar 
    Yi = Y1
    Y = Yi
    for(i in 2:(max_eip-1)){
      Yi <- Omega*Yi*(1-G[i])
      Y = Y+Yi
    }
    Yn <- Omega*Yi
    Y = Y + Yn
    Zbar <- Omega*Yn/(1-p) 
    return(c(M = Mbar, P=Pbar, U=Ubar, Y=Y, Z= Zbar)) 
})}

compute_MYZ_equil(rm1, 1000, .1) 

c(
M = tail(rm1$outputs$orbits$MYZ[[1]]$M, 1), 
P = tail(rm1$outputs$orbits$MYZ[[1]]$P, 1), 
U = tail(rm1$outputs$orbits$MYZ[[1]]$U, 1),
Y = tail(rm1$outputs$orbits$MYZ[[1]]$Y, 1),
Z = tail(rm1$outputs$orbits$MYZ[[1]]$Z, 1)
)

rm10 <- dts_setup(Xname = "trivial", nPatches = 10, membership = 1:10)
rm10$Lpar[[1]]$scale = 1.5^c(1:10) 

dts_solve(rm10, 200) -> rm10

with(rm10$outputs$orbits$MYZ[[1]], {
  plot(time, M[,10], type = "l")
  for(i in 2:10)
    lines(time, M[,i]) 
#  lines(time, U, col = "darkblue") 
#  lines(time, Y, col = "purple") 
#  lines(time, Z, col = "darkred") 
})