A Ross-Macdonald Adult Mosquito Model

Source: vignettes/adult_RM.Rmd

The Ross-Macdonald adult mosquito model fulfills the generic interface of the adult mosquito component.

Here, we use a Ross-Macdonald model based on a model first published in 1982 by Joan Aron and Robert May1. It includes state variables for total mosquito density MM, infected mosquito density YY, and infectious mosquito density ZZ. Because of the interest in parity, the model has been extended to include a variable that tracks parous mosquitoes, PP.

Differential Equations

Delay Differential Equations

These equations are naturally implemented by ramp.xds::dMYZdt.RM_dde, but they can also be implemented in a closely related set of odes using ramp.xds::dMYZdt.RM_ode (see below).

Ṁ=Λ−Ω⋅M \dot{M} = \Lambda - \Omega\cdot M Ṗ=diag(f)⋅(M−P)−Ω⋅P \dot{P} = \mbox{diag}(f) \cdot (M-P) - \Omega \cdot P

Ẏ=diag(fqκ)⋅(M−Y)−Ω⋅Y \dot{Y} = \mbox{diag}(fq\kappa) \cdot (M-Y) - \Omega \cdot Y

Ż=e−Ωτ⋅diag(fqκt−τ)⋅(Mt−τ−Yt−τ)−Ω⋅Z \dot{Z} = e^{-\Omega \tau} \cdot \mbox{diag}(fq\kappa_{t-\tau}) \cdot (M_{t-\tau}-Y_{t-\tau}) - \Omega \cdot Z

Recall that the mosquito demography matrix describing mortality and dispersal is given by:

$$ \Omega = \mbox{diag(g)} + \left(I- {\cal K}\right) \cdot \mbox{diag}(\sigma) $$

Ordinary Differential Equations

In the following, the equation are solved by ramp.xds::dMYZdt.RM_ode.

The system of ODEs is the same as above except for the equation giving the rate of change in infectious mosquito density, which becomes:

Ż=e−Ωτ⋅diag(fqκ)⋅(M−Y)−Ω⋅Z \dot{Z} = e^{-\Omega \tau} \cdot \mbox{diag}(fq\kappa) \cdot (M-Y) - \Omega \cdot Z The resulting set of equations is similar in spirit to the simple model presented in Smith & McKenzie (2004)2. in that mortality and dispersal over the EIP is accounted for, but the time lag is not. While transient dynamics of the ODE model will not equal the DDE model, they have the same equilibrium values, and so for numerical work requiring finding equilibrium points, the faster ODE model can be safely substituted.

Equilibrium solutions

There are two logical ways to begin solving the non-trivial equilibrium. The first assumes Λ\Lambda is known, which implies good knowledge of mosquito ecology. The second assumes ZZ is known, which implies knowledge of the biting rate on the human population. We show both below.

Starting with Λ\Lambda

Given Λ\Lambda we can solve:

M=Ω−1⋅Λ M = \Omega^{-1} \cdot \Lambda Then given MM we set Ẏ\dot{Y} to zero and factor out YY to get:

Y=(diag(fqκ)+Ω)−1⋅diag(fqκ)⋅M Y = (\mbox{diag}(fq\kappa) + \Omega)^{-1} \cdot \mbox{diag}(fq\kappa) \cdot M We set Ż\dot{Z} to zero to get:

Z=Ω−1⋅e−Ωτ⋅diag(fqκ)⋅(M−Y) Z = \Omega^{-1} \cdot e^{-\Omega \tau} \cdot \mbox{diag}(fq\kappa) \cdot (M-Y)

Because the dynamics of PP are independent of the infection dynamics, we can solve it given MM as:

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

Starting with ZZ

It is more common that we start from an estimate of ZZ, perhaps derived from an estimated EIR (entomological inoculation rate). Given ZZ, we can calculate the other state variables and Λ\Lambda. For numerical implementation, note that (e−Ωτ)−1=eΩτ(e^{-\Omega\tau})^{-1} = e^{\Omega\tau}.

M−Y=diag(1/fqκ)⋅(e−Ωτ)−1⋅Ω⋅Z M-Y = \mbox{diag}(1/fq\kappa) \cdot (e^{-\Omega\tau})^{-1} \cdot \Omega \cdot Z

