Hybrid MoI (Multiplicity of Infection) Human Model • ramp.xds

This is a hybrid model which tracks the mean multiplicity of infection (superinfection) in two compartments. The first, $m 1$ is all infections, and the second $m 2$ are apparent (patent) infections. Therefore $m 2$ is “nested” within $m 1$ . It is a “hybrid” model in the sense of Nåsell (1985).

Differential Equations

The equations are as follows:

$\dot{m_{1}} = h - r_{1} m_{1}$ $\dot{m_{2}} = h - r_{2} m_{2}$ Where $h = b E I R$ , is the force of infection. Prevalence can be calculated from these MoI values by:

$x_{1} = 1 - e^{- m_{1}}$ $x_{2} = 1 - e^{- m_{2}}$ The net infectious probability to mosquitoes is therefore given by:

$x = c_{2} x_{2} + c_{1} (x_{1} - x_{2})$

Where $c_{1}$ is the infectiousness of inapparent infections, and $c_{2}$ is the infectiousness of patent infections.

Equilibrium solutions

One way to proceed is assume that $m_{2}$ is known, as it models the MoI of patent (observable) infections. Then we have:

$h = r_{2} / m_{2}$ $m_{1} = h / r_{1}$ We can use this to calculate the net infectious probability, and then $κ = β^{⊤} \cdot x$ , allowing the equilibrium solutions of this model to feed into the other components.

Example

library(ramp.xds)
library(deSolve)
library(viridisLite)

Here we run a simple example with 3 population strata at equilibrium. We use ramp.xds::make_parameters_X_hMoI_xde to set up parameters. Please note that this only runs the human population component and that most users should read our fully worked example to run a full simulation.

We use the null (constant) model of human demography ( $H$ constant for all time). While $H$ does not appear in the equations above, it must be specified as the model of bloodfeeding ( $β$ ) relies on $H$ to compute consistent values.

The Long Way

Here, we build a model step-by-step.

nStrata <- 3
H <- c(100, 500, 250)
residence = rep(1,3) 
params <- make_xds_object_template("ode", "cohort", 1, 1, residence)

b <- 0.55
c1 <- 0.05
c2 <- 0.25
r1 <- 1/250
r2 <- 1/50

m2 <- rlnorm(3, .2, .5) 
h <- r2*m2
m1 <- h/r1
Xo = list(b = b, c1 = c1, c2 = c2, r1 = r1, r2 = r2, m1=m1, m2=m2)

eir <- h/b

params = setup_Xpar("hMoI", params, 1, Xo)
params = setup_Xinits(params, H, 1, Xo)

F_season = function(t){0*t+1}
F_trend = function(t){0*t+1}
F_age = function(a){0*a+1}

params$EIRpar = list() 
params$EIRpar$eir <- as.vector(eir)
params$EIRpar$scale <- 1 
params$EIRpar$F_season <- F_season
params$EIRpar$F_trend <- F_trend
params$EIRpar$F_age <- F_age
params$EIRpar$season_par <- list() 
params$EIRpar$trend_par<- list() 
params$EIRpar$age_par<- list()

params = make_indices(params)

y0 <- get_inits(params)
y0
#> $L
#> list()
#> 
#> $MYZ
#> NULL
#> 
#> $X
#> $X$m1
#> [1] 3.032587 6.938575 1.805443
#> 
#> $X$m2
#> [1] 0.6065175 1.3877150 0.3610886

params <- xds_solve_cohort(params) 
out1 <- params$outputs$orbits

XH <- out1$XH[[1]] 
age <- out1$age
clrs = turbo(5)
plot(age, XH$m1[,1], col = clrs[1], ylim = range(XH$m1), type = "l")
lines(age, XH$m1[,2], col = clrs[2])
lines(age, XH$m1[,3], col =clrs[5])
  
lines(age, XH$m2[,1], col = clrs[1], lty = 2)
lines(age, XH$m2[,2], col = clrs[2], lty=2)
lines(age, XH$m2[,3], col =clrs[5], lty=2)

Using Setup

xds_setup_cohort(eir, Xname="hMoI", HPop=H, Xopts = Xo) ->test_hMoI

xds_solve_cohort(test_hMoI)-> test_hMoI 
XH2 <- params$outputs$orbits$XH[[1]]
sum((XH$m1-XH$m1)^2)
#> [1] 0