Hybrid MoI (Multiplicity of Infection) Human Model

Source: vignettes/human_hmoi.Rmd

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

Differential Equations

The equations are as follows:

m1̇=h−r1m1 \dot{m_{1}} = h - r_{1}m_{1} m2̇=h−r2m2 \dot{m_{2}} = h - r_{2}m_{2} Where h=bEIRh = b EIR, is the force of infection. Prevalence can be calculated from these MoI values by:

x1=1−e−m1 x_{1} = 1-e^{-m_{1}} x2=1−e−m2 x_{2} = 1-e^{-m_{2}} The net infectious probability to mosquitoes is therefore given by:

x=c2x2+c1(x1−x2) x = c_{2}x_{2} + c_{1}(x_{1}-x_{2})

Where c1c_{1} is the infectiousness of inapparent infections, and c2c_{2} is the infectiousness of patent infections.

Equilibrium solutions

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

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

Example

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 (HH constant for all time). While HH does not appear in the equations above, it must be specified as the model of bloodfeeding (β\beta) relies on HH 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_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 = make_Xpar("hMoI", params, 1, Xo)
params = make_Xinits(params, H, 1, Xo) 
F_season = function(t){1}
F_trend = function(t){1}
F_age = function(a){1}
params$EIRpar = list() 
params$EIRpar$eir <- as.vector(eir)
params$EIRpar$F_season <- F_season
params$EIRpar$F_trend <- F_trend
params$EIRpar$F_age <- F_age
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
out1 <- xds_solve_cohort(params)$outputs$cohort
clrs = turbo(5)
with(out1,{
  plot(age, m1[,1], col = clrs[1], ylim = range(m1), type = "l")
  lines(age, m1[,2], col = clrs[2])
  lines(age, m1[,3], col =clrs[5])
  
  lines(age, m2[,1], col = clrs[1], lty = 2)
  lines(age, m2[,2], col = clrs[2], lty=2)
  lines(age, 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)$outputs$cohort -> out2
sum((out2$XH[[1]]$m1-out1$XH[[1]]$m1)^2)
#> [1] 0