Y=Ω−1⋅diag(fqκ)⋅(M−Y) Y = \Omega^{-1} \cdot \mbox{diag}(fq\kappa) \cdot (M-Y)

M=(M−Y)+Y M = (M - Y) + Y

Λ=Ω⋅M \Lambda = \Omega \cdot M We can use the same equation for PP as above.

Example

Here we show an example of starting and solving a model at equilibrium. Please note that this only runs this adult mosquito model and that most users should read our fully worked example to run a full simulation.

Ross-Macdonald

dydt=aκ(1−y)−gydzdt=e−gτqκτ[1−yτ]−gz \begin{array}{rl} \frac{\textstyle{dy}}{\textstyle{dt}} & = a \kappa \left(1 -y \right) - g y \\ \frac{\textstyle{dz}}{\textstyle{dt}} & = e^{-g\tau} q \kappa_\tau \left[1 - y_\tau \right] - g z \end{array}

The long way

Here we set up some parameters for a simulation with 3 patches.

HPop = rep(1, 3)
nPatches <- 3
f <- rep(0.3, nPatches)
q <- rep(0.9, nPatches)
g <- rep(1/20, nPatches)
sigma <- rep(1/10, nPatches)
mu <- rep(0, nPatches)
eip <- 12
nu <- 1/2
eggsPerBatch <- 30
MYZo = list(f=f,q=q,g=g,sigma=sigma,mu=mu,eip=eip,nu=nu,eggsPerBatch=eggsPerBatch)
calK <- matrix(0, nPatches, nPatches)
calK[1, 2:3] <- c(0.2, 0.8)
calK[2, c(1,3)] <- c(0.5, 0.5)
calK[3, 1:2] <- c(0.7, 0.3)
calK <- t(calK)

Omega <- compute_Omega_xde(g, sigma, mu, calK)
Upsilon <- expm::expm(-Omega * eip)

Now we set up the parameter environment with the correct class using ramp.xds::make_parameters_MYZ_RM_xde, noting that we will be solving as an ode.

Now we set the values of κ\kappa and Λ\Lambda and solve for the equilibrium values.

kappa <- c(0.1, 0.075, 0.025)
Xo = list(kappa=kappa)
Lambda <- c(5, 10, 8)
Lo = list(Lambda=Lambda)
Omega_inv <- solve(Omega)
M_eq <- as.vector(Omega_inv %*% Lambda)
P_eq <- as.vector(solve(diag(f, nPatches) + Omega) %*% diag(f, nPatches) %*% M_eq)
Y_eq <- as.vector(solve(diag(f*q*kappa) + Omega) %*% diag(f*q*kappa) %*% M_eq)
Z_eq <- as.vector(Omega_inv %*% Upsilon %*% diag(f*q*kappa) %*% (M_eq - Y_eq))
MYZo$M=M_eq
MYZo$P=P_eq
MYZo$Y=Y_eq
MYZo$Z=Z_eq

We use ramp.xds::make_inits_MYZ_RM_xde to store the initial values. These equations have been implemented to compute Υ\Upsilon dynamically, so we attach Upsilon as initial values:

params <- make_xds_template("dde", "mosy", nPatches, 1:3, 1:3)
params <- make_MYZpar("macdonald", params, 1, MYZo)  
params <- make_MYZinits(params, 1, MYZo)
params <- make_Xpar("trivial", params, 1, Xo) 
params <- make_Xinits(params, HPop, 1)
params <- make_Lpar("trivial", params, 1, Lo)
params <- make_Linits(params, 1, Lo)
params <- setup_Hpar_static(params, 1)

We set the indices with ramp.xds::make_indices.

params = make_indices(params)
params <- change_calK(calK, params, 1)
params$Lambda[[1]] = Lambda
params$kappa[[1]] = kappa

Then we can set up the initial conditions vector and use deSolve::ode to solve the model. Normally these values would be computed within ramp.xds::xDE_diffeqn. Here, we set up a local version:

