Basic Mosquito Infection Dynamics in Discrete Time

Discrete-Time Dynamics

Source: vignettes/RM-dts-Mosquito.Rmd

In this vignette, we present a model for the dynamics of malaria infection in adult mosquito populations. We call it a Ross-Macdonald model because it makes the same basic assumptions as the model that Macdonald analyzed in 1952. While that model was developed using differential equations, the equations we describe herein are difference equations.

We begin by discrete-time model in a simple population, with one patch, and then we present the equations for the multi-patch model.

In a Single Patch

We define several variables describing adult female mosquitoes. Let MtM_t denote the total density of uninfected mosquitoes at time t.t. We assume that mosquitoes emerge from aquatic habitats at the rate Λt\Lambda_t and that a proportion of mosquitoes, pp, survives each day. The total density of mosquitoes is described by an equation:

Mt+1=Λt+pMtM_{t+1} = \Lambda_t + p M_t

We let ff denote the proportion of mosquitoes that blood feeds each day, and let PtP_t denote the density of parous mosquitoes:

Pt+1=f(Mt−Pt)+pPtP_{t+1} = f(M_t - P_t) + p P_t

To model the infection process we need two additional parameters:

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

  • κt\kappa_t - fraction of human blood meals that infect the mosquito

Let UU denote the density of uninfected mosquitoes:

Ut+1=Λt+pe−fqκtUtU_{t+1} = \Lambda_t + p e^{-f q \kappa_t} U_t Let Yi,tY_{i,t} denote a the density of cohorts of infected mosquitoes infected ii days ago:

Y1,t+1=p(1−e−fqκt)UtY_{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 GiG_i be the fraction that becomes infectious the next day. If we modeled the EIP with a fixed delay, then Gi=1G_i = 1 for i<τi < \tau, and Gτ=1G_\tau = 1. Mosquitoes would become infectious on day τ+1.\tau+1.

Yi+1,t+1=p(1−Gi)Yi,tY_{i+1, t+1} = p \left(1-G_i \right) Y_{i, t}

and

Yτ,t+1=p(1−Gτ−1)Yτ−1,t+p(1−Gτ)Yτ,tY_{\tau, t+1} = p \left(1-G_{\tau-1} \right) Y_{\tau-1, t} + p \left(1-G_\tau \right) Y_{\tau, t}

ZtZ_t denote the density of infectious mosquitoes at time tt

Zt+1=∑ipGiYi+pZtZ_{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−∑iYi−ZU = M- \sum_i Y_i - Z

The Multi-Patch Model

In a model with nn patches, the model is defined as before, but now we let Mt,M_t,Pt,P_t,UtU_t, and ZtZ_t be vectors of length p,p, and we let GG be a vector of length τ,\tau, (where G1=0G_1 = 0). We let YtY_t denote a p×τp \times \tau matrix, where each row is a cohort of mosquitoes, and where column represents the mosquitoes in a location. We let Yi,tY_{i,t} denote the ithi^{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:

Mt+1=Λt+Ω⋅MtPt+1=f(Mt−Pt)+Ω⋅PtUt+1=Λt+Ω⋅(e−fqκtUt)Y1,t+1=Ω⋅((1−e−fqκt)Ut)Yi+1,t+1=Ω⋅Yi,t(1−Gi)Yτ,t+1=Ω⋅(Yτ−1,t(1−Gτ−1))+Ω⋅(Yτ,t(1−Gτ))Zt+1=Ω⋅(Yt⋅diag(G))+Ω⋅Zt \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;Y; using modulo arithmetic, each columns tracks a single cohort for τ\tau days. On the τ+1\tau +1 day, the column is added to the one before, and then we compute Y1,t.Y_{1,t}. The computation thus does this:

Mt+1=Λt+Ω⋅MtPt+1=f(Mt−Pt)+Ω⋅PtUt+1=Λt+Ω⋅(e−fqκtUt)Yt+1=Ω⋅(Yi,t⋅diag(1−G))Zt+1=Ω⋅(Yt⋅diagG)+Ω⋅Zt\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_{t+1} &= \Omega \cdot \left(Y_{i, t} \cdot \mbox{diag} \left(1-G \right) \right) \\ Z_{t+1} &= \Omega \cdot \left( Y_t \cdot \mbox{diag}\; G \right)+ \Omega \cdot Z_t \\ \end{array}

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

Yτ,t+1=Yτ,t+1+Yτ+1,tY1,t+1=Ω⋅((1−e−fqκt)Ut) \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.

M‾=Λ⋅(1−Ω)−1\bar M = \Lambda \cdot (1-\Omega)^{-1}

P‾=fM‾⋅(1−Ω+diagf)−1\bar P = f \bar M \cdot (1-\Omega+ \mbox{diag}\;f )^{-1}

U‾=Λ⋅(1−Ω⋅diage−fqκ)−1\bar U = \Lambda \cdot (1-\Omega \cdot \mbox{diag}\; e^{-fq\kappa})^{-1}

Y‾1=Ω⋅(1−e−fqκ)⋅U‾\bar Y_1 = \Omega \cdot \left( 1- e^{-fq\kappa}\right) \cdot \bar U For 1<i<τ,1 < i < \tau, we get a recursive relationship:

Y‾i+1=Ω⋅Y‾i(1−Gi)\bar Y_{i+1} = \Omega \cdot \; \bar Y_i (1-G_i) and for Yτ,Y_{\tau}, we get:

Y‾τ=Ω⋅(diag(1−G))Y‾τ−1⋅(1−Ω(1−Gτ))−1\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

Z‾=Ω⋅G⋅Y(1−Ω)−1\bar Z = \Omega \cdot G \cdot Y \left( 1- \Omega\right)^{-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") 
})