Ross-Macdonald Mosquito Infection Dynamics

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 $\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} = \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} = \Lambda_t + p e^{-f q \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^{-f q \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 \left(1-G_i \right) Y_{i, t}$$

and

$$Y_{\tau, t+1} = p \left(1-G_{\tau-1} \right) Y_{\tau-1, t} + p \left(1-G_\tau \right) 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 $\cal K$ denote a matrix describing the fraction of emigrating mosquitoes that end up in every other patch, where $$\mbox{diag}\; {\cal 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 $\Omega$ denote the demographic matrix where:

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

$$\begin{array}{rl} M_{t+1} &= \Lambda_t + \Omega \cdot M_t \\ P_{t+1} &= f (M_t - P_t) + \Omega \cdot P_t \\ \hline U_{t+1} &= \Lambda_t + \Omega \cdot \left(e^{-fq\kappa_t} U_t \right) \\ Y_{1, t+1} &= \Omega \cdot \left(\left(1-e^{-fq\kappa_t}\right) U_t \right) \\ Y_{i+1, t+1} &= \Omega \cdot Y_{i,t} \left(1-G_i \right) \\ Y_{\tau, t+1} &= \Omega \cdot \left( Y_{\tau-1, t}\;\left(1-G_{\tau-1} \right) \right) + \Omega \cdot \left( Y_{\tau, t} \; \left(1-G_\tau \right) \right) \\ Z_{t+1} &= \Omega \cdot \left(Y_t \cdot \mbox{diag}(G) \right)+ \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} &= \Omega \cdot \left(\left(1-e^{-fq\kappa_t}\right) U_t \right) \\ \end{array}$$

Verification

We can verify our model by solving for steady states.

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

$$\bar P = f \bar M \cdot (1-\Omega+ \mbox{diag}\;f )^{-1}$$

$$\bar U = \Lambda \cdot (1-\Omega \cdot \mbox{diag}\; e^{-fq\kappa})^{-1}$$

$$\bar Y_1 = \Omega \cdot \left( 1- e^{-fq\kappa}\right) \cdot \bar U$$ For $1 < i < \tau,$ we get a recursive relationship:

$$\bar Y_{i+1} = \Omega \cdot \; \bar Y_i (1-G_i)$$ and for $Y_{\tau},$ we get:

$$\bar Y_{\tau} = \Omega \cdot \left( \mbox{diag} \left(1-G\right) \right) \bar Y_{\tau-1} \cdot \left(1 - \Omega \left (1-G_\tau \right) \right)^{-1}$$ and

$$\bar Z = \Omega \cdot G \cdot Y \left( 1- \Omega\right)^{-1}$$

Demo

devtools::load_all()

rm1 <- dts_setup(Xname = "trace")

dts_solve(rm1, Tmax=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") 
})

dts_plot_YZ(rm1)

compute_MYZ_equil = function(pars, Lambda, kappa, i=1){
  Omega = make_Omega(0, pars$MYZpar[[i]])
  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) 
#>          M          P          U          Y          Z 
#> 12000.0000  9391.3043  9166.7797  1835.9391   997.2812

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)
)
#>          M          P          U          Y          Z 
#> 11999.9997  9391.3040  9166.7797  1835.9391   997.2808

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

dts_solve(rm10, Tmax=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") 
})