y0 = get_MYZinits(params, 1) 
y0 = as.vector(unlist(y0))
params <- MBionomics(0,y0,params, 1)

dMYZdt_local = func=function(t, y, pars) {
  list(dMYZdt(t, y, pars, 1))
}

out <- deSolve::dede(y = y0, times = 0:50, dMYZdt_local, parms=params, 
                    method = 'lsoda') 
out1 <- out

The output is plotted below. The flat lines shown here is a verification that the steady state solutions that we computed above match the steady states computed by solving the equations:

out = out[, 1:13]
colnames(out)[params$MYZpar$M_ix+1] <- paste0('M_', 1:params$nPatches)
colnames(out)[params$MYZpar$P_ix+1] <- paste0('P_', 1:params$nPatches)
colnames(out)[params$MYZpar$Y_ix+1] <- paste0('Y_', 1:params$nPatches)
colnames(out)[params$MYZpar$Z_ix+1] <- paste0('Z_', 1:params$nPatches)

out <- as.data.table(out)
out <- melt(out, id.vars = 'time')
out[, c("Component", "Patch") := tstrsplit(variable, '_', fixed = TRUE)]
out[, variable := NULL]

ggplot(data = out, mapping = aes(x = time, y = value, color = Patch)) +
  geom_line() +
  facet_wrap(. ~ Component, scales = 'free') +
  theme_bw()

Setup Utilities

In the vignette above, we set up a function to solve the differential equation. We hope this helps the end user to understand how ramp.xds works under the hood, but the point of ramp.xds is to lower the costs of building, analyzing, and using models. The functionality in ramp.xds can handle this case – we can set up and solve the same model using built-in setup utilities. Learning to use these utilities makes it very easy to set up other models without having to learn some internals.

To set up a model with the parameters above, we make three list with the named parameters and their values. We also attach to the list the initial values we want to use, if applicable. For the Ross-Macdonald adult mosquito model, we attach the parameter values:

Each one of the dynamical components has a configurable trivial algorithm that computes no derivatives, but passes its output as a parameter (see human-trivial.R. To configure Xpar, we attach the values of kappa to a list:

Xo = list(kappa =  c(0.1, 0.075, 0.025))

Similarly, we configure the trivial algorithm for aquatic mosquitoes (see aquatic-trivial.R).

Lo = list(Lambda = c(5, 10, 8))

To set up the model, we call xds_setup with

  • nPatches is set to 3

  • MYZname is set to “macdonald” to run the Ross-Macdonald model for adult mosquitoes; to pass our options, we write MYZopts = MYZo; and finally, the dispersal matrix calK is passed as calK=calK

  • Xname is set to “trivial” to run the trivial module for human infection dynamics

  • Lname is set to “trivial” to run the trivial module aquatic mosquitoes

Otherwise, setup takes care of all the internals:

xds_setup(MYZname = "macdonald", Xname = "trivial", Lname = "trivial",  
    nPatches=3, calK=calK, membership = c(1:3), 
    residence = c(1:3), HPop = HPop,
    MYZopts = MYZo, Xopts = Xo, Lopts = Lo) -> MYZeg

Now, we can solve the equations using xds_solve and compare the output to what we got above. If they are identical, the two objects should be identical, so can simply add the absolute value of their differences:

xds_solve(MYZeg, Tmax=50, dt=1) -> MYZeg 
out2 <- MYZeg$outputs$orbits$deout
sum(abs(out1-out2))==0 
#> [1] TRUE
rbind(tail(out2,1)[1 + 1:12],
c(M_eq, P_eq, Y_eq, Z_eq))
#>         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#> [1,] 157.868 123.4518 178.6802 140.3517 98.87622 155.0578 48.09598 27.56429
#>          [,9]    [,10]    [,11]    [,12]
#> [1,] 41.03339 24.96933 13.31699 25.75651

  1. The population dynamics of malaria. In The Population Dynamics of Infectious Diseases: Theory and Applications, R. M. Anderson, ed. (Springer US), pp. 139–179. online↩︎

  2. Smith, D.L., Ellis McKenzie, F. Statics and dynamics of malaria infection in Anopheles mosquitoes. Malar J 3, 13 (2004). online↩